svn-gvsig-desktop / tags / v2_0_0_Build_2027 / libraries / libjni-proj4 / src / PJ_laea.c @ 39128
History | View | Annotate | Download (5 KB)
| 1 |
#ifndef lint
|
|---|---|
| 2 |
static const char SCCSID[]="@(#)PJ_laea.c 4.1 94/02/15 GIE REL"; |
| 3 |
#endif
|
| 4 |
#define PROJ_PARMS__ \
|
| 5 |
double sinb1; \
|
| 6 |
double cosb1; \
|
| 7 |
double xmf; \
|
| 8 |
double ymf; \
|
| 9 |
double mmf; \
|
| 10 |
double qp; \
|
| 11 |
double dd; \
|
| 12 |
double rq; \
|
| 13 |
double *apa; \
|
| 14 |
int mode;
|
| 15 |
#define PJ_LIB__
|
| 16 |
#include <projects.h> |
| 17 |
PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell"; |
| 18 |
#define sinph0 P->sinb1
|
| 19 |
#define cosph0 P->cosb1
|
| 20 |
#define EPS10 1.e-10 |
| 21 |
#define NITER 20 |
| 22 |
#define CONV 1.e-10 |
| 23 |
#define N_POLE 0 |
| 24 |
#define S_POLE 1 |
| 25 |
#define EQUIT 2 |
| 26 |
#define OBLIQ 3 |
| 27 |
FORWARD(e_forward); /* ellipsoid */
|
| 28 |
double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0; |
| 29 |
|
| 30 |
coslam = cos(lp.lam); |
| 31 |
sinlam = sin(lp.lam); |
| 32 |
sinphi = sin(lp.phi); |
| 33 |
q = pj_qsfn(sinphi, P->e, P->one_es); |
| 34 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
| 35 |
sinb = q / P->qp; |
| 36 |
cosb = sqrt(1. - sinb * sinb);
|
| 37 |
} |
| 38 |
switch (P->mode) {
|
| 39 |
case OBLIQ:
|
| 40 |
b = 1. + P->sinb1 * sinb + P->cosb1 * cosb * coslam;
|
| 41 |
break;
|
| 42 |
case EQUIT:
|
| 43 |
b = 1. + cosb * coslam;
|
| 44 |
break;
|
| 45 |
case N_POLE:
|
| 46 |
b = HALFPI + lp.phi; |
| 47 |
q = P->qp - q; |
| 48 |
break;
|
| 49 |
case S_POLE:
|
| 50 |
b = lp.phi - HALFPI; |
| 51 |
q = P->qp + q; |
| 52 |
break;
|
| 53 |
} |
| 54 |
if (fabs(b) < EPS10) F_ERROR;
|
| 55 |
switch (P->mode) {
|
| 56 |
case OBLIQ:
|
| 57 |
xy.y = P->ymf * ( b = sqrt(2. / b) )
|
| 58 |
* (P->cosb1 * sinb - P->sinb1 * cosb * coslam); |
| 59 |
goto eqcon;
|
| 60 |
break;
|
| 61 |
case EQUIT:
|
| 62 |
xy.y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * P->ymf; |
| 63 |
eqcon:
|
| 64 |
xy.x = P->xmf * b * cosb * sinlam; |
| 65 |
break;
|
| 66 |
case N_POLE:
|
| 67 |
case S_POLE:
|
| 68 |
if (q >= 0.) { |
| 69 |
xy.x = (b = sqrt(q)) * sinlam; |
| 70 |
xy.y = coslam * (P->mode == S_POLE ? b : -b); |
| 71 |
} else
|
| 72 |
xy.x = xy.y = 0.;
|
| 73 |
break;
|
| 74 |
} |
| 75 |
return (xy);
|
| 76 |
} |
| 77 |
FORWARD(s_forward); /* spheroid */
|
| 78 |
double coslam, cosphi, sinphi;
|
| 79 |
|
| 80 |
sinphi = sin(lp.phi); |
| 81 |
cosphi = cos(lp.phi); |
| 82 |
coslam = cos(lp.lam); |
| 83 |
switch (P->mode) {
|
| 84 |
case EQUIT:
|
| 85 |
xy.y = 1. + cosphi * coslam;
|
| 86 |
goto oblcon;
|
| 87 |
case OBLIQ:
|
| 88 |
xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
|
| 89 |
oblcon:
|
| 90 |
if (xy.y <= EPS10) F_ERROR;
|
| 91 |
xy.x = (xy.y = sqrt(2. / xy.y)) * cosphi * sin(lp.lam);
|
| 92 |
xy.y *= P->mode == EQUIT ? sinphi : |
| 93 |
cosph0 * sinphi - sinph0 * cosphi * coslam; |
| 94 |
break;
|
| 95 |
case N_POLE:
|
| 96 |
coslam = -coslam; |
| 97 |
case S_POLE:
|
| 98 |
if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
|
| 99 |
xy.y = FORTPI - lp.phi * .5;
|
| 100 |
xy.y = 2. * (P->mode == S_POLE ? cos(xy.y) : sin(xy.y));
|
| 101 |
xy.x = xy.y * sin(lp.lam); |
| 102 |
xy.y *= coslam; |
| 103 |
break;
|
| 104 |
} |
| 105 |
return (xy);
|
| 106 |
} |
| 107 |
INVERSE(e_inverse); /* ellipsoid */
|
| 108 |
double cCe, sCe, q, rho, ab=0.0; |
| 109 |
|
| 110 |
switch (P->mode) {
|
| 111 |
case EQUIT:
|
| 112 |
case OBLIQ:
|
| 113 |
if ((rho = hypot(xy.x /= P->dd, xy.y *= P->dd)) < EPS10) {
|
| 114 |
lp.lam = 0.;
|
| 115 |
lp.phi = P->phi0; |
| 116 |
return (lp);
|
| 117 |
} |
| 118 |
cCe = cos(sCe = 2. * asin(.5 * rho / P->rq)); |
| 119 |
xy.x *= (sCe = sin(sCe)); |
| 120 |
if (P->mode == OBLIQ) {
|
| 121 |
q = P->qp * (ab = cCe * P->sinb1 + xy.y * sCe * P->cosb1 / rho); |
| 122 |
xy.y = rho * P->cosb1 * cCe - xy.y * P->sinb1 * sCe; |
| 123 |
} else {
|
| 124 |
q = P->qp * (ab = xy.y * sCe / rho); |
| 125 |
xy.y = rho * cCe; |
| 126 |
} |
| 127 |
break;
|
| 128 |
case N_POLE:
|
| 129 |
xy.y = -xy.y; |
| 130 |
case S_POLE:
|
| 131 |
if (!(q = (xy.x * xy.x + xy.y * xy.y)) ) {
|
| 132 |
lp.lam = 0.;
|
| 133 |
lp.phi = P->phi0; |
| 134 |
return (lp);
|
| 135 |
} |
| 136 |
/*
|
| 137 |
q = P->qp - q;
|
| 138 |
*/
|
| 139 |
ab = 1. - q / P->qp;
|
| 140 |
if (P->mode == S_POLE)
|
| 141 |
ab = - ab; |
| 142 |
break;
|
| 143 |
} |
| 144 |
lp.lam = atan2(xy.x, xy.y); |
| 145 |
lp.phi = pj_authlat(asin(ab), P->apa); |
| 146 |
return (lp);
|
| 147 |
} |
| 148 |
INVERSE(s_inverse); /* spheroid */
|
| 149 |
double cosz=0.0, rh, sinz=0.0; |
| 150 |
|
| 151 |
rh = hypot(xy.x, xy.y); |
| 152 |
if ((lp.phi = rh * .5 ) > 1.) I_ERROR; |
| 153 |
lp.phi = 2. * asin(lp.phi);
|
| 154 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
| 155 |
sinz = sin(lp.phi); |
| 156 |
cosz = cos(lp.phi); |
| 157 |
} |
| 158 |
switch (P->mode) {
|
| 159 |
case EQUIT:
|
| 160 |
lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
|
| 161 |
xy.x *= sinz; |
| 162 |
xy.y = cosz * rh; |
| 163 |
break;
|
| 164 |
case OBLIQ:
|
| 165 |
lp.phi = fabs(rh) <= EPS10 ? P->phi0 : |
| 166 |
asin(cosz * sinph0 + xy.y * sinz * cosph0 / rh); |
| 167 |
xy.x *= sinz * cosph0; |
| 168 |
xy.y = (cosz - sin(lp.phi) * sinph0) * rh; |
| 169 |
break;
|
| 170 |
case N_POLE:
|
| 171 |
xy.y = -xy.y; |
| 172 |
lp.phi = HALFPI - lp.phi; |
| 173 |
break;
|
| 174 |
case S_POLE:
|
| 175 |
lp.phi -= HALFPI; |
| 176 |
break;
|
| 177 |
} |
| 178 |
lp.lam = (xy.y == 0. && (P->mode == EQUIT || P->mode == OBLIQ)) ?
|
| 179 |
0. : atan2(xy.x, xy.y);
|
| 180 |
return (lp);
|
| 181 |
} |
| 182 |
FREEUP; if (P) pj_dalloc(P); }
|
| 183 |
ENTRY0(laea) |
| 184 |
double t;
|
| 185 |
|
| 186 |
if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
|
| 187 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
| 188 |
else if (fabs(t) < EPS10) |
| 189 |
P->mode = EQUIT; |
| 190 |
else
|
| 191 |
P->mode = OBLIQ; |
| 192 |
if (P->es) {
|
| 193 |
double sinphi;
|
| 194 |
|
| 195 |
P->e = sqrt(P->es); |
| 196 |
P->qp = pj_qsfn(1., P->e, P->one_es);
|
| 197 |
P->mmf = .5 / (1. - P->es); |
| 198 |
P->apa = pj_authset(P->es); |
| 199 |
switch (P->mode) {
|
| 200 |
case N_POLE:
|
| 201 |
case S_POLE:
|
| 202 |
P->dd = 1.;
|
| 203 |
break;
|
| 204 |
case EQUIT:
|
| 205 |
P->dd = 1. / (P->rq = sqrt(.5 * P->qp)); |
| 206 |
P->xmf = 1.;
|
| 207 |
P->ymf = .5 * P->qp;
|
| 208 |
break;
|
| 209 |
case OBLIQ:
|
| 210 |
P->rq = sqrt(.5 * P->qp);
|
| 211 |
sinphi = sin(P->phi0); |
| 212 |
P->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / P->qp; |
| 213 |
P->cosb1 = sqrt(1. - P->sinb1 * P->sinb1);
|
| 214 |
P->dd = cos(P->phi0) / (sqrt(1. - P->es * sinphi * sinphi) *
|
| 215 |
P->rq * P->cosb1); |
| 216 |
P->ymf = (P->xmf = P->rq) / P->dd; |
| 217 |
P->xmf *= P->dd; |
| 218 |
break;
|
| 219 |
} |
| 220 |
P->inv = e_inverse; |
| 221 |
P->fwd = e_forward; |
| 222 |
} else {
|
| 223 |
if (P->mode == OBLIQ) {
|
| 224 |
sinph0 = sin(P->phi0); |
| 225 |
cosph0 = cos(P->phi0); |
| 226 |
} |
| 227 |
P->inv = s_inverse; |
| 228 |
P->fwd = s_forward; |
| 229 |
} |
| 230 |
ENDENTRY(P) |