1 #include <u.h>
2 #include <libc.h>
3 #include "map.h"
4 
5 /*
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)
20  *
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
28  *	being tan(acos...)
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
32 */
33 
34 #define TFUZZ 0.00001
35 
36 static struct place tpole[4];	/* point of tangency of tetrahedron face*/
37 static double tpoleinit[4][2] = {
38 	1.,	0.,
39 	1.,	180.,
40 	-1.,	90.,
41 	-1.,	-90.
42 };
43 static struct tproj {
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*/
50 } tproj[4][4] = {
51 {/*00*/	{0.},
52  /*01*/	{90.,	0.,	90.,	-90.},
53  /*02*/	{0.,	45.,	-45.,	150.},
54  /*03*/	{0.,	-45.,	-135.,	30.}
55 },
56 {/*10*/	{90.,	0.,	-90.,	90.},
57  /*11*/ {0.},
58  /*12*/ {0.,	135.,	-135.,	-150.},
59  /*13*/	{0.,	-135.,	-45.,	-30.}
60 },
61 {/*20*/	{0.,	45.,	135.,	-30.},
62  /*21*/	{0.,	135.,	45.,	-150.},
63  /*22*/	{0.},
64  /*23*/	{-90.,	0.,	180.,	90.}
65 },
66 {/*30*/	{0.,	-45.,	45.,	-150.},
67  /*31*/ {0.,	-135.,	135.,	-30.},
68  /*32*/	{-90.,	0.,	0.,	90.},
69  /*33*/ {0.} 
70 }};
71 static double tx[4] = {	/*where to move facet after final rotation*/
72 	0.,	0.,	-1.,	1.	/*-1,1 to be sqrt(3)*/
73 };
74 static double ty[4] = {
75 	0.,	2.,	-1.,	-1.
76 };
77 static double root3;
78 static double rt3inv;
79 static double two_rt3;
80 static double tkc,tk,tcon;
81 static double f0r,f0i,fpir,fpii;
82 
83 static void
84 twhichp(struct place *g, int *p, int *q)
85 {
86 	int i,j,k;
87 	double cosdist[4];
88 	struct place *tp;
89 	for(i=0;i<4;i++) {
90 		tp = &tpole[i];
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);
95 	}
96 	j = 0;
97 	for(i=1;i<4;i++)
98 		if(cosdist[i] > cosdist[j])
99 			j = i;
100 	*p = j;
101 	k = j==0?1:0;
102 	for(i=0;i<4;i++)
103 		if(i!=j&&cosdist[i]>cosdist[k])
104 			k = i;
105 	*q = k;
106 }
107 
108 int
109 Xtetra(struct place *place, double *x, double *y)
110 {
111 	int i,j;
112 	struct place pl;
113 	register struct tproj *tpp;
114 	double vr, vi;
115 	double br, bi;
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);
121 	zr = vr/2;
122 	zi = vi/2;
123 	if(zr<=TFUZZ)
124 		zr = TFUZZ;
125 	csq(zr,zi,&z2r,&z2i);
126 	csq(z2r,z2i,&z4r,&z4i);
127 	z2r *= two_rt3;
128 	z2i *= two_rt3;
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);
132 	if(br<0) {
133 		br = -br;
134 		bi = -bi;
135 		if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
136 			return 0;
137 		vr = fpir - vr;
138 		vi = fpii - vi;
139 	} else 
140 		if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
141 			return 0;
142 	if(si>=0) {
143 		tr = f0r - vi;
144 		ti = f0i + vr;
145 	} else {
146 		tr = f0r + vi;
147 		ti = f0i - vr;
148 	}
149 	tpp = &tproj[i][j];
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];
154 	return(1);
155 }
156 
157 int
158 tetracut(struct place *g, struct place *og, double *cutlon)
159 {
160 	int i,j,k;
161 	if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) && 
162 	   (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
163 		return(2);
164 	twhichp(g,&i,&k);
165 	twhichp(og,&j,&k);
166 	if(i==j||i==0||j==0)
167 		return(1);
168 	return(0);
169 }
170 
171 proj
172 tetra(void)
173 {
174 	int i;
175 	int j;
176 	register struct place *tp;
177 	register struct tproj *tpp;
178 	double t;
179 	root3 = sqrt(3.);
180 	rt3inv = 1/root3;
181 	two_rt3 = 2*root3;
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);
187 	fpir *= 2;
188 	fpii *= 2;
189 	for(i=0;i<4;i++) {
190 		tx[i] *= f0r*root3;
191 		ty[i] *= f0r;
192 		tp = &tpole[i];
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);
197 		for(j=0;j<4;j++) {
198 			tpp = &tproj[i][j];
199 			latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
200 			deg2rad(tpp->ttwist,&tpp->projtw);
201 			deg2rad(tpp->trot,&tpp->postrot);
202 		}
203 	}
204 	return(Xtetra);
205 }
206