mirror of
https://github.com/9fans/plan9port.git
synced 2025-01-15 11:20:03 +00:00
206 lines
4.8 KiB
C
206 lines
4.8 KiB
C
#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)
|
|
{
|
|
int 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);
|
|
}
|
|
|