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