diff options
Diffstat (limited to 'src/cmd/map/libmap/homing.c')
-rw-r--r-- | src/cmd/map/libmap/homing.c | 121 |
1 files changed, 121 insertions, 0 deletions
diff --git a/src/cmd/map/libmap/homing.c b/src/cmd/map/libmap/homing.c new file mode 100644 index 00000000..366f69fe --- /dev/null +++ b/src/cmd/map/libmap/homing.c @@ -0,0 +1,121 @@ +#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; +} |