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