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