svn-gvsig-desktop / trunk / libraries / libjni-proj4 / src / PJ_stere.c @ 7098
History | View | Annotate | Download (5.74 KB)
| 1 |
#ifndef lint
|
|---|---|
| 2 |
static const char SCCSID[]="@(#)PJ_stere.c 4.1 94/02/15 GIE REL"; |
| 3 |
#endif
|
| 4 |
#define PROJ_PARMS__ \
|
| 5 |
double phits; \
|
| 6 |
double sinX1; \
|
| 7 |
double cosX1; \
|
| 8 |
double akm1; \
|
| 9 |
int mode;
|
| 10 |
#define PJ_LIB__
|
| 11 |
#include <projects.h> |
| 12 |
PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts="; |
| 13 |
PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth"; |
| 14 |
#define sinph0 P->sinX1
|
| 15 |
#define cosph0 P->cosX1
|
| 16 |
#define EPS10 1.e-10 |
| 17 |
#define TOL 1.e-8 |
| 18 |
#define NITER 8 |
| 19 |
#define CONV 1.e-10 |
| 20 |
#define S_POLE 0 |
| 21 |
#define N_POLE 1 |
| 22 |
#define OBLIQ 2 |
| 23 |
#define EQUIT 3 |
| 24 |
static double |
| 25 |
ssfn_(double phit, double sinphi, double eccen) { |
| 26 |
sinphi *= eccen; |
| 27 |
return (tan (.5 * (HALFPI + phit)) * |
| 28 |
pow((1. - sinphi) / (1. + sinphi), .5 * eccen)); |
| 29 |
} |
| 30 |
FORWARD(e_forward); /* ellipsoid */
|
| 31 |
double coslam, sinlam, sinX=0.0, cosX=0.0, X, A, sinphi; |
| 32 |
|
| 33 |
coslam = cos(lp.lam); |
| 34 |
sinlam = sin(lp.lam); |
| 35 |
sinphi = sin(lp.phi); |
| 36 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
| 37 |
sinX = sin(X = 2. * atan(ssfn_(lp.phi, sinphi, P->e)) - HALFPI);
|
| 38 |
cosX = cos(X); |
| 39 |
} |
| 40 |
switch (P->mode) {
|
| 41 |
case OBLIQ:
|
| 42 |
A = P->akm1 / (P->cosX1 * (1. + P->sinX1 * sinX +
|
| 43 |
P->cosX1 * cosX * coslam)); |
| 44 |
xy.y = A * (P->cosX1 * sinX - P->sinX1 * cosX * coslam); |
| 45 |
goto xmul;
|
| 46 |
case EQUIT:
|
| 47 |
A = 2. * P->akm1 / (1. + cosX * coslam); |
| 48 |
xy.y = A * sinX; |
| 49 |
xmul:
|
| 50 |
xy.x = A * cosX; |
| 51 |
break;
|
| 52 |
case S_POLE:
|
| 53 |
lp.phi = -lp.phi; |
| 54 |
coslam = - coslam; |
| 55 |
sinphi = -sinphi; |
| 56 |
case N_POLE:
|
| 57 |
xy.x = P->akm1 * pj_tsfn(lp.phi, sinphi, P->e); |
| 58 |
xy.y = - xy.x * coslam; |
| 59 |
break;
|
| 60 |
} |
| 61 |
xy.x = xy.x * sinlam; |
| 62 |
return (xy);
|
| 63 |
} |
| 64 |
FORWARD(s_forward); /* spheroid */
|
| 65 |
double sinphi, cosphi, coslam, sinlam;
|
| 66 |
|
| 67 |
sinphi = sin(lp.phi); |
| 68 |
cosphi = cos(lp.phi); |
| 69 |
coslam = cos(lp.lam); |
| 70 |
sinlam = sin(lp.lam); |
| 71 |
switch (P->mode) {
|
| 72 |
case EQUIT:
|
| 73 |
xy.y = 1. + cosphi * coslam;
|
| 74 |
goto oblcon;
|
| 75 |
case OBLIQ:
|
| 76 |
xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
|
| 77 |
oblcon:
|
| 78 |
if (xy.y <= EPS10) F_ERROR;
|
| 79 |
xy.x = (xy.y = P->akm1 / xy.y) * cosphi * sinlam; |
| 80 |
xy.y *= (P->mode == EQUIT) ? sinphi : |
| 81 |
cosph0 * sinphi - sinph0 * cosphi * coslam; |
| 82 |
break;
|
| 83 |
case N_POLE:
|
| 84 |
coslam = - coslam; |
| 85 |
lp.phi = - lp.phi; |
| 86 |
case S_POLE:
|
| 87 |
if (fabs(lp.phi - HALFPI) < TOL) F_ERROR;
|
| 88 |
xy.x = sinlam * ( xy.y = P->akm1 * tan(FORTPI + .5 * lp.phi) );
|
| 89 |
xy.y *= coslam; |
| 90 |
break;
|
| 91 |
} |
| 92 |
return (xy);
|
| 93 |
} |
| 94 |
INVERSE(e_inverse); /* ellipsoid */
|
| 95 |
double cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0; |
| 96 |
int i;
|
| 97 |
|
| 98 |
rho = hypot(xy.x, xy.y); |
| 99 |
switch (P->mode) {
|
| 100 |
case OBLIQ:
|
| 101 |
case EQUIT:
|
| 102 |
cosphi = cos( tp = 2. * atan2(rho * P->cosX1 , P->akm1) );
|
| 103 |
sinphi = sin(tp); |
| 104 |
if( rho == 0.0 ) |
| 105 |
phi_l = asin(cosphi * P->sinX1); |
| 106 |
else
|
| 107 |
phi_l = asin(cosphi * P->sinX1 + (xy.y * sinphi * P->cosX1 / rho)); |
| 108 |
|
| 109 |
tp = tan(.5 * (HALFPI + phi_l));
|
| 110 |
xy.x *= sinphi; |
| 111 |
xy.y = rho * P->cosX1 * cosphi - xy.y * P->sinX1* sinphi; |
| 112 |
halfpi = HALFPI; |
| 113 |
halfe = .5 * P->e;
|
| 114 |
break;
|
| 115 |
case N_POLE:
|
| 116 |
xy.y = -xy.y; |
| 117 |
case S_POLE:
|
| 118 |
phi_l = HALFPI - 2. * atan(tp = - rho / P->akm1);
|
| 119 |
halfpi = -HALFPI; |
| 120 |
halfe = -.5 * P->e;
|
| 121 |
break;
|
| 122 |
} |
| 123 |
for (i = NITER; i--; phi_l = lp.phi) {
|
| 124 |
sinphi = P->e * sin(phi_l); |
| 125 |
lp.phi = 2. * atan(tp * pow((1.+sinphi)/(1.-sinphi), |
| 126 |
halfe)) - halfpi; |
| 127 |
if (fabs(phi_l - lp.phi) < CONV) {
|
| 128 |
if (P->mode == S_POLE)
|
| 129 |
lp.phi = -lp.phi; |
| 130 |
lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y); |
| 131 |
return (lp);
|
| 132 |
} |
| 133 |
} |
| 134 |
I_ERROR; |
| 135 |
} |
| 136 |
INVERSE(s_inverse); /* spheroid */
|
| 137 |
double c, rh, sinc, cosc;
|
| 138 |
|
| 139 |
sinc = sin(c = 2. * atan((rh = hypot(xy.x, xy.y)) / P->akm1));
|
| 140 |
cosc = cos(c); |
| 141 |
lp.lam = 0.;
|
| 142 |
switch (P->mode) {
|
| 143 |
case EQUIT:
|
| 144 |
if (fabs(rh) <= EPS10)
|
| 145 |
lp.phi = 0.;
|
| 146 |
else
|
| 147 |
lp.phi = asin(xy.y * sinc / rh); |
| 148 |
if (cosc != 0. || xy.x != 0.) |
| 149 |
lp.lam = atan2(xy.x * sinc, cosc * rh); |
| 150 |
break;
|
| 151 |
case OBLIQ:
|
| 152 |
if (fabs(rh) <= EPS10)
|
| 153 |
lp.phi = P->phi0; |
| 154 |
else
|
| 155 |
lp.phi = asin(cosc * sinph0 + xy.y * sinc * cosph0 / rh); |
| 156 |
if ((c = cosc - sinph0 * sin(lp.phi)) != 0. || xy.x != 0.) |
| 157 |
lp.lam = atan2(xy.x * sinc * cosph0, c * rh); |
| 158 |
break;
|
| 159 |
case N_POLE:
|
| 160 |
xy.y = -xy.y; |
| 161 |
case S_POLE:
|
| 162 |
if (fabs(rh) <= EPS10)
|
| 163 |
lp.phi = P->phi0; |
| 164 |
else
|
| 165 |
lp.phi = asin(P->mode == S_POLE ? - cosc : cosc); |
| 166 |
lp.lam = (xy.x == 0. && xy.y == 0.) ? 0. : atan2(xy.x, xy.y); |
| 167 |
break;
|
| 168 |
} |
| 169 |
return (lp);
|
| 170 |
} |
| 171 |
FREEUP; if (P) pj_dalloc(P); }
|
| 172 |
static PJ *
|
| 173 |
setup(PJ *P) { /* general initialization */
|
| 174 |
double t;
|
| 175 |
|
| 176 |
if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
|
| 177 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
| 178 |
else
|
| 179 |
P->mode = t > EPS10 ? OBLIQ : EQUIT; |
| 180 |
P->phits = fabs(P->phits); |
| 181 |
if (P->es) {
|
| 182 |
double X;
|
| 183 |
|
| 184 |
switch (P->mode) {
|
| 185 |
case N_POLE:
|
| 186 |
case S_POLE:
|
| 187 |
if (fabs(P->phits - HALFPI) < EPS10)
|
| 188 |
P->akm1 = 2. * P->k0 /
|
| 189 |
sqrt(pow(1+P->e,1+P->e)*pow(1-P->e,1-P->e)); |
| 190 |
else {
|
| 191 |
P->akm1 = cos(P->phits) / |
| 192 |
pj_tsfn(P->phits, t = sin(P->phits), P->e); |
| 193 |
t *= P->e; |
| 194 |
P->akm1 /= sqrt(1. - t * t);
|
| 195 |
} |
| 196 |
break;
|
| 197 |
case EQUIT:
|
| 198 |
P->akm1 = 2. * P->k0;
|
| 199 |
break;
|
| 200 |
case OBLIQ:
|
| 201 |
t = sin(P->phi0); |
| 202 |
X = 2. * atan(ssfn_(P->phi0, t, P->e)) - HALFPI;
|
| 203 |
t *= P->e; |
| 204 |
P->akm1 = 2. * P->k0 * cos(P->phi0) / sqrt(1. - t * t); |
| 205 |
P->sinX1 = sin(X); |
| 206 |
P->cosX1 = cos(X); |
| 207 |
break;
|
| 208 |
} |
| 209 |
P->inv = e_inverse; |
| 210 |
P->fwd = e_forward; |
| 211 |
} else {
|
| 212 |
switch (P->mode) {
|
| 213 |
case OBLIQ:
|
| 214 |
sinph0 = sin(P->phi0); |
| 215 |
cosph0 = cos(P->phi0); |
| 216 |
case EQUIT:
|
| 217 |
P->akm1 = 2. * P->k0;
|
| 218 |
break;
|
| 219 |
case S_POLE:
|
| 220 |
case N_POLE:
|
| 221 |
P->akm1 = fabs(P->phits - HALFPI) >= EPS10 ? |
| 222 |
cos(P->phits) / tan(FORTPI - .5 * P->phits) :
|
| 223 |
2. * P->k0 ;
|
| 224 |
break;
|
| 225 |
} |
| 226 |
P->inv = s_inverse; |
| 227 |
P->fwd = s_forward; |
| 228 |
} |
| 229 |
return P;
|
| 230 |
} |
| 231 |
ENTRY0(stere) |
| 232 |
P->phits = pj_param(P->params, "tlat_ts").i ?
|
| 233 |
P->phits = pj_param(P->params, "rlat_ts").f : HALFPI;
|
| 234 |
ENDENTRY(setup(P)) |
| 235 |
ENTRY0(ups) |
| 236 |
/* International Ellipsoid */
|
| 237 |
P->phi0 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI;
|
| 238 |
if (!P->es) E_ERROR(-34); |
| 239 |
P->k0 = .994;
|
| 240 |
P->x0 = 2000000.; |
| 241 |
P->y0 = 2000000.; |
| 242 |
P->phits = HALFPI; |
| 243 |
P->lam0 = 0.;
|
| 244 |
ENDENTRY(setup(P)) |