diff options
Diffstat (limited to 'src/cmd/map/libmap/gilbert.c')
-rw-r--r-- | src/cmd/map/libmap/gilbert.c | 51 |
1 files changed, 51 insertions, 0 deletions
diff --git a/src/cmd/map/libmap/gilbert.c b/src/cmd/map/libmap/gilbert.c new file mode 100644 index 00000000..173ffcd3 --- /dev/null +++ b/src/cmd/map/libmap/gilbert.c @@ -0,0 +1,51 @@ +#include <u.h> +#include <libc.h> +#include "map.h" + +int +Xgilbert(struct place *p, double *x, double *y) +{ +/* the interesting part - map the sphere onto a hemisphere */ + struct place q; + q.nlat.s = tan(0.5*(p->nlat.l)); + if(q.nlat.s > 1) q.nlat.s = 1; + if(q.nlat.s < -1) q.nlat.s = -1; + q.nlat.c = sqrt(1 - q.nlat.s*q.nlat.s); + q.wlon.l = p->wlon.l/2; + sincos(&q.wlon); +/* the dull part: present the hemisphere orthogrpahically */ + *y = q.nlat.s; + *x = -q.wlon.s*q.nlat.c; + return(1); +} + +proj +gilbert(void) +{ + return(Xgilbert); +} + +/* derivation of the interesting part: + map the sphere onto the plane by stereographic projection; + map the plane onto a half plane by sqrt; + map the half plane back to the sphere by stereographic + projection + + n,w are original lat and lon + r is stereographic radius + primes are transformed versions + + r = cos(n)/(1+sin(n)) + r' = sqrt(r) = cos(n')/(1+sin(n')) + + r'^2 = (1-sin(n')^2)/(1+sin(n')^2) = cos(n)/(1+sin(n)) + + this is a linear equation for sin n', with solution + + sin n' = (1+sin(n)-cos(n))/(1+sin(n)+cos(n)) + + use standard formula: tan x/2 = (1-cos x)/sin x = sin x/(1+cos x) + to show that the right side of the last equation is tan(n/2) +*/ + + |