tgilbert.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tgilbert.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 */