1 28994509 2004-04-21 devnull #include <u.h>
2 28994509 2004-04-21 devnull #include <libc.h>
3 28994509 2004-04-21 devnull #include "map.h"
6 28994509 2004-04-21 devnull * conformal map of earth onto tetrahedron
7 28994509 2004-04-21 devnull * the stages of mapping are
8 28994509 2004-04-21 devnull * (a) stereo projection of tetrahedral face onto
9 28994509 2004-04-21 devnull * isosceles curvilinear triangle with 3 120-degree
10 28994509 2004-04-21 devnull * angles and one straight side
11 28994509 2004-04-21 devnull * (b) map of this triangle onto half plane cut along
12 28994509 2004-04-21 devnull * 3 rays from the roots of unity to infinity
13 28994509 2004-04-21 devnull * formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1)
14 28994509 2004-04-21 devnull * (c) do 3 times for each sector of plane:
15 28994509 2004-04-21 devnull * map of |arg z|<=pi/6, cut along z>1 into
16 28994509 2004-04-21 devnull * triangle |arg z|<=pi/6, Re z<=const,
17 28994509 2004-04-21 devnull * with upper side of cut going into upper half of
18 28994509 2004-04-21 devnull * of vertical side of triangle and lowere into lower
19 28994509 2004-04-21 devnull * formula int from 0 to z dz/sqrt(1-z^3)
21 28994509 2004-04-21 devnull * int from u to 1 3^.25*du/sqrt(1-u^3) =
22 28994509 2004-04-21 devnull F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4))
23 28994509 2004-04-21 devnull * int from 1 to u 3^.25*du/sqrt(u^3-1) =
24 28994509 2004-04-21 devnull * F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4))
25 28994509 2004-04-21 devnull * this latter formula extends analytically down to
26 28994509 2004-04-21 devnull * u=0 and is the basis of this routine, with the
27 28994509 2004-04-21 devnull * argument of complex elliptic integral elco2
28 28994509 2004-04-21 devnull * being tan(acos...)
29 28994509 2004-04-21 devnull * the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is
30 28994509 2004-04-21 devnull * used to cross over into the region where Re(acos...)>pi/2
31 28994509 2004-04-21 devnull * f0 and fpi are suitably scaled complete integrals
34 28994509 2004-04-21 devnull #define TFUZZ 0.00001
36 28994509 2004-04-21 devnull static struct place tpole[4]; /* point of tangency of tetrahedron face*/
37 28994509 2004-04-21 devnull static double tpoleinit[4][2] = {
39 28994509 2004-04-21 devnull 1., 180.,
40 28994509 2004-04-21 devnull -1., 90.,
41 28994509 2004-04-21 devnull -1., -90.
43 28994509 2004-04-21 devnull static struct tproj {
44 28994509 2004-04-21 devnull double tlat,tlon; /* center of stereo projection*/
45 28994509 2004-04-21 devnull double ttwist; /* rotatn before stereo*/
46 28994509 2004-04-21 devnull double trot; /*rotate after projection*/
47 28994509 2004-04-21 devnull struct place projpl; /*same as tlat,tlon*/
48 28994509 2004-04-21 devnull struct coord projtw; /*same as ttwist*/
49 28994509 2004-04-21 devnull struct coord postrot; /*same as trot*/
50 28994509 2004-04-21 devnull } tproj[4][4] = {
51 28994509 2004-04-21 devnull {/*00*/ {0.},
52 28994509 2004-04-21 devnull /*01*/ {90., 0., 90., -90.},
53 28994509 2004-04-21 devnull /*02*/ {0., 45., -45., 150.},
54 28994509 2004-04-21 devnull /*03*/ {0., -45., -135., 30.}
56 28994509 2004-04-21 devnull {/*10*/ {90., 0., -90., 90.},
57 28994509 2004-04-21 devnull /*11*/ {0.},
58 28994509 2004-04-21 devnull /*12*/ {0., 135., -135., -150.},
59 28994509 2004-04-21 devnull /*13*/ {0., -135., -45., -30.}
61 28994509 2004-04-21 devnull {/*20*/ {0., 45., 135., -30.},
62 28994509 2004-04-21 devnull /*21*/ {0., 135., 45., -150.},
63 28994509 2004-04-21 devnull /*22*/ {0.},
64 28994509 2004-04-21 devnull /*23*/ {-90., 0., 180., 90.}
66 28994509 2004-04-21 devnull {/*30*/ {0., -45., 45., -150.},
67 28994509 2004-04-21 devnull /*31*/ {0., -135., 135., -30.},
68 28994509 2004-04-21 devnull /*32*/ {-90., 0., 0., 90.},
69 28994509 2004-04-21 devnull /*33*/ {0.}
71 28994509 2004-04-21 devnull static double tx[4] = { /*where to move facet after final rotation*/
72 28994509 2004-04-21 devnull 0., 0., -1., 1. /*-1,1 to be sqrt(3)*/
74 28994509 2004-04-21 devnull static double ty[4] = {
75 28994509 2004-04-21 devnull 0., 2., -1., -1.
77 28994509 2004-04-21 devnull static double root3;
78 28994509 2004-04-21 devnull static double rt3inv;
79 28994509 2004-04-21 devnull static double two_rt3;
80 28994509 2004-04-21 devnull static double tkc,tk,tcon;
81 28994509 2004-04-21 devnull static double f0r,f0i,fpir,fpii;
83 28994509 2004-04-21 devnull static void
84 28994509 2004-04-21 devnull twhichp(struct place *g, int *p, int *q)
86 28994509 2004-04-21 devnull int i,j,k;
87 28994509 2004-04-21 devnull double cosdist[4];
88 28994509 2004-04-21 devnull struct place *tp;
89 28994509 2004-04-21 devnull for(i=0;i<4;i++) {
90 28994509 2004-04-21 devnull tp = &tpole[i];
91 28994509 2004-04-21 devnull cosdist[i] = g->nlat.s*tp->nlat.s +
92 28994509 2004-04-21 devnull g->nlat.c*tp->nlat.c*(
93 28994509 2004-04-21 devnull g->wlon.s*tp->wlon.s +
94 28994509 2004-04-21 devnull g->wlon.c*tp->wlon.c);
97 28994509 2004-04-21 devnull for(i=1;i<4;i++)
98 28994509 2004-04-21 devnull if(cosdist[i] > cosdist[j])
101 28994509 2004-04-21 devnull k = j==0?1:0;
102 28994509 2004-04-21 devnull for(i=0;i<4;i++)
103 28994509 2004-04-21 devnull if(i!=j&&cosdist[i]>cosdist[k])
109 28994509 2004-04-21 devnull Xtetra(struct place *place, double *x, double *y)
111 28994509 2004-04-21 devnull int i,j;
112 28994509 2004-04-21 devnull struct place pl;
113 28994509 2004-04-21 devnull register struct tproj *tpp;
114 28994509 2004-04-21 devnull double vr, vi;
115 28994509 2004-04-21 devnull double br, bi;
116 28994509 2004-04-21 devnull double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti;
117 28994509 2004-04-21 devnull twhichp(place,&i,&j);
118 28994509 2004-04-21 devnull copyplace(place,&pl);
119 28994509 2004-04-21 devnull norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw);
120 28994509 2004-04-21 devnull Xstereographic(&pl,&vr,&vi);
121 28994509 2004-04-21 devnull zr = vr/2;
122 28994509 2004-04-21 devnull zi = vi/2;
123 28994509 2004-04-21 devnull if(zr<=TFUZZ)
124 28994509 2004-04-21 devnull zr = TFUZZ;
125 28994509 2004-04-21 devnull csq(zr,zi,&z2r,&z2i);
126 28994509 2004-04-21 devnull csq(z2r,z2i,&z4r,&z4i);
127 28994509 2004-04-21 devnull z2r *= two_rt3;
128 28994509 2004-04-21 devnull z2i *= two_rt3;
129 28994509 2004-04-21 devnull cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si);
130 28994509 2004-04-21 devnull csqrt(sr-1,si,&tr,&ti);
131 28994509 2004-04-21 devnull cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi);
132 28994509 2004-04-21 devnull if(br<0) {
133 28994509 2004-04-21 devnull br = -br;
134 28994509 2004-04-21 devnull bi = -bi;
135 28994509 2004-04-21 devnull if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
136 28994509 2004-04-21 devnull return 0;
137 28994509 2004-04-21 devnull vr = fpir - vr;
138 28994509 2004-04-21 devnull vi = fpii - vi;
140 28994509 2004-04-21 devnull if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
141 28994509 2004-04-21 devnull return 0;
142 28994509 2004-04-21 devnull if(si>=0) {
143 28994509 2004-04-21 devnull tr = f0r - vi;
144 28994509 2004-04-21 devnull ti = f0i + vr;
145 28994509 2004-04-21 devnull } else {
146 28994509 2004-04-21 devnull tr = f0r + vi;
147 28994509 2004-04-21 devnull ti = f0i - vr;
149 28994509 2004-04-21 devnull tpp = &tproj[i][j];
150 28994509 2004-04-21 devnull *x = tr*tpp->postrot.c +
151 28994509 2004-04-21 devnull ti*tpp->postrot.s + tx[i];
152 28994509 2004-04-21 devnull *y = ti*tpp->postrot.c -
153 28994509 2004-04-21 devnull tr*tpp->postrot.s + ty[i];
154 28994509 2004-04-21 devnull return(1);
158 28994509 2004-04-21 devnull tetracut(struct place *g, struct place *og, double *cutlon)
160 28994509 2004-04-21 devnull int i,j,k;
161 28994509 2004-04-21 devnull if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) &&
162 28994509 2004-04-21 devnull (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
163 28994509 2004-04-21 devnull return(2);
164 28994509 2004-04-21 devnull twhichp(g,&i,&k);
165 28994509 2004-04-21 devnull twhichp(og,&j,&k);
166 28994509 2004-04-21 devnull if(i==j||i==0||j==0)
167 28994509 2004-04-21 devnull return(1);
168 28994509 2004-04-21 devnull return(0);
172 28994509 2004-04-21 devnull tetra(void)
176 28994509 2004-04-21 devnull register struct place *tp;
177 28994509 2004-04-21 devnull register struct tproj *tpp;
178 28994509 2004-04-21 devnull double t;
179 28994509 2004-04-21 devnull root3 = sqrt(3.);
180 28994509 2004-04-21 devnull rt3inv = 1/root3;
181 28994509 2004-04-21 devnull two_rt3 = 2*root3;
182 28994509 2004-04-21 devnull tkc = sqrt(.5-.25*root3);
183 28994509 2004-04-21 devnull tk = sqrt(.5+.25*root3);
184 28994509 2004-04-21 devnull tcon = 2*sqrt(root3);
185 28994509 2004-04-21 devnull elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i);
186 28994509 2004-04-21 devnull elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii);
187 28994509 2004-04-21 devnull fpir *= 2;
188 28994509 2004-04-21 devnull fpii *= 2;
189 28994509 2004-04-21 devnull for(i=0;i<4;i++) {
190 28994509 2004-04-21 devnull tx[i] *= f0r*root3;
191 28994509 2004-04-21 devnull ty[i] *= f0r;
192 28994509 2004-04-21 devnull tp = &tpole[i];
193 28994509 2004-04-21 devnull t = tp->nlat.s = tpoleinit[i][0]/root3;
194 28994509 2004-04-21 devnull tp->nlat.c = sqrt(1 - t*t);
195 28994509 2004-04-21 devnull tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c);
196 28994509 2004-04-21 devnull deg2rad(tpoleinit[i][1],&tp->wlon);
197 28994509 2004-04-21 devnull for(j=0;j<4;j++) {
198 28994509 2004-04-21 devnull tpp = &tproj[i][j];
199 28994509 2004-04-21 devnull latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
200 28994509 2004-04-21 devnull deg2rad(tpp->ttwist,&tpp->projtw);
201 28994509 2004-04-21 devnull deg2rad(tpp->trot,&tpp->postrot);
204 28994509 2004-04-21 devnull return(Xtetra);
