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