Return to twocirc.c CVS log

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

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

1.1     ! root        1: #include "map.h"
        !             2: 
        !             3: static double
        !             4: quadratic(double a, double b, double c)
        !             5: {
        !             6:        double disc = b*b - 4*a*c;
        !             7:        return disc<0? 0: (-b-sqrt(disc))/(2*a);
        !             8: }
        !             9: 
        !            10: /* for projections with meridians being circles centered
        !            11: on the x axis and parallels being circles centered on the y
        !            12: axis.  Find the intersection of the meridian thru (m,0), (90,0),
        !            13: with the parallel thru (0,p), (p1,p2) */
        !            14: 
        !            15: static int
        !            16: twocircles(double m, double p, double p1, double p2, double *x, double *y)
        !            17: {
        !            18:        double a;       /* center of meridian circle, a>0 */
        !            19:        double b;       /* center of parallel circle, b>0 */
        !            20:        double t,bb;
        !            21:        if(m > 0) {
        !            22:                twocircles(-m,p,p1,p2,x,y);
        !            23:                *x = -*x;
        !            24:        } else if(p < 0) {
        !            25:                twocircles(m,-p,p1,-p2,x,y);
        !            26:                *y = -*y;
        !            27:        } else if(p < .01) {
        !            28:                *x = m;
        !            29:                t = m/p1;
        !            30:                *y = p + (p2-p)*t*t;
        !            31:        } else if(m > -.01) {
        !            32:                *y = p;
        !            33:                *x = m - m*p*p;
        !            34:        } else {
        !            35:                b = p>=1? 1: p>.99? 0.5*(p+1 + p1*p1/(1-p)):
        !            36:                        0.5*(p*p-p1*p1-p2*p2)/(p-p2);
        !            37:                a = .5*(m - 1/m);
        !            38:                t = m*m-p*p+2*(b*p-a*m);
        !            39:                bb = b*b;
        !            40:                *x = quadratic(1+a*a/bb, -2*a + a*t/bb,
        !            41:                        t*t/(4*bb) - m*m + 2*a*m);
        !            42:                *y = (*x*a+t/2)/b;
        !            43:        }
        !            44:        return 1;
        !            45: }              
        !            46: 
        !            47: static int
        !            48: Xglobular(struct place *place, double *x, double *y)
        !            49: {
        !            50:        twocircles(-2*place->wlon.l/PI,
        !            51:                2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y);
        !            52:        return 1;
        !            53: }      
        !            54: 
        !            55: proj
        !            56: globular(void)
        !            57: {
        !            58:        return Xglobular;
        !            59: }
        !            60: 
        !            61: static int
        !            62: Xvandergrinten(struct place *place, double *x, double *y)
        !            63: {
        !            64:        double t = 2*place->nlat.l/PI;
        !            65:        double abst = fabs(t);
        !            66:        double pval = abst>=1? 1: abst/(1+sqrt(1-t*t));
        !            67:        double p2 = 2*pval/(1+pval);
        !            68:        twocircles(-place->wlon.l/PI, pval, sqrt(1-p2*p2), p2, x, y);
        !            69:        if(t < 0) 
        !            70:                *y = -*y;
        !            71:        return 1;
        !            72: }
        !            73: 
        !            74: proj
        !            75: vandergrinten(void)
        !            76: {
        !            77:        return Xvandergrinten;
        !            78: }

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