albers.c (2401B)
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 int southpole; 11 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 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 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 75 sp_albers(double lat1, double lat2) 76 { 77 return(albinit(lat1,lat2,EC2)); 78 } 79 80 proj 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 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 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 }