1 #include "map.h" 2 3 static struct coord p0; /* standard parallel */ 4 5 int first; 6 7 static double 8 trigclamp(double x) 9 { 10 return x>1? 1: x<-1? -1: x; 11 } 12 13 static struct coord az; /* azimuth of p0 seen from place */ 14 static struct coord rad; /* angular dist from place to p0 */ 15 16 static int 17 azimuth(struct place *place) 18 { 19 if(place->nlat.c < FUZZ) { 20 az.l = PI/2 + place->nlat.l - place->wlon.l; 21 sincos(&az); 22 rad.l = fabs(place->nlat.l - p0.l); 23 if(rad.l > PI) 24 rad.l = 2*PI - rad.l; 25 sincos(&rad); 26 return 1; 27 } 28 rad.c = trigclamp(p0.s*place->nlat.s + /* law of cosines */ 29 p0.c*place->nlat.c*place->wlon.c); 30 rad.s = sqrt(1 - rad.c*rad.c); 31 if(fabs(rad.s) < .001) { 32 az.s = 0; 33 az.c = 1; 34 } else { 35 az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */ 36 az.c = trigclamp((p0.s - rad.c*place->nlat.s) 37 /(rad.s*place->nlat.c)); 38 } 39 rad.l = atan2(rad.s, rad.c); 40 return 1; 41 } 42 43 static int 44 Xmecca(struct place *place, double *x, double *y) 45 { 46 if(!azimuth(place)) 47 return 0; 48 *x = -place->wlon.l; 49 *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s; 50 return fabs(*y)>2? -1: 51 rad.c<0? 0: 52 1; 53 } 54 55 proj 56 mecca(double par) 57 { 58 first = 1; 59 if(fabs(par)>80.) 60 return(0); 61 deg2rad(par,&p0); 62 return(Xmecca); 63 } 64 65 static int 66 Xhoming(struct place *place, double *x, double *y) 67 { 68 if(!azimuth(place)) 69 return 0; 70 *x = -rad.l*az.s; 71 *y = -rad.l*az.c; 72 return place->wlon.c<0? 0: 1; 73 } 74 75 proj 76 homing(double par) 77 { 78 first = 1; 79 if(fabs(par)>80.) 80 return(0); 81 deg2rad(par,&p0); 82 return(Xhoming); 83 } 84 85 int 86 hlimb(double *lat, double *lon, double res) 87 { 88 if(first) { 89 *lon = -90; 90 *lat = -90; 91 first = 0; 92 return 0; 93 } 94 *lat += res; 95 if(*lat <= 90) 96 return 1; 97 if(*lon == 90) 98 return -1; 99 *lon = 90; 100 *lat = -90; 101 return 0; 102 } 103 104 int 105 mlimb(double *lat, double *lon, double res) 106 { 107 int ret = !first; 108 if(fabs(p0.s) < .01) 109 return -1; 110 if(first) { 111 *lon = -180; 112 first = 0; 113 } else 114 *lon += res; 115 if(*lon > 180) 116 return -1; 117 *lat = atan(-cos(*lon*RAD)/p0.s*p0.c)/RAD; 118 return ret; 119 } 120