Subversion Repositories planix.SVN

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