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