13e12c5d1SDavid du Colombier #include "map.h" 23e12c5d1SDavid du Colombier 33e12c5d1SDavid du Colombier static double 43e12c5d1SDavid du Colombier quadratic(double a, double b, double c) 53e12c5d1SDavid du Colombier { 63e12c5d1SDavid du Colombier double disc = b*b - 4*a*c; 73e12c5d1SDavid du Colombier return disc<0? 0: (-b-sqrt(disc))/(2*a); 83e12c5d1SDavid du Colombier } 93e12c5d1SDavid du Colombier 103e12c5d1SDavid du Colombier /* for projections with meridians being circles centered 113e12c5d1SDavid du Colombier on the x axis and parallels being circles centered on the y 123e12c5d1SDavid du Colombier axis. Find the intersection of the meridian thru (m,0), (90,0), 133e12c5d1SDavid du Colombier with the parallel thru (0,p), (p1,p2) */ 143e12c5d1SDavid du Colombier 153e12c5d1SDavid du Colombier static int 16*219b2ee8SDavid du Colombier twocircles(double m, double p, double p1, double p2, double *x, double *y) 173e12c5d1SDavid du Colombier { 18*219b2ee8SDavid du Colombier double a; /* center of meridian circle, a>0 */ 19*219b2ee8SDavid du Colombier double b; /* center of parallel circle, b>0 */ 20*219b2ee8SDavid du Colombier double t,bb; 213e12c5d1SDavid du Colombier if(m > 0) { 223e12c5d1SDavid du Colombier twocircles(-m,p,p1,p2,x,y); 233e12c5d1SDavid du Colombier *x = -*x; 243e12c5d1SDavid du Colombier } else if(p < 0) { 253e12c5d1SDavid du Colombier twocircles(m,-p,p1,-p2,x,y); 263e12c5d1SDavid du Colombier *y = -*y; 273e12c5d1SDavid du Colombier } else if(p < .01) { 283e12c5d1SDavid du Colombier *x = m; 293e12c5d1SDavid du Colombier t = m/p1; 303e12c5d1SDavid du Colombier *y = p + (p2-p)*t*t; 313e12c5d1SDavid du Colombier } else if(m > -.01) { 323e12c5d1SDavid du Colombier *y = p; 333e12c5d1SDavid du Colombier *x = m - m*p*p; 343e12c5d1SDavid du Colombier } else { 353e12c5d1SDavid du Colombier b = p>=1? 1: p>.99? 0.5*(p+1 + p1*p1/(1-p)): 363e12c5d1SDavid du Colombier 0.5*(p*p-p1*p1-p2*p2)/(p-p2); 373e12c5d1SDavid du Colombier a = .5*(m - 1/m); 383e12c5d1SDavid du Colombier t = m*m-p*p+2*(b*p-a*m); 393e12c5d1SDavid du Colombier bb = b*b; 403e12c5d1SDavid du Colombier *x = quadratic(1+a*a/bb, -2*a + a*t/bb, 413e12c5d1SDavid du Colombier t*t/(4*bb) - m*m + 2*a*m); 423e12c5d1SDavid du Colombier *y = (*x*a+t/2)/b; 433e12c5d1SDavid du Colombier } 443e12c5d1SDavid du Colombier return 1; 453e12c5d1SDavid du Colombier } 463e12c5d1SDavid du Colombier 473e12c5d1SDavid du Colombier static int 48*219b2ee8SDavid du Colombier Xglobular(struct place *place, double *x, double *y) 493e12c5d1SDavid du Colombier { 503e12c5d1SDavid du Colombier twocircles(-2*place->wlon.l/PI, 513e12c5d1SDavid du Colombier 2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y); 523e12c5d1SDavid du Colombier return 1; 533e12c5d1SDavid du Colombier } 543e12c5d1SDavid du Colombier 553e12c5d1SDavid du Colombier proj 563e12c5d1SDavid du Colombier globular(void) 573e12c5d1SDavid du Colombier { 583e12c5d1SDavid du Colombier return Xglobular; 593e12c5d1SDavid du Colombier } 603e12c5d1SDavid du Colombier 613e12c5d1SDavid du Colombier static int 62*219b2ee8SDavid du Colombier Xvandergrinten(struct place *place, double *x, double *y) 633e12c5d1SDavid du Colombier { 643e12c5d1SDavid du Colombier double t = 2*place->nlat.l/PI; 653e12c5d1SDavid du Colombier double abst = fabs(t); 663e12c5d1SDavid du Colombier double pval = abst>=1? 1: abst/(1+sqrt(1-t*t)); 673e12c5d1SDavid du Colombier double p2 = 2*pval/(1+pval); 683e12c5d1SDavid du Colombier twocircles(-place->wlon.l/PI, pval, sqrt(1-p2*p2), p2, x, y); 693e12c5d1SDavid du Colombier if(t < 0) 703e12c5d1SDavid du Colombier *y = -*y; 713e12c5d1SDavid du Colombier return 1; 723e12c5d1SDavid du Colombier } 733e12c5d1SDavid du Colombier 743e12c5d1SDavid du Colombier proj 753e12c5d1SDavid du Colombier vandergrinten(void) 763e12c5d1SDavid du Colombier { 773e12c5d1SDavid du Colombier return Xvandergrinten; 783e12c5d1SDavid du Colombier } 79