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