tetra.c (4935B)
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 }