blob: 366f69fe440bbf95918ce28734151675c986599b [file] [log] [blame]
#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;
}