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