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