| #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; |
| } |