ref: dfb4b522830edab8f3856289d326d6cf5e930644
dir: /sys/src/cmd/map/libmap/homing.c/
#include <u.h> #include <libc.h> #include "map.h" static struct coord p0; /* standard parallel */ int first; static double trigclamp(double x) { return x>1? 1: x<-1? -1: x; } static struct coord az; /* azimuth of p0 seen from place */ static struct coord rad; /* angular dist from place to p0 */ static int azimuth(struct place *place) { if(place->nlat.c < FUZZ) { az.l = PI/2 + place->nlat.l - place->wlon.l; sincos(&az); rad.l = fabs(place->nlat.l - p0.l); if(rad.l > PI) rad.l = 2*PI - rad.l; sincos(&rad); return 1; } rad.c = trigclamp(p0.s*place->nlat.s + /* law of cosines */ p0.c*place->nlat.c*place->wlon.c); rad.s = sqrt(1 - rad.c*rad.c); if(fabs(rad.s) < .001) { az.s = 0; az.c = 1; } else { az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */ az.c = trigclamp((p0.s - rad.c*place->nlat.s) /(rad.s*place->nlat.c)); } rad.l = atan2(rad.s, rad.c); return 1; } static int Xmecca(struct place *place, double *x, double *y) { if(!azimuth(place)) return 0; *x = -place->wlon.l; *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s; return fabs(*y)>2? -1: rad.c<0? 0: 1; } proj mecca(double par) { first = 1; if(fabs(par)>80.) return(0); deg2rad(par,&p0); return(Xmecca); } static int Xhoming(struct place *place, double *x, double *y) { if(!azimuth(place)) return 0; *x = -rad.l*az.s; *y = -rad.l*az.c; return place->wlon.c<0? 0: 1; } proj homing(double par) { first = 1; if(fabs(par)>80.) return(0); deg2rad(par,&p0); return(Xhoming); } int hlimb(double *lat, double *lon, double res) { if(first) { *lon = -90; *lat = -90; first = 0; return 0; } *lat += res; if(*lat <= 90) return 1; if(*lon == 90) return -1; *lon = 90; *lat = -90; return 0; } int mlimb(double *lat, double *lon, double res) { int ret = !first; if(fabs(p0.s) < .01) return -1; if(first) { *lon = -180; first = 0; } else *lon += res; if(*lon > 180) return -1; *lat = atan(-cos(*lon*RAD)/p0.s*p0.c)/RAD; return ret; }