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