1 28994509 2004-04-21 devnull #include <u.h>
2 28994509 2004-04-21 devnull #include <libc.h>
3 28994509 2004-04-21 devnull #include "map.h"
5 28994509 2004-04-21 devnull static double
6 28994509 2004-04-21 devnull quadratic(double a, double b, double c)
8 28994509 2004-04-21 devnull double disc = b*b - 4*a*c;
9 28994509 2004-04-21 devnull return disc<0? 0: (-b-sqrt(disc))/(2*a);
12 28994509 2004-04-21 devnull /* for projections with meridians being circles centered
13 28994509 2004-04-21 devnull on the x axis and parallels being circles centered on the y
14 28994509 2004-04-21 devnull axis. Find the intersection of the meridian thru (m,0), (90,0),
15 28994509 2004-04-21 devnull with the parallel thru (0,p), (p1,p2) */
17 28994509 2004-04-21 devnull static int
18 28994509 2004-04-21 devnull twocircles(double m, double p, double p1, double p2, double *x, double *y)
20 28994509 2004-04-21 devnull double a; /* center of meridian circle, a>0 */
21 28994509 2004-04-21 devnull double b; /* center of parallel circle, b>0 */
22 28994509 2004-04-21 devnull double t,bb;
23 28994509 2004-04-21 devnull if(m > 0) {
24 28994509 2004-04-21 devnull twocircles(-m,p,p1,p2,x,y);
25 28994509 2004-04-21 devnull *x = -*x;
26 28994509 2004-04-21 devnull } else if(p < 0) {
27 28994509 2004-04-21 devnull twocircles(m,-p,p1,-p2,x,y);
28 28994509 2004-04-21 devnull *y = -*y;
29 28994509 2004-04-21 devnull } else if(p < .01) {
31 28994509 2004-04-21 devnull t = m/p1;
32 28994509 2004-04-21 devnull *y = p + (p2-p)*t*t;
33 28994509 2004-04-21 devnull } else if(m > -.01) {
35 28994509 2004-04-21 devnull *x = m - m*p*p;
37 28994509 2004-04-21 devnull b = p>=1? 1: p>.99? 0.5*(p+1 + p1*p1/(1-p)):
38 28994509 2004-04-21 devnull 0.5*(p*p-p1*p1-p2*p2)/(p-p2);
39 28994509 2004-04-21 devnull a = .5*(m - 1/m);
40 28994509 2004-04-21 devnull t = m*m-p*p+2*(b*p-a*m);
41 28994509 2004-04-21 devnull bb = b*b;
42 28994509 2004-04-21 devnull *x = quadratic(1+a*a/bb, -2*a + a*t/bb,
43 28994509 2004-04-21 devnull t*t/(4*bb) - m*m + 2*a*m);
44 28994509 2004-04-21 devnull *y = (*x*a+t/2)/b;
46 28994509 2004-04-21 devnull return 1;
49 28994509 2004-04-21 devnull static int
50 28994509 2004-04-21 devnull Xglobular(struct place *place, double *x, double *y)
52 28994509 2004-04-21 devnull twocircles(-2*place->wlon.l/PI,
53 28994509 2004-04-21 devnull 2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y);
54 28994509 2004-04-21 devnull return 1;
58 28994509 2004-04-21 devnull globular(void)
60 28994509 2004-04-21 devnull return Xglobular;
63 28994509 2004-04-21 devnull static int
64 28994509 2004-04-21 devnull Xvandergrinten(struct place *place, double *x, double *y)
66 28994509 2004-04-21 devnull double t = 2*place->nlat.l/PI;
67 28994509 2004-04-21 devnull double abst = fabs(t);
68 28994509 2004-04-21 devnull double pval = abst>=1? 1: abst/(1+sqrt(1-t*t));
69 28994509 2004-04-21 devnull double p2 = 2*pval/(1+pval);
70 28994509 2004-04-21 devnull twocircles(-place->wlon.l/PI, pval, sqrt(1-p2*p2), p2, x, y);
71 28994509 2004-04-21 devnull if(t < 0)
72 28994509 2004-04-21 devnull *y = -*y;
73 28994509 2004-04-21 devnull return 1;
77 28994509 2004-04-21 devnull vandergrinten(void)
79 28994509 2004-04-21 devnull return Xvandergrinten;
