Return to tetra.c CVS log

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

Annotation of researchv10no/cmd/map/libmap/tetra.c, revision 1.1

1.1     ! root        1: #include "map.h"
        !             2: 
        !             3: /*
        !             4:  *     conformal map of earth onto tetrahedron
        !             5:  *     the stages of mapping are
        !             6:  *     (a) stereo projection of tetrahedral face onto
        !             7:  *         isosceles curvilinear triangle with 3 120-degree
        !             8:  *         angles and one straight side
        !             9:  *     (b) map of this triangle onto half plane cut along
        !            10:  *         3 rays from the roots of unity to infinity
        !            11:  *             formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1)
        !            12:  *     (c) do 3 times for  each sector of plane:
        !            13:  *         map of |arg z|<=pi/6, cut along z>1 into
        !            14:  *         triangle |arg z|<=pi/6, Re z<=const,
        !            15:  *         with upper side of cut going into upper half of
        !            16:  *         of vertical side of triangle and lowere into lower
        !            17:  *             formula int from 0 to z dz/sqrt(1-z^3)
        !            18:  *
        !            19:  *     int from u to 1 3^.25*du/sqrt(1-u^3) =
        !            20:                F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4))
        !            21:  *     int from 1 to u 3^.25*du/sqrt(u^3-1) =
        !            22:  *             F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4))
        !            23:  *     this latter formula extends analytically down to
        !            24:  *     u=0 and is the basis of this routine, with the
        !            25:  *     argument of complex elliptic integral elco2
        !            26:  *     being tan(acos...)
        !            27:  *     the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is
        !            28:  *     used to cross over into the region where Re(acos...)>pi/2
        !            29:  *             f0 and fpi are suitably scaled complete integrals
        !            30: */
        !            31: 
        !            32: #define TFUZZ 0.00001
        !            33: 
        !            34: static struct place tpole[4];  /* point of tangency of tetrahedron face*/
        !            35: static double tpoleinit[4][2] = {
        !            36:        1.,     0.,
        !            37:        1.,     180.,
        !            38:        -1.,    90.,
        !            39:        -1.,    -90.
        !            40: };
        !            41: static struct tproj {
        !            42:        double tlat,tlon;       /* center of stereo projection*/
        !            43:        double ttwist;          /* rotatn before stereo*/
        !            44:        double trot;            /*rotate after projection*/
        !            45:        struct place projpl;    /*same as tlat,tlon*/
        !            46:        struct coord projtw;    /*same as ttwist*/
        !            47:        struct coord postrot;   /*same as trot*/
        !            48: } tproj[4][4] = {
        !            49: {/*00*/        {0.},
        !            50:  /*01*/        {90.,   0.,     90.,    -90.},
        !            51:  /*02*/        {0.,    45.,    -45.,   150.},
        !            52:  /*03*/        {0.,    -45.,   -135.,  30.}
        !            53: },
        !            54: {/*10*/        {90.,   0.,     -90.,   90.},
        !            55:  /*11*/ {0.},
        !            56:  /*12*/ {0.,   135.,   -135.,  -150.},
        !            57:  /*13*/        {0.,    -135.,  -45.,   -30.}
        !            58: },
        !            59: {/*20*/        {0.,    45.,    135.,   -30.},
        !            60:  /*21*/        {0.,    135.,   45.,    -150.},
        !            61:  /*22*/        {0.},
        !            62:  /*23*/        {-90.,  0.,     180.,   90.}
        !            63: },
        !            64: {/*30*/        {0.,    -45.,   45.,    -150.},
        !            65:  /*31*/ {0.,   -135.,  135.,   -30.},
        !            66:  /*32*/        {-90.,  0.,     0.,     90.},
        !            67:  /*33*/ {0.} 
        !            68: }};
        !            69: static double tx[4] = {        /*where to move facet after final rotation*/
        !            70:        0.,     0.,     -1.,    1.      /*-1,1 to be sqrt(3)*/
        !            71: };
        !            72: static double ty[4] = {
        !            73:        0.,     2.,     -1.,    -1.
        !            74: };
        !            75: static double root3;
        !            76: static double rt3inv;
        !            77: static double two_rt3;
        !            78: static double tkc,tk,tcon;
        !            79: static double f0r,f0i,fpir,fpii;
        !            80: 
        !            81: static void
        !            82: twhichp(struct place *g, int *p, int *q)
        !            83: {
        !            84:        int i,j,k;
        !            85:        double cosdist[4];
        !            86:        struct place *tp;
        !            87:        for(i=0;i<4;i++) {
        !            88:                tp = &tpole[i];
        !            89:                cosdist[i] = g->nlat.s*tp->nlat.s +
        !            90:                          g->nlat.c*tp->nlat.c*(
        !            91:                          g->wlon.s*tp->wlon.s +
        !            92:                          g->wlon.c*tp->wlon.c);
        !            93:        }
        !            94:        j = 0;
        !            95:        for(i=1;i<4;i++)
        !            96:                if(cosdist[i] > cosdist[j])
        !            97:                        j = i;
        !            98:        *p = j;
        !            99:        k = j==0?1:0;
        !           100:        for(i=0;i<4;i++)
        !           101:                if(i!=j&&cosdist[i]>cosdist[k])
        !           102:                        k = i;
        !           103:        *q = k;
        !           104: }
        !           105: 
        !           106: int
        !           107: Xtetra(struct place *place, double *x, double *y)
        !           108: {
        !           109:        int i,j;
        !           110:        struct place pl;
        !           111:        register struct tproj *tpp;
        !           112:        double vr, vi;
        !           113:        double br, bi;
        !           114:        double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti;
        !           115:        twhichp(place,&i,&j);
        !           116:        copyplace(place,&pl);
        !           117:        norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw);
        !           118:        Xstereographic(&pl,&vr,&vi);
        !           119:        zr = vr/2;
        !           120:        zi = vi/2;
        !           121:        if(zr<=TFUZZ)
        !           122:                zr = TFUZZ;
        !           123:        csq(zr,zi,&z2r,&z2i);
        !           124:        csq(z2r,z2i,&z4r,&z4i);
        !           125:        z2r *= two_rt3;
        !           126:        z2i *= two_rt3;
        !           127:        cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si);
        !           128:        csqrt(sr-1,si,&tr,&ti);
        !           129:        cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi);
        !           130:        if(br<0) {
        !           131:                br = -br;
        !           132:                bi = -bi;
        !           133:                if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
        !           134:                        return 0;
        !           135:                vr = fpir - vr;
        !           136:                vi = fpii - vi;
        !           137:        } else 
        !           138:                if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
        !           139:                        return 0;
        !           140:        if(si>=0) {
        !           141:                tr = f0r - vi;
        !           142:                ti = f0i + vr;
        !           143:        } else {
        !           144:                tr = f0r + vi;
        !           145:                ti = f0i - vr;
        !           146:        }
        !           147:        tpp = &tproj[i][j];
        !           148:        *x = tr*tpp->postrot.c +
        !           149:             ti*tpp->postrot.s + tx[i];
        !           150:        *y = ti*tpp->postrot.c -
        !           151:             tr*tpp->postrot.s + ty[i];
        !           152:        return(1);
        !           153: }
        !           154: 
        !           155: int
        !           156: tetracut(struct place *g, struct place *og, double *cutlon)
        !           157: {
        !           158:        int i,j,k;
        !           159:        if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) && 
        !           160:           (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
        !           161:                return(2);
        !           162:        twhichp(g,&i,&k);
        !           163:        twhichp(og,&j,&k);
        !           164:        if(i==j||i==0||j==0)
        !           165:                return(1);
        !           166:        return(0);
        !           167: }
        !           168: 
        !           169: proj
        !           170: tetra(void)
        !           171: {
        !           172:        register i;
        !           173:        int j;
        !           174:        register struct place *tp;
        !           175:        register struct tproj *tpp;
        !           176:        double t;
        !           177:        root3 = sqrt(3.);
        !           178:        rt3inv = 1/root3;
        !           179:        two_rt3 = 2*root3;
        !           180:        tkc = sqrt(.5-.25*root3);
        !           181:        tk = sqrt(.5+.25*root3);
        !           182:        tcon = 2*sqrt(root3);
        !           183:        elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i);
        !           184:        elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii);
        !           185:        fpir *= 2;
        !           186:        fpii *= 2;
        !           187:        for(i=0;i<4;i++) {
        !           188:                tx[i] *= f0r*root3;
        !           189:                ty[i] *= f0r;
        !           190:                tp = &tpole[i];
        !           191:                t = tp->nlat.s = tpoleinit[i][0]/root3;
        !           192:                tp->nlat.c = sqrt(1 - t*t);
        !           193:                tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c);
        !           194:                deg2rad(tpoleinit[i][1],&tp->wlon);
        !           195:                for(j=0;j<4;j++) {
        !           196:                        tpp = &tproj[i][j];
        !           197:                        latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
        !           198:                        deg2rad(tpp->ttwist,&tpp->projtw);
        !           199:                        deg2rad(tpp->trot,&tpp->postrot);
        !           200:                }
        !           201:        }
        !           202:        return(Xtetra);
        !           203: }
        !           204: 

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