1*59cc4ca5SDavid du Colombier #include <u.h> 2*59cc4ca5SDavid du Colombier #include <libc.h> 3219b2ee8SDavid du Colombier #include "map.h" 4219b2ee8SDavid du Colombier 5219b2ee8SDavid du Colombier int Xgilbert(struct place * p,double * x,double * y)6219b2ee8SDavid du ColombierXgilbert(struct place *p, double *x, double *y) 7219b2ee8SDavid du Colombier { 8219b2ee8SDavid du Colombier /* the interesting part - map the sphere onto a hemisphere */ 9219b2ee8SDavid du Colombier struct place q; 10219b2ee8SDavid du Colombier q.nlat.s = tan(0.5*(p->nlat.l)); 11219b2ee8SDavid du Colombier if(q.nlat.s > 1) q.nlat.s = 1; 12219b2ee8SDavid du Colombier if(q.nlat.s < -1) q.nlat.s = -1; 13219b2ee8SDavid du Colombier q.nlat.c = sqrt(1 - q.nlat.s*q.nlat.s); 14219b2ee8SDavid du Colombier q.wlon.l = p->wlon.l/2; 15219b2ee8SDavid du Colombier sincos(&q.wlon); 16219b2ee8SDavid du Colombier /* the dull part: present the hemisphere orthogrpahically */ 17219b2ee8SDavid du Colombier *y = q.nlat.s; 18219b2ee8SDavid du Colombier *x = -q.wlon.s*q.nlat.c; 19219b2ee8SDavid du Colombier return(1); 20219b2ee8SDavid du Colombier } 21219b2ee8SDavid du Colombier 22219b2ee8SDavid du Colombier proj gilbert(void)23219b2ee8SDavid du Colombiergilbert(void) 24219b2ee8SDavid du Colombier { 25219b2ee8SDavid du Colombier return(Xgilbert); 26219b2ee8SDavid du Colombier } 27219b2ee8SDavid du Colombier 28219b2ee8SDavid du Colombier /* derivation of the interesting part: 29219b2ee8SDavid du Colombier map the sphere onto the plane by stereographic projection; 30219b2ee8SDavid du Colombier map the plane onto a half plane by sqrt; 31219b2ee8SDavid du Colombier map the half plane back to the sphere by stereographic 32219b2ee8SDavid du Colombier projection 33219b2ee8SDavid du Colombier 34219b2ee8SDavid du Colombier n,w are original lat and lon 35219b2ee8SDavid du Colombier r is stereographic radius 36219b2ee8SDavid du Colombier primes are transformed versions 37219b2ee8SDavid du Colombier 38219b2ee8SDavid du Colombier r = cos(n)/(1+sin(n)) 39219b2ee8SDavid du Colombier r' = sqrt(r) = cos(n')/(1+sin(n')) 40219b2ee8SDavid du Colombier 41219b2ee8SDavid du Colombier r'^2 = (1-sin(n')^2)/(1+sin(n')^2) = cos(n)/(1+sin(n)) 42219b2ee8SDavid du Colombier 43219b2ee8SDavid du Colombier this is a linear equation for sin n', with solution 44219b2ee8SDavid du Colombier 45219b2ee8SDavid du Colombier sin n' = (1+sin(n)-cos(n))/(1+sin(n)+cos(n)) 46219b2ee8SDavid du Colombier 47219b2ee8SDavid du Colombier use standard formula: tan x/2 = (1-cos x)/sin x = sin x/(1+cos x) 48219b2ee8SDavid du Colombier to show that the right side of the last equation is tan(n/2) 49219b2ee8SDavid du Colombier */ 50219b2ee8SDavid du Colombier 51219b2ee8SDavid du Colombier 52