Return to lox.c CVS log

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

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

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

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