Return to loxodromic.c CVS log

Up to [Research Unix] / researchv10no / cmd / map / lox

Annotation of researchv10no/cmd/map/lox/loxodromic.c, revision 1.1

1.1     ! root        1: #include "map.h"
        !             2: #include <stdio.h>
        !             3: 
        !             4: /* development of a loxodrome from a unit sphere onto a plane
        !             5:    the loxodrome passes through 0,0 on the sphere,
        !             6:    inclined a (alpha) from the equator
        !             7: 
        !             8:    latitude p (phi) positive north
        !             9:    longitude t (theta) positive east - opposite to "map"
        !            10:    all angles in radians
        !            11: 
        !            12:    x,y are coords of developed curve
        !            13:    c (chi) is angle of inclination of developed curve from x axis
        !            14: 
        !            15:    dx/dp = cos(c)/sin(a)
        !            16:    dy/dp = sin(c)/sin(a)
        !            17:    dc/dp = tan(p)/tan(a)
        !            18:    dt/dp = 1/(cos(p)*tan(a))
        !            19: 
        !            20:    cc lox.c -lport -lF77 -lI77
        !            21:    a.out alpha         here alpha is in degrees
        !            22: */
        !            23: extern char *realloc();
        !            24: static double renorm();
        !            25: static bsearch(), getxy(), calcf(), solve();
        !            26: #define DEG (180/PI)
        !            27: #define MAX 85
        !            28: 
        !            29: /* curlicue is an ordered sequence of points
        !            30: along the midline of the spiral, in both geographic
        !            31: and developed coordinates */
        !            32: struct curl {
        !            33:        double phi,x,y,theta,chi;
        !            34: } *curlicue;   /* the midline */
        !            35: int nc;                /* number of points on the midline */
        !            36: 
        !            37: struct coord a;        /* alpha, angle from equator */
        !            38: static signa = 1;
        !            39: static double z[4];            /* x==z[0],y==z[1],t=z[2],c==z[3] */
        !            40: 
        !            41: /* inputs to dodes, port library integration routine */
        !            42: 
        !            43: static int n = 4;
        !            44: static double z0[4];           /* initial values = 0,0,0,a */
        !            45: static double p0;              /* initial value of p = 0 */
        !            46: static double p1 = 86*PI/180;  /* final value, 89.99+ degrees */
        !            47: static double dt;              /* current time step */
        !            48: static double deg = PI/180;    /* 1 degree */
        !            49: static float errpar[2] = {1.e-7,1.e-7};
        !            50: 
        !            51: static
        !            52: f(p,z,n,r)
        !            53: double *p, *z, *r;
        !            54: int *n;
        !            55: {
        !            56:        r[0] = cos(z[3])/a.s;
        !            57:        r[1] = sin(z[3])/a.s;
        !            58:        r[2] = a.c/(cos(*p)*a.s);
        !            59:        r[3] = sin(*p)*r[2];
        !            60: }
        !            61: static
        !            62: handle(t0,z0,t1,z1,n,dt,ts,e)
        !            63: double *t0,*z0,*t1,*z1,*dt,*ts;
        !            64: float *e;
        !            65: int *n;
        !            66: {
        !            67:        if(*dt > deg) *dt = deg;
        !            68:        if(fabs(z1[3]-z0[3])>.1 && *dt >.00001)
        !            69:                *dt *= .1/fabs(z1[3]-z0[3]);
        !            70:        if(*t1 + *dt >= PI/2)
        !            71:                *dt = (PI/2 - *t1)/2;
        !            72:        if(*t1 + *dt > *ts)
        !            73:                *dt = *ts - *t1 + .00001;
        !            74:        curlicue = (struct curl*)realloc((char*)curlicue,
        !            75:                        (nc+1)*sizeof(*curlicue));
        !            76:        curlicue[nc].phi = *t1;
        !            77:        curlicue[nc].x = z1[0];
        !            78:        curlicue[nc].y = z1[1];
        !            79:        curlicue[nc].theta = z1[2];
        !            80:        curlicue[nc].chi = z1[3];
        !            81:        nc++;
        !            82:        /*printf("%f %f %f %f %f\n",*t1*DEG,
        !            83:                z1[0]*DEG,z1[1]*DEG,z1[2]*DEG,z1[3]*DEG);
        !            84:        fflush(stdout);*/
        !            85:        if(*t0==*t1) abort();
        !            86: }
        !            87: static
        !            88: struct place negplace(place)
        !            89: struct place *place;
        !            90: {
        !            91:        struct place iplace;
        !            92:        iplace = *place;
        !            93:        iplace.nlat.l = -place->nlat.l;
        !            94:        iplace.nlat.s = -place->nlat.s;
        !            95:        iplace.wlon.l = -place->wlon.l;
        !            96:        iplace.wlon.s = -place->wlon.s;
        !            97:        return iplace;
        !            98: }
        !            99: static
        !           100: Xloxodromic(place, x, y)
        !           101: struct place *place;
        !           102: float *x, *y;
        !           103: {
        !           104:        double u1 = place->nlat.l;
        !           105:        double v1 = -place->wlon.l;
        !           106:        double d,p,t;
        !           107:        if(fabs(u1) > MAX*RAD) return 0;
        !           108:        if(u1 < 0) {
        !           109:                struct place iplace;
        !           110:                int ret;
        !           111:                iplace = negplace(place);
        !           112:                ret = Xloxodromic(&iplace,x,y);
        !           113:                *x = -*x;
        !           114:                *y = -*y;
        !           115:                return ret;
        !           116:        }
        !           117:        v1 = renorm(u1,bsearch(u1,0),v1);
        !           118:        solve(u1,v1,&d,&p,&t);
        !           119:        /*printf("%f %f %f %f %f %f\n",
        !           120:                u,v,v1*DEG,d*DEG,p*DEG,t*DEG);
        !           121:        fflush(stdout);*/
        !           122:        getxy(d,p,x,y);
        !           123:        /*printf("%f %f\n",x*DEG,y*DEG);
        !           124:        fflush(stdout);*/
        !           125:        return 1;
        !           126: }
        !           127: int (*loxodromic(par))()
        !           128: float par;
        !           129: {
        !           130:        deg2rad(par, &a);
        !           131:        if(par < 0)
        !           132:                signa = -1;
        !           133:        z0[3] = a.l;
        !           134:        dt = deg;
        !           135:        handle(&deg,z0,&p0,z0,&n,&deg,&p1,errpar);
        !           136:        dodes_(f,z0,&n,&p0,&p1,&dt,errpar,handle);
        !           137:        return Xloxodromic;
        !           138: }
        !           139: 
        !           140: /* find the index for latitude (i=0) or longitude (i=1)
        !           141: in array curlicue.  In case of longitude, u must have
        !           142: been normalized by renorm() */
        !           143: static
        !           144: bsearch(u,i)
        !           145: double u;
        !           146: {
        !           147:        int top = nc;
        !           148:        int bot = 0;
        !           149:        int mid;
        !           150:        while(bot<top-1) {
        !           151:                mid = (bot+top)/2;
        !           152:                if(i==0&&curlicue[mid].phi<=u||
        !           153:                   i==1&&(signa>0&&curlicue[mid].theta<=u||
        !           154:                          signa<0&&curlicue[mid].theta>=u))
        !           155:                        bot = mid;
        !           156:                else
        !           157:                        top = mid;
        !           158:        }
        !           159:        return bot;
        !           160: }
        !           161: /* wind v thru multiples of 2*PI until it lies within PI
        !           162: of the (linearly interpolated) midline */
        !           163: static double
        !           164: renorm(u,n,v)
        !           165: double u,v;
        !           166: {
        !           167:        double t,fract;
        !           168:        fract = (u-curlicue[n].phi)/
        !           169:                (curlicue[n+1].phi-curlicue[n].phi);
        !           170:        t = curlicue[n].theta*(1-fract) +
        !           171:                curlicue[n+1].theta*fract;
        !           172:        if(signa > 0) while(v<t-PI)
        !           173:                v += 2*PI;
        !           174:        else while(v>t+PI)
        !           175:                v -= 2*PI;
        !           176:        return v;
        !           177: }
        !           178: /* given the height (d) and foot (p,t) of the perpendicular projection
        !           179: of a point onto the midline, return the x-y coordinates of the point*/
        !           180: static
        !           181: getxy(d,p,xp,yp)
        !           182: float d,p,*xp,*yp;
        !           183: {
        !           184:        int n, sign = 1;
        !           185:        double fract, x0, y0, c0;
        !           186:        if(p < 0) {
        !           187:                p = -p;
        !           188:                sign = -1;
        !           189:        }
        !           190:        n = bsearch(p,0);
        !           191:        fract = (p-curlicue[n].phi)/
        !           192:                (curlicue[n+1].phi-curlicue[n].phi);
        !           193:        x0 = sign*(curlicue[n].x*(1-fract) +
        !           194:                curlicue[n+1].x*fract);
        !           195:        y0 = sign*(curlicue[n].y*(1-fract) +
        !           196:                curlicue[n+1].y*fract);
        !           197:        c0 = curlicue[n].chi*(1-fract) +
        !           198:                curlicue[n+1].chi*fract;
        !           199:        *xp = x0 - d*sin(c0);
        !           200:        *yp = y0 + d*cos(c0);
        !           201: }      
        !           202: 
        !           203: /*
        !           204:        solve simultaneously for d,p,t (dist to midline,
        !           205:                                        phi lat on midline,
        !           206:                                        theta lon on midline)
        !           207:        given a, u, v (rhumb angle from parallels, lat, east lon)
        !           208: 
        !           209:        sin d sin a = cos u sin(t-v)            law of sines
        !           210:        sin u = cos d sin p + sin d cos p cos a law of cosines
        !           211:        t tan a = log tan (p/2 + PI/4)          eqn of rhumb line
        !           212: 
        !           213:        initial guess: d=u-p, t=v, p from table lookup in curlicue
        !           214: */
        !           215: static double su,cu,v0;
        !           216: static int n3 = 3, jmax = 50;
        !           217: static double acc = 1e-6;
        !           218: static
        !           219: calcf(n,x,nf,g)
        !           220: int *n, *nf;
        !           221: double *x,*g;
        !           222: {
        !           223:        double f[3], sd;
        !           224:        if(x[1]>=PI/2) {
        !           225:                *nf = 0;
        !           226:                return;
        !           227:        }
        !           228:        sd = sin(x[0]);
        !           229:        f[0] = sd*a.s-cu*sin(x[2]-v0);
        !           230:        f[1] = su - cos(x[0])*sin(x[1]) - sd*cos(x[1])*a.c;
        !           231:        f[2] = x[2]*a.s/a.c - log(tan(x[1]/2+PI/4));
        !           232:        /*printf("calcf: d p t %g %f %f",
        !           233:                x[0]*DEG,x[1]*DEG,x[2]*DEG);
        !           234:        printf(" resid %f %f %f\n",f[0],f[1],f[2]);
        !           235:        fflush(stdout);*/
        !           236:        *g = f[0]*f[0]+f[1]*f[1]+f[2]*f[2];
        !           237: }
        !           238: /* given renormed lat-lon (u,v) find the height (*dp) and foot
        !           239: (*pp,*tp) of a perpendicular to the midline. precondition: u>=0 */
        !           240: static
        !           241: solve(u,v,dp,pp,tp)
        !           242: double u,v,*dp,*pp,*tp;
        !           243: {
        !           244:        double z[3];    /* d, p, t */
        !           245:        int n;
        !           246:        if(u < 0) {
        !           247:                abort();
        !           248:                solve(-u,-v,dp,pp,tp);
        !           249:                *dp = -*dp;
        !           250:                *pp = -*pp;
        !           251:                *tp = -*tp;
        !           252:                return;
        !           253:        }
        !           254:        n = bsearch(signa*fabs(v),1);
        !           255:        su = sin(u);
        !           256:        cu = cos(u);
        !           257:        z[1] = curlicue[n].phi; /* initial guess */
        !           258:        if(signa*v < 0) z[1] = -z[1];
        !           259:        z[0] = u - z[1];
        !           260:        z[2] = v0 = v;
        !           261:        /*printf("inputs (deg) %f %f\n",u*DEG,v*DEG);
        !           262:        printf("interval %d guess (deg) %f %f %f\n",n,
        !           263:                z[0]*DEG,z[1]*DEG,z[2]*DEG);
        !           264:        fflush(stdout);*/
        !           265:        dsmnf_(&n3,z,calcf,&jmax,&acc);
        !           266:        *dp = z[0];
        !           267:        *pp = z[1];
        !           268:        *tp = z[2];
        !           269: }
        !           270: loxcut(g,og,cutlon)
        !           271: struct place *g, *og;
        !           272: 
        !           273: {
        !           274:        double glon, oglon;
        !           275:        double sign = g->nlat.l>=0? 1: -1;
        !           276:        double nlat = sign*g->nlat.l;
        !           277:        if(nlat > MAX*RAD)
        !           278:                return 0;
        !           279:        glon = sign*renorm(nlat, bsearch(nlat,0),
        !           280:                                -sign*g->wlon.l);
        !           281:        sign = og->nlat.l>=0? 1: -1;
        !           282:        nlat = sign*og->nlat.l;
        !           283:        if(nlat > MAX*RAD)
        !           284:                return 0;
        !           285:        oglon = sign*renorm(nlat, bsearch(nlat,0),
        !           286:                                -sign*og->wlon.l);
        !           287:        return fabs(glon-oglon) < PI;
        !           288: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.