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