blob: f7ee7b6f8497b04db92db4b15d6812c43d451b29 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
|
#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)
*/
|