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"
 4 28994509 2004-04-21 devnull 
 5 28994509 2004-04-21 devnull int
 6 28994509 2004-04-21 devnull Xgilbert(struct place *p, double *x, double *y)
 7 28994509 2004-04-21 devnull {
 8 28994509 2004-04-21 devnull /* the interesting part - map the sphere onto a hemisphere */
 9 28994509 2004-04-21 devnull 	struct place q;
10 28994509 2004-04-21 devnull 	q.nlat.s = tan(0.5*(p->nlat.l));
11 28994509 2004-04-21 devnull 	if(q.nlat.s > 1) q.nlat.s = 1;
12 28994509 2004-04-21 devnull 	if(q.nlat.s < -1) q.nlat.s = -1;
13 28994509 2004-04-21 devnull 	q.nlat.c = sqrt(1 - q.nlat.s*q.nlat.s);
14 28994509 2004-04-21 devnull 	q.wlon.l = p->wlon.l/2;
15 28994509 2004-04-21 devnull 	sincos(&q.wlon);
16 28994509 2004-04-21 devnull /* the dull part: present the hemisphere orthogrpahically */
17 28994509 2004-04-21 devnull 	*y = q.nlat.s;
18 28994509 2004-04-21 devnull 	*x = -q.wlon.s*q.nlat.c;
19 28994509 2004-04-21 devnull 	return(1);
20 28994509 2004-04-21 devnull }
21 28994509 2004-04-21 devnull 
22 28994509 2004-04-21 devnull proj
23 28994509 2004-04-21 devnull gilbert(void)
24 28994509 2004-04-21 devnull {
25 28994509 2004-04-21 devnull 	return(Xgilbert);
26 28994509 2004-04-21 devnull }
27 28994509 2004-04-21 devnull 
28 28994509 2004-04-21 devnull /* derivation of the interesting part:
29 28994509 2004-04-21 devnull    map the sphere onto the plane by stereographic projection;
30 28994509 2004-04-21 devnull    map the plane onto a half plane by sqrt;
31 28994509 2004-04-21 devnull    map the half plane back to the sphere by stereographic
32 28994509 2004-04-21 devnull    projection
33 28994509 2004-04-21 devnull 
34 28994509 2004-04-21 devnull    n,w are original lat and lon
35 28994509 2004-04-21 devnull    r is stereographic radius
36 28994509 2004-04-21 devnull    primes are transformed versions
37 28994509 2004-04-21 devnull 
38 28994509 2004-04-21 devnull    r = cos(n)/(1+sin(n))
39 28994509 2004-04-21 devnull    r' = sqrt(r) = cos(n')/(1+sin(n'))
40 28994509 2004-04-21 devnull 
41 28994509 2004-04-21 devnull    r'^2 = (1-sin(n')^2)/(1+sin(n')^2) = cos(n)/(1+sin(n))
42 28994509 2004-04-21 devnull 
43 28994509 2004-04-21 devnull    this is a linear equation for sin n', with solution
44 28994509 2004-04-21 devnull 
45 28994509 2004-04-21 devnull    sin n' = (1+sin(n)-cos(n))/(1+sin(n)+cos(n))
46 28994509 2004-04-21 devnull 
47 28994509 2004-04-21 devnull    use standard formula: tan x/2 = (1-cos x)/sin x = sin x/(1+cos x)
48 28994509 2004-04-21 devnull    to show that the right side of the last equation is tan(n/2)
49 28994509 2004-04-21 devnull */