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

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.