6 * conformal map of earth onto tetrahedron
7 * the stages of mapping are
8 * (a) stereo projection of tetrahedral face onto
9 * isosceles curvilinear triangle with 3 120-degree
10 * angles and one straight side
11 * (b) map of this triangle onto half plane cut along
12 * 3 rays from the roots of unity to infinity
13 * formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1)
14 * (c) do 3 times for each sector of plane:
15 * map of |arg z|<=pi/6, cut along z>1 into
16 * triangle |arg z|<=pi/6, Re z<=const,
17 * with upper side of cut going into upper half of
18 * of vertical side of triangle and lowere into lower
19 * formula int from 0 to z dz/sqrt(1-z^3)
21 * int from u to 1 3^.25*du/sqrt(1-u^3) =
22 F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4))
23 * int from 1 to u 3^.25*du/sqrt(u^3-1) =
24 * F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4))
25 * this latter formula extends analytically down to
26 * u=0 and is the basis of this routine, with the
27 * argument of complex elliptic integral elco2
29 * the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is
30 * used to cross over into the region where Re(acos...)>pi/2
31 * f0 and fpi are suitably scaled complete integrals
36 static struct place tpole[4]; /* point of tangency of tetrahedron face*/
37 static double tpoleinit[4][2] = {
44 double tlat,tlon; /* center of stereo projection*/
45 double ttwist; /* rotatn before stereo*/
46 double trot; /*rotate after projection*/
47 struct place projpl; /*same as tlat,tlon*/
48 struct coord projtw; /*same as ttwist*/
49 struct coord postrot; /*same as trot*/
52 /*01*/ {90., 0., 90., -90.},
53 /*02*/ {0., 45., -45., 150.},
54 /*03*/ {0., -45., -135., 30.}
56 {/*10*/ {90., 0., -90., 90.},
58 /*12*/ {0., 135., -135., -150.},
59 /*13*/ {0., -135., -45., -30.}
61 {/*20*/ {0., 45., 135., -30.},
62 /*21*/ {0., 135., 45., -150.},
64 /*23*/ {-90., 0., 180., 90.}
66 {/*30*/ {0., -45., 45., -150.},
67 /*31*/ {0., -135., 135., -30.},
68 /*32*/ {-90., 0., 0., 90.},
71 static double tx[4] = { /*where to move facet after final rotation*/
72 0., 0., -1., 1. /*-1,1 to be sqrt(3)*/
74 static double ty[4] = {
79 static double two_rt3;
80 static double tkc,tk,tcon;
81 static double f0r,f0i,fpir,fpii;
84 twhichp(struct place *g, int *p, int *q)
91 cosdist[i] = g->nlat.s*tp->nlat.s +
92 g->nlat.c*tp->nlat.c*(
93 g->wlon.s*tp->wlon.s +
94 g->wlon.c*tp->wlon.c);
98 if(cosdist[i] > cosdist[j])
103 if(i!=j&&cosdist[i]>cosdist[k])
109 Xtetra(struct place *place, double *x, double *y)
113 register struct tproj *tpp;
116 double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti;
117 twhichp(place,&i,&j);
118 copyplace(place,&pl);
119 norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw);
120 Xstereographic(&pl,&vr,&vi);
125 csq(zr,zi,&z2r,&z2i);
126 csq(z2r,z2i,&z4r,&z4i);
129 cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si);
130 csqrt(sr-1,si,&tr,&ti);
131 cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi);
135 if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
140 if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
150 *x = tr*tpp->postrot.c +
151 ti*tpp->postrot.s + tx[i];
152 *y = ti*tpp->postrot.c -
153 tr*tpp->postrot.s + ty[i];
158 tetracut(struct place *g, struct place *og, double *cutlon)
161 if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) &&
162 (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
176 register struct place *tp;
177 register struct tproj *tpp;
182 tkc = sqrt(.5-.25*root3);
183 tk = sqrt(.5+.25*root3);
184 tcon = 2*sqrt(root3);
185 elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i);
186 elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii);
193 t = tp->nlat.s = tpoleinit[i][0]/root3;
194 tp->nlat.c = sqrt(1 - t*t);
195 tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c);
196 deg2rad(tpoleinit[i][1],&tp->wlon);
199 latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
200 deg2rad(tpp->ttwist,&tpp->projtw);
201 deg2rad(tpp->trot,&tpp->postrot);