gilbert.c (1208B)
1 #include <u.h> 2 #include <libc.h> 3 #include "map.h" 4 5 int 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 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 */