1 28994509 2004-04-21 devnull #include <u.h>
2 28994509 2004-04-21 devnull #include <libc.h>
3 28994509 2004-04-21 devnull #include "map.h"
5 28994509 2004-04-21 devnull /* For Albers formulas see Deetz and Adams "Elements of Map Projection", */
6 28994509 2004-04-21 devnull /* USGS Special Publication No. 68, GPO 1921 */
8 28994509 2004-04-21 devnull static double r0sq, r1sq, d2, n, den, sinb1, sinb2;
9 28994509 2004-04-21 devnull static struct coord plat1, plat2;
10 28994509 2004-04-21 devnull static int southpole;
12 28994509 2004-04-21 devnull static double num(double s)
14 28994509 2004-04-21 devnull if(d2==0)
15 28994509 2004-04-21 devnull return(1);
16 28994509 2004-04-21 devnull s = d2*s*s;
17 28994509 2004-04-21 devnull return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));
20 28994509 2004-04-21 devnull /* Albers projection for a spheroid, good only when N pole is fixed */
22 28994509 2004-04-21 devnull static int
23 28994509 2004-04-21 devnull Xspalbers(struct place *place, double *x, double *y)
25 28994509 2004-04-21 devnull double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);
26 28994509 2004-04-21 devnull double t = n*place->wlon.l;
27 28994509 2004-04-21 devnull *y = r*cos(t);
28 28994509 2004-04-21 devnull *x = -r*sin(t);
29 28994509 2004-04-21 devnull if(!southpole)
30 28994509 2004-04-21 devnull *y = -*y;
32 28994509 2004-04-21 devnull *x = -*x;
33 28994509 2004-04-21 devnull return(1);
36 28994509 2004-04-21 devnull /* lat1, lat2: std parallels; e2: squared eccentricity */
38 28994509 2004-04-21 devnull static proj albinit(double lat1, double lat2, double e2)
40 28994509 2004-04-21 devnull double r1;
41 28994509 2004-04-21 devnull double t;
42 28994509 2004-04-21 devnull for(;;) {
43 28994509 2004-04-21 devnull if(lat1 < -90)
44 28994509 2004-04-21 devnull lat1 = -180 - lat1;
45 28994509 2004-04-21 devnull if(lat2 > 90)
46 28994509 2004-04-21 devnull lat2 = 180 - lat2;
47 28994509 2004-04-21 devnull if(lat1 <= lat2)
49 28994509 2004-04-21 devnull t = lat1; lat1 = lat2; lat2 = t;
51 28994509 2004-04-21 devnull if(lat2-lat1 < 1) {
52 28994509 2004-04-21 devnull if(lat1 > 89)
53 28994509 2004-04-21 devnull return(azequalarea());
54 28994509 2004-04-21 devnull return(0);
56 28994509 2004-04-21 devnull if(fabs(lat2+lat1) < 1)
57 28994509 2004-04-21 devnull return(cylequalarea(lat1));
59 28994509 2004-04-21 devnull den = num(1.);
60 28994509 2004-04-21 devnull deg2rad(lat1,&plat1);
61 28994509 2004-04-21 devnull deg2rad(lat2,&plat2);
62 28994509 2004-04-21 devnull sinb1 = plat1.s*num(plat1.s)/den;
63 28994509 2004-04-21 devnull sinb2 = plat2.s*num(plat2.s)/den;
64 28994509 2004-04-21 devnull n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) -
65 28994509 2004-04-21 devnull plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) /
66 28994509 2004-04-21 devnull (2*(1-e2)*den*(sinb2-sinb1));
67 28994509 2004-04-21 devnull r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s));
68 28994509 2004-04-21 devnull r1sq = r1*r1;
69 28994509 2004-04-21 devnull r0sq = r1sq + 2*(1-e2)*den*sinb1/n;
70 28994509 2004-04-21 devnull southpole = lat1<0 && plat2.c>plat1.c;
71 28994509 2004-04-21 devnull return(Xspalbers);
75 28994509 2004-04-21 devnull sp_albers(double lat1, double lat2)
77 28994509 2004-04-21 devnull return(albinit(lat1,lat2,EC2));
81 28994509 2004-04-21 devnull albers(double lat1, double lat2)
83 28994509 2004-04-21 devnull return(albinit(lat1,lat2,0.));
86 28994509 2004-04-21 devnull static double scale = 1;
87 28994509 2004-04-21 devnull static double twist = 0;
90 28994509 2004-04-21 devnull albscale(double x, double y, double lat, double lon)
92 28994509 2004-04-21 devnull struct place place;
93 28994509 2004-04-21 devnull double alat, alon, x1,y1;
94 28994509 2004-04-21 devnull scale = 1;
95 28994509 2004-04-21 devnull twist = 0;
96 28994509 2004-04-21 devnull invalb(x,y,&alat,&alon);
97 28994509 2004-04-21 devnull twist = lon - alon;
98 28994509 2004-04-21 devnull deg2rad(lat,&place.nlat);
99 28994509 2004-04-21 devnull deg2rad(lon,&place.wlon);
100 28994509 2004-04-21 devnull Xspalbers(&place,&x1,&y1);
101 28994509 2004-04-21 devnull scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));
105 28994509 2004-04-21 devnull invalb(double x, double y, double *lat, double *lon)
108 28994509 2004-04-21 devnull double sinb_den, sinp;
109 28994509 2004-04-21 devnull x *= scale;
110 28994509 2004-04-21 devnull y *= scale;
111 28994509 2004-04-21 devnull *lon = atan2(-x,fabs(y))/(RAD*n) + twist;
112 28994509 2004-04-21 devnull sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));
113 28994509 2004-04-21 devnull sinp = sinb_den;
114 28994509 2004-04-21 devnull for(i=0; i<5; i++)
115 28994509 2004-04-21 devnull sinp = sinb_den/num(sinp);
116 28994509 2004-04-21 devnull *lat = asin(sinp)/RAD;
