Logo Search packages:      
Sourcecode: r-cran-mapproj version File versions  Download package

tetra.c

/* RSB #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.
};
00043 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);
      map_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 int i; /* RSB */
      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);
}


Generated by  Doxygen 1.6.0   Back to index