6 quadratic(double a, double b, double c)
8 double disc = b*b - 4*a*c;
9 return disc<0? 0: (-b-sqrt(disc))/(2*a);
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) */
18 twocircles(double m, double p, double p1, double p2, double *x, double *y)
20 double a; /* center of meridian circle, a>0 */
21 double b; /* center of parallel circle, b>0 */
24 twocircles(-m,p,p1,p2,x,y);
27 twocircles(m,-p,p1,-p2,x,y);
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);
40 t = m*m-p*p+2*(b*p-a*m);
42 *x = quadratic(1+a*a/bb, -2*a + a*t/bb,
43 t*t/(4*bb) - m*m + 2*a*m);
50 Xglobular(struct place *place, double *x, double *y)
52 twocircles(-2*place->wlon.l/PI,
53 2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y);
64 Xvandergrinten(struct place *place, double *x, double *y)
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);
79 return Xvandergrinten;