Return to albers.c CVS log

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

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

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