Return to albers.c CVS log

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

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

1.1     ! root        1: #include "map.h"
        !             2: 
        !             3: /* For Albers formulas see Deetz and Adams "Elements of Map Projection", */
        !             4: /* USGS Special Publication No. 68, GPO 1921 */
        !             5: 
        !             6: static float r0sq, r1sq, d2, n, den, sinb1, sinb2;
        !             7: static struct place plat1, plat2;
        !             8: static southpole;
        !             9: 
        !            10: static float num(s)
        !            11: float s;
        !            12: {
        !            13:        if(d2==0)
        !            14:                return(1);
        !            15:        s = d2*s*s;
        !            16:        return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));
        !            17: }
        !            18: 
        !            19: /* Albers projection for a spheroid, good only when N pole is fixed */
        !            20: 
        !            21: static int
        !            22: Xspalbers(place,x,y)
        !            23: struct place *place;
        !            24: float *x, *y;
        !            25: {
        !            26:        float r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);
        !            27:        float t = n*place->wlon.l;
        !            28:        *y = r*cos(t);
        !            29:        if(!southpole)
        !            30:                *y = -*y;
        !            31:        *x = -r*sin(t);
        !            32:        return(1);
        !            33: }
        !            34: 
        !            35: /* lat1, lat2: std parallels; e2: squared eccentricity */
        !            36: 
        !            37: int (*albinit(lat1,lat2,e2))()
        !            38: float lat1,lat2,e2;
        !            39: {
        !            40:        float r1,r2;
        !            41:        extern (*azequalarea())();
        !            42:        extern (*cylequalarea())();
        !            43:        float t;
        !            44:        for(;;) {
        !            45:                if(lat1 < -90)
        !            46:                        lat1 = -180 - lat1;
        !            47:                if(lat2 > 90)
        !            48:                        lat2 = 180 - lat2;
        !            49:                if(lat1 <= lat2)
        !            50:                        break;
        !            51:                t = lat1; lat1 = lat2; lat2 = t;
        !            52:        }
        !            53:        if(lat2-lat1 < 1) {
        !            54:                if(lat1 > 89)
        !            55:                        return(azequalarea());
        !            56:                return(0);
        !            57:        }
        !            58:        if(fabs(lat2+lat1) < 1)
        !            59:                return(cylequalarea(lat1));
        !            60:        d2 = e2;
        !            61:        den = num(1.);
        !            62:        deg2rad(lat1,&plat1);
        !            63:        deg2rad(lat2,&plat2);
        !            64:        sinb1 = plat1.nlat.s*num(plat1.nlat.s)/den;
        !            65:        sinb2 = plat2.nlat.s*num(plat2.nlat.s)/den;
        !            66:        n = (plat1.nlat.c*plat1.nlat.c/(1-e2*plat1.nlat.s*plat1.nlat.s) -
        !            67:            plat2.nlat.c*plat2.nlat.c/(1-e2*plat2.nlat.s*plat2.nlat.s)) /
        !            68:            (2*(1-e2)*den*(sinb2-sinb1));
        !            69:        r1 = plat1.nlat.c/(n*sqrt(1-e2*plat1.nlat.s*plat1.nlat.s));
        !            70:        r2 = plat2.nlat.c/(n*sqrt(2-e2*plat2.nlat.s*plat2.nlat.s));
        !            71:        r1sq = r1*r1;
        !            72:        r0sq = r1sq + 2*(1-e2)*den*sinb1/n;
        !            73:        southpole = lat1<0 && plat2.nlat.c>plat1.nlat.c;
        !            74:        return(Xspalbers);
        !            75: }
        !            76: 
        !            77: int (*sp_albers(lat1,lat2))()
        !            78: float lat1, lat2;
        !            79: {
        !            80:        return(albinit(lat1,lat2,EC2));
        !            81: }
        !            82: 
        !            83: int (*albers(lat1,lat2))()
        !            84: float lat1,lat2;
        !            85: {
        !            86:        return(albinit(lat1,lat2,0.));
        !            87: }
        !            88: 
        !            89: static float scale = 1;
        !            90: static float twist = 0;
        !            91: 
        !            92: albscale(x,y,lat,lon)
        !            93: float x,y,lat,lon;
        !            94: {
        !            95:        struct place place;
        !            96:        float alat, alon, x1,y1;
        !            97:        scale = 1;
        !            98:        twist = 0;
        !            99:        invalb(x,y,&alat,&alon);
        !           100:        twist = lon - alon;
        !           101:        deg2rad(lat,&place.nlat);
        !           102:        deg2rad(lon,&place.wlon);
        !           103:        Xspalbers(&place,&x1,&y1);
        !           104:        scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));
        !           105: }
        !           106: 
        !           107: invalb(x,y,lat,lon)
        !           108: float x,y,*lat,*lon;
        !           109: {
        !           110:        int i;
        !           111:        float sinb_den, sinp;
        !           112:        x *= scale;
        !           113:        y *= scale;
        !           114:        *lon = atan2(-x,fabs(y))/(RAD*n) + twist;
        !           115:        sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));
        !           116:        sinp = sinb_den;
        !           117:        for(i=0; i<5; i++)
        !           118:                sinp = sinb_den/num(sinp);
        !           119:        *lat = asin(sinp)/RAD;
        !           120: }

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