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 }