123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206 |
- #include <u.h>
- #include <libc.h>
- #include "map.h"
- /*
- * conformal map of earth onto tetrahedron
- * the stages of mapping are
- * (a) stereo projection of tetrahedral face onto
- * isosceles curvilinear triangle with 3 120-degree
- * angles and one straight side
- * (b) map of this triangle onto half plane cut along
- * 3 rays from the roots of unity to infinity
- * formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1)
- * (c) do 3 times for each sector of plane:
- * map of |arg z|<=pi/6, cut along z>1 into
- * triangle |arg z|<=pi/6, Re z<=const,
- * with upper side of cut going into upper half of
- * of vertical side of triangle and lowere into lower
- * formula int from 0 to z dz/sqrt(1-z^3)
- *
- * int from u to 1 3^.25*du/sqrt(1-u^3) =
- F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4))
- * int from 1 to u 3^.25*du/sqrt(u^3-1) =
- * F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4))
- * this latter formula extends analytically down to
- * u=0 and is the basis of this routine, with the
- * argument of complex elliptic integral elco2
- * being tan(acos...)
- * the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is
- * used to cross over into the region where Re(acos...)>pi/2
- * f0 and fpi are suitably scaled complete integrals
- */
- #define TFUZZ 0.00001
- static struct place tpole[4]; /* point of tangency of tetrahedron face*/
- static double tpoleinit[4][2] = {
- 1., 0.,
- 1., 180.,
- -1., 90.,
- -1., -90.
- };
- static struct tproj {
- double tlat,tlon; /* center of stereo projection*/
- double ttwist; /* rotatn before stereo*/
- double trot; /*rotate after projection*/
- struct place projpl; /*same as tlat,tlon*/
- struct coord projtw; /*same as ttwist*/
- struct coord postrot; /*same as trot*/
- } tproj[4][4] = {
- {/*00*/ {0.},
- /*01*/ {90., 0., 90., -90.},
- /*02*/ {0., 45., -45., 150.},
- /*03*/ {0., -45., -135., 30.}
- },
- {/*10*/ {90., 0., -90., 90.},
- /*11*/ {0.},
- /*12*/ {0., 135., -135., -150.},
- /*13*/ {0., -135., -45., -30.}
- },
- {/*20*/ {0., 45., 135., -30.},
- /*21*/ {0., 135., 45., -150.},
- /*22*/ {0.},
- /*23*/ {-90., 0., 180., 90.}
- },
- {/*30*/ {0., -45., 45., -150.},
- /*31*/ {0., -135., 135., -30.},
- /*32*/ {-90., 0., 0., 90.},
- /*33*/ {0.}
- }};
- static double tx[4] = { /*where to move facet after final rotation*/
- 0., 0., -1., 1. /*-1,1 to be sqrt(3)*/
- };
- static double ty[4] = {
- 0., 2., -1., -1.
- };
- static double root3;
- static double rt3inv;
- static double two_rt3;
- static double tkc,tk,tcon;
- static double f0r,f0i,fpir,fpii;
- static void
- twhichp(struct place *g, int *p, int *q)
- {
- int i,j,k;
- double cosdist[4];
- struct place *tp;
- for(i=0;i<4;i++) {
- tp = &tpole[i];
- cosdist[i] = g->nlat.s*tp->nlat.s +
- g->nlat.c*tp->nlat.c*(
- g->wlon.s*tp->wlon.s +
- g->wlon.c*tp->wlon.c);
- }
- j = 0;
- for(i=1;i<4;i++)
- if(cosdist[i] > cosdist[j])
- j = i;
- *p = j;
- k = j==0?1:0;
- for(i=0;i<4;i++)
- if(i!=j&&cosdist[i]>cosdist[k])
- k = i;
- *q = k;
- }
- int
- Xtetra(struct place *place, double *x, double *y)
- {
- int i,j;
- struct place pl;
- register struct tproj *tpp;
- double vr, vi;
- double br, bi;
- double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti;
- twhichp(place,&i,&j);
- copyplace(place,&pl);
- norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw);
- Xstereographic(&pl,&vr,&vi);
- zr = vr/2;
- zi = vi/2;
- if(zr<=TFUZZ)
- zr = TFUZZ;
- csq(zr,zi,&z2r,&z2i);
- csq(z2r,z2i,&z4r,&z4i);
- z2r *= two_rt3;
- z2i *= two_rt3;
- cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si);
- csqrt(sr-1,si,&tr,&ti);
- cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi);
- if(br<0) {
- br = -br;
- bi = -bi;
- if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
- return 0;
- vr = fpir - vr;
- vi = fpii - vi;
- } else
- if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
- return 0;
- if(si>=0) {
- tr = f0r - vi;
- ti = f0i + vr;
- } else {
- tr = f0r + vi;
- ti = f0i - vr;
- }
- tpp = &tproj[i][j];
- *x = tr*tpp->postrot.c +
- ti*tpp->postrot.s + tx[i];
- *y = ti*tpp->postrot.c -
- tr*tpp->postrot.s + ty[i];
- return(1);
- }
- int
- tetracut(struct place *g, struct place *og, double *cutlon)
- {
- int i,j,k;
- if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) &&
- (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
- return(2);
- twhichp(g,&i,&k);
- twhichp(og,&j,&k);
- if(i==j||i==0||j==0)
- return(1);
- return(0);
- }
- proj
- tetra(void)
- {
- register i;
- int j;
- register struct place *tp;
- register struct tproj *tpp;
- double t;
- root3 = sqrt(3.);
- rt3inv = 1/root3;
- two_rt3 = 2*root3;
- tkc = sqrt(.5-.25*root3);
- tk = sqrt(.5+.25*root3);
- tcon = 2*sqrt(root3);
- elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i);
- elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii);
- fpir *= 2;
- fpii *= 2;
- for(i=0;i<4;i++) {
- tx[i] *= f0r*root3;
- ty[i] *= f0r;
- tp = &tpole[i];
- t = tp->nlat.s = tpoleinit[i][0]/root3;
- tp->nlat.c = sqrt(1 - t*t);
- tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c);
- deg2rad(tpoleinit[i][1],&tp->wlon);
- for(j=0;j<4;j++) {
- tpp = &tproj[i][j];
- latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
- deg2rad(tpp->ttwist,&tpp->projtw);
- deg2rad(tpp->trot,&tpp->postrot);
- }
- }
- return(Xtetra);
- }
|