1*59cc4ca5SDavid du Colombier #include <u.h>
2*59cc4ca5SDavid du Colombier #include <libc.h>
33e12c5d1SDavid du Colombier #include "map.h"
43e12c5d1SDavid du Colombier
53e12c5d1SDavid du Colombier /* For Albers formulas see Deetz and Adams "Elements of Map Projection", */
63e12c5d1SDavid du Colombier /* USGS Special Publication No. 68, GPO 1921 */
73e12c5d1SDavid du Colombier
8219b2ee8SDavid du Colombier static double r0sq, r1sq, d2, n, den, sinb1, sinb2;
93e12c5d1SDavid du Colombier static struct coord plat1, plat2;
103e12c5d1SDavid du Colombier static southpole;
113e12c5d1SDavid du Colombier
num(double s)12219b2ee8SDavid du Colombier static double num(double s)
133e12c5d1SDavid du Colombier {
143e12c5d1SDavid du Colombier if(d2==0)
153e12c5d1SDavid du Colombier return(1);
163e12c5d1SDavid du Colombier s = d2*s*s;
173e12c5d1SDavid du Colombier return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));
183e12c5d1SDavid du Colombier }
193e12c5d1SDavid du Colombier
203e12c5d1SDavid du Colombier /* Albers projection for a spheroid, good only when N pole is fixed */
213e12c5d1SDavid du Colombier
223e12c5d1SDavid du Colombier static int
Xspalbers(struct place * place,double * x,double * y)23219b2ee8SDavid du Colombier Xspalbers(struct place *place, double *x, double *y)
243e12c5d1SDavid du Colombier {
25219b2ee8SDavid du Colombier double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);
26219b2ee8SDavid du Colombier double t = n*place->wlon.l;
273e12c5d1SDavid du Colombier *y = r*cos(t);
283e12c5d1SDavid du Colombier *x = -r*sin(t);
293e12c5d1SDavid du Colombier if(!southpole)
303e12c5d1SDavid du Colombier *y = -*y;
313e12c5d1SDavid du Colombier else
323e12c5d1SDavid du Colombier *x = -*x;
333e12c5d1SDavid du Colombier return(1);
343e12c5d1SDavid du Colombier }
353e12c5d1SDavid du Colombier
363e12c5d1SDavid du Colombier /* lat1, lat2: std parallels; e2: squared eccentricity */
373e12c5d1SDavid du Colombier
albinit(double lat1,double lat2,double e2)38219b2ee8SDavid du Colombier static proj albinit(double lat1, double lat2, double e2)
393e12c5d1SDavid du Colombier {
407dd7cddfSDavid du Colombier double r1;
41219b2ee8SDavid du Colombier double t;
423e12c5d1SDavid du Colombier for(;;) {
433e12c5d1SDavid du Colombier if(lat1 < -90)
443e12c5d1SDavid du Colombier lat1 = -180 - lat1;
453e12c5d1SDavid du Colombier if(lat2 > 90)
463e12c5d1SDavid du Colombier lat2 = 180 - lat2;
473e12c5d1SDavid du Colombier if(lat1 <= lat2)
483e12c5d1SDavid du Colombier break;
493e12c5d1SDavid du Colombier t = lat1; lat1 = lat2; lat2 = t;
503e12c5d1SDavid du Colombier }
513e12c5d1SDavid du Colombier if(lat2-lat1 < 1) {
523e12c5d1SDavid du Colombier if(lat1 > 89)
533e12c5d1SDavid du Colombier return(azequalarea());
543e12c5d1SDavid du Colombier return(0);
553e12c5d1SDavid du Colombier }
563e12c5d1SDavid du Colombier if(fabs(lat2+lat1) < 1)
573e12c5d1SDavid du Colombier return(cylequalarea(lat1));
583e12c5d1SDavid du Colombier d2 = e2;
593e12c5d1SDavid du Colombier den = num(1.);
603e12c5d1SDavid du Colombier deg2rad(lat1,&plat1);
613e12c5d1SDavid du Colombier deg2rad(lat2,&plat2);
623e12c5d1SDavid du Colombier sinb1 = plat1.s*num(plat1.s)/den;
633e12c5d1SDavid du Colombier sinb2 = plat2.s*num(plat2.s)/den;
643e12c5d1SDavid du Colombier n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) -
653e12c5d1SDavid du Colombier plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) /
663e12c5d1SDavid du Colombier (2*(1-e2)*den*(sinb2-sinb1));
673e12c5d1SDavid du Colombier r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s));
683e12c5d1SDavid du Colombier r1sq = r1*r1;
693e12c5d1SDavid du Colombier r0sq = r1sq + 2*(1-e2)*den*sinb1/n;
703e12c5d1SDavid du Colombier southpole = lat1<0 && plat2.c>plat1.c;
713e12c5d1SDavid du Colombier return(Xspalbers);
723e12c5d1SDavid du Colombier }
733e12c5d1SDavid du Colombier
743e12c5d1SDavid du Colombier proj
sp_albers(double lat1,double lat2)75219b2ee8SDavid du Colombier sp_albers(double lat1, double lat2)
763e12c5d1SDavid du Colombier {
773e12c5d1SDavid du Colombier return(albinit(lat1,lat2,EC2));
783e12c5d1SDavid du Colombier }
793e12c5d1SDavid du Colombier
803e12c5d1SDavid du Colombier proj
albers(double lat1,double lat2)81219b2ee8SDavid du Colombier albers(double lat1, double lat2)
823e12c5d1SDavid du Colombier {
833e12c5d1SDavid du Colombier return(albinit(lat1,lat2,0.));
843e12c5d1SDavid du Colombier }
853e12c5d1SDavid du Colombier
86219b2ee8SDavid du Colombier static double scale = 1;
87219b2ee8SDavid du Colombier static double twist = 0;
883e12c5d1SDavid du Colombier
893e12c5d1SDavid du Colombier void
albscale(double x,double y,double lat,double lon)90219b2ee8SDavid du Colombier albscale(double x, double y, double lat, double lon)
913e12c5d1SDavid du Colombier {
923e12c5d1SDavid du Colombier struct place place;
93219b2ee8SDavid du Colombier double alat, alon, x1,y1;
943e12c5d1SDavid du Colombier scale = 1;
953e12c5d1SDavid du Colombier twist = 0;
963e12c5d1SDavid du Colombier invalb(x,y,&alat,&alon);
973e12c5d1SDavid du Colombier twist = lon - alon;
983e12c5d1SDavid du Colombier deg2rad(lat,&place.nlat);
993e12c5d1SDavid du Colombier deg2rad(lon,&place.wlon);
1003e12c5d1SDavid du Colombier Xspalbers(&place,&x1,&y1);
1013e12c5d1SDavid du Colombier scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));
1023e12c5d1SDavid du Colombier }
1033e12c5d1SDavid du Colombier
1043e12c5d1SDavid du Colombier void
invalb(double x,double y,double * lat,double * lon)105219b2ee8SDavid du Colombier invalb(double x, double y, double *lat, double *lon)
1063e12c5d1SDavid du Colombier {
1073e12c5d1SDavid du Colombier int i;
108219b2ee8SDavid du Colombier double sinb_den, sinp;
1093e12c5d1SDavid du Colombier x *= scale;
1103e12c5d1SDavid du Colombier y *= scale;
1113e12c5d1SDavid du Colombier *lon = atan2(-x,fabs(y))/(RAD*n) + twist;
1123e12c5d1SDavid du Colombier sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));
1133e12c5d1SDavid du Colombier sinp = sinb_den;
1143e12c5d1SDavid du Colombier for(i=0; i<5; i++)
1153e12c5d1SDavid du Colombier sinp = sinb_den/num(sinp);
1163e12c5d1SDavid du Colombier *lat = asin(sinp)/RAD;
1173e12c5d1SDavid du Colombier }
118