/*  * This file is part of the UCB release of Plan 9. It is subject to the license * terms in the LICENSE file found in the top-level directory of this * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No * part of the UCB release of Plan 9, including this file, may be copied, * modified, propagated, or distributed except according to the terms contained * in the LICENSE file. */ #include <u.h>#include <libc.h>#include "map.h" /* For Albers formulas see Deetz and Adams "Elements of Map Projection", *//* USGS Special Publication No. 68, GPO 1921 */ static double r0sq, r1sq, d2, n, den, sinb1, sinb2;static struct coord plat1, plat2;static southpole; static double num(double s){	if(d2==0)		return(1);	s = d2*s*s;	return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));} /* Albers projection for a spheroid, good only when N pole is fixed */ static intXspalbers(struct place *place, double *x, double *y){	double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);	double t = n*place->wlon.l;	*y = r*cos(t);	*x = -r*sin(t);	if(!southpole)		*y = -*y;	else		*x = -*x;	return(1);} /* lat1, lat2: std parallels; e2: squared eccentricity */ static proj albinit(double lat1, double lat2, double e2){	double r1;	double t;	for(;;) {		if(lat1 < -90)			lat1 = -180 - lat1;		if(lat2 > 90)			lat2 = 180 - lat2;		if(lat1 <= lat2)			break;		t = lat1; lat1 = lat2; lat2 = t;	}	if(lat2-lat1 < 1) {		if(lat1 > 89)			return(azequalarea());		return(0);	}	if(fabs(lat2+lat1) < 1)		return(cylequalarea(lat1));	d2 = e2;	den = num(1.);	deg2rad(lat1,&plat1);	deg2rad(lat2,&plat2);	sinb1 = plat1.s*num(plat1.s)/den;	sinb2 = plat2.s*num(plat2.s)/den;	n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) -	    plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) /	    (2*(1-e2)*den*(sinb2-sinb1));	r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s));	r1sq = r1*r1;	r0sq = r1sq + 2*(1-e2)*den*sinb1/n;	southpole = lat1<0 && plat2.c>plat1.c;	return(Xspalbers);} projsp_albers(double lat1, double lat2){	return(albinit(lat1,lat2,EC2));} projalbers(double lat1, double lat2){	return(albinit(lat1,lat2,0.));} static double scale = 1;static double twist = 0; voidalbscale(double x, double y, double lat, double lon){	struct place place;	double alat, alon, x1,y1;	scale = 1;	twist = 0;	invalb(x,y,&alat,&alon);	twist = lon - alon;	deg2rad(lat,&place.nlat);	deg2rad(lon,&place.wlon);	Xspalbers(&place,&x1,&y1);	scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));} voidinvalb(double x, double y, double *lat, double *lon){	int i;	double sinb_den, sinp;	x *= scale;	y *= scale;	*lon = atan2(-x,fabs(y))/(RAD*n) + twist;	sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));	sinp = sinb_den;	for(i=0; i<5; i++)		sinp = sinb_den/num(sinp);	*lat = asin(sinp)/RAD;}