xref: /plan9/sys/src/cmd/map/libmap/homing.c (revision 59cc4ca53493a3c6d2349fe2b7f7c40f7dce7294)

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