twocirc.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 }