Return to tetra.c CVS log

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

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

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

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