aboutsummaryrefslogtreecommitdiff
path: root/src/PJ_aea.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/PJ_aea.c')
-rw-r--r--src/PJ_aea.c316
1 files changed, 247 insertions, 69 deletions
diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 0881ba78..8285e68f 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -1,7 +1,9 @@
/******************************************************************************
* Project: PROJ.4
* Purpose: Implementation of the aea (Albers Equal Area) projection.
- * Author: Gerald Evenden
+ * and the leac (Lambert Equal Area Conic) projection
+ * Author: Gerald Evenden (1995)
+ * Thomas Knudsen (2016) - revise/add regression tests
*
******************************************************************************
* Copyright (c) 1995, Gerald Evenden
@@ -25,35 +27,21 @@
* DEALINGS IN THE SOFTWARE.
*****************************************************************************/
-#define PROJ_PARMS__ \
- double ec; \
- double n; \
- double c; \
- double dd; \
- double n2; \
- double rho0; \
- double rho; \
- double phi1; \
- double phi2; \
- double *en; \
- int ellips;
-
#define PJ_LIB__
#include <projects.h>
# define EPS10 1.e-10
# define TOL7 1.e-7
-PROJ_HEAD(aea, "Albers Equal Area")
- "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
-PROJ_HEAD(leac, "Lambert Equal Area Conic")
- "\n\tConic, Sph&Ell\n\tlat_1= south";
+PROJ_HEAD(aea, "Albers Equal Area") "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
+PROJ_HEAD(leac, "Lambert Equal Area Conic") "\n\tConic, Sph&Ell\n\tlat_1= south";
+
+
/* determine latitude angle phi-1 */
# define N_ITER 15
# define EPSILON 1.0e-7
# define TOL 1.0e-10
- static double
-phi1_(double qs, double Te, double Tone_es) {
+static double phi1_(double qs, double Te, double Tone_es) {
int i;
double Phi, sinpi, cospi, con, com, dphi;
@@ -73,86 +61,276 @@ phi1_(double qs, double Te, double Tone_es) {
} while (fabs(dphi) > TOL && --i);
return( i ? Phi : HUGE_VAL );
}
-FORWARD(e_forward); /* ellipsoid & spheroid */
- if ((P->rho = P->c - (P->ellips ? P->n * pj_qsfn(sin(lp.phi),
- P->e, P->one_es) : P->n2 * sin(lp.phi))) < 0.) F_ERROR
- P->rho = P->dd * sqrt(P->rho);
- xy.x = P->rho * sin( lp.lam *= P->n );
- xy.y = P->rho0 - P->rho * cos(lp.lam);
- return (xy);
+
+
+struct pj_opaque {
+ double ec;
+ double n;
+ double c;
+ double dd;
+ double n2;
+ double rho0;
+ double rho;
+ double phi1;
+ double phi2;
+ double *en;
+ int ellips;
+};
+
+
+
+static XY e_forward (LP lp, PJ *P) { /* Ellipsoid/spheroid, forward */
+ XY xy = {0.0,0.0};
+ struct pj_opaque *Q = P->opaque;
+ if ((Q->rho = Q->c - (Q->ellips ? Q->n * pj_qsfn(sin(lp.phi),
+ P->e, P->one_es) : Q->n2 * sin(lp.phi))) < 0.) F_ERROR
+ Q->rho = Q->dd * sqrt(Q->rho);
+ xy.x = Q->rho * sin( lp.lam *= Q->n );
+ xy.y = Q->rho0 - Q->rho * cos(lp.lam);
+ return xy;
}
-INVERSE(e_inverse) /* ellipsoid & spheroid */;
- if( (P->rho = hypot(xy.x, xy.y = P->rho0 - xy.y)) != 0.0 ) {
- if (P->n < 0.) {
- P->rho = -P->rho;
+
+
+static LP e_inverse (XY xy, PJ *P) { /* Ellipsoid/spheroid, inverse */
+ LP lp = {0.0,0.0};
+ struct pj_opaque *Q = P->opaque;
+ if( (Q->rho = hypot(xy.x, xy.y = Q->rho0 - xy.y)) != 0.0 ) {
+ if (Q->n < 0.) {
+ Q->rho = -Q->rho;
xy.x = -xy.x;
xy.y = -xy.y;
}
- lp.phi = P->rho / P->dd;
- if (P->ellips) {
- lp.phi = (P->c - lp.phi * lp.phi) / P->n;
- if (fabs(P->ec - fabs(lp.phi)) > TOL7) {
+ lp.phi = Q->rho / Q->dd;
+ if (Q->ellips) {
+ lp.phi = (Q->c - lp.phi * lp.phi) / Q->n;
+ if (fabs(Q->ec - fabs(lp.phi)) > TOL7) {
if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
I_ERROR
} else
lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
- } else if (fabs(lp.phi = (P->c - lp.phi * lp.phi) / P->n2) <= 1.)
+ } else if (fabs(lp.phi = (Q->c - lp.phi * lp.phi) / Q->n2) <= 1.)
lp.phi = asin(lp.phi);
else
lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
- lp.lam = atan2(xy.x, xy.y) / P->n;
+ lp.lam = atan2(xy.x, xy.y) / Q->n;
} else {
lp.lam = 0.;
- lp.phi = P->n > 0. ? HALFPI : - HALFPI;
+ lp.phi = Q->n > 0. ? HALFPI : - HALFPI;
}
- return (lp);
+ return lp;
}
-FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
- static PJ *
-setup(PJ *P) {
+
+
+static void *freeup_new (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+
+ if (0==P->opaque)
+ return pj_dealloc (P);
+
+ pj_dealloc (P->opaque->en);
+ pj_dealloc (P->opaque);
+ return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+ freeup_new (P);
+ return;
+}
+
+
+static PJ *setup(PJ *P) {
double cosphi, sinphi;
int secant;
+ struct pj_opaque *Q = P->opaque;
- if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21);
- P->n = sinphi = sin(P->phi1);
- cosphi = cos(P->phi1);
- secant = fabs(P->phi1 - P->phi2) >= EPS10;
- if( (P->ellips = (P->es > 0.))) {
+ P->inv = e_inverse;
+ P->fwd = e_forward;
+
+ if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21);
+ Q->n = sinphi = sin(Q->phi1);
+ cosphi = cos(Q->phi1);
+ secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
+ if( (Q->ellips = (P->es > 0.))) {
double ml1, m1;
- if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
+ if (!(Q->en = pj_enfn(P->es))) E_ERROR_0;
m1 = pj_msfn(sinphi, cosphi, P->es);
ml1 = pj_qsfn(sinphi, P->e, P->one_es);
if (secant) { /* secant cone */
double ml2, m2;
- sinphi = sin(P->phi2);
- cosphi = cos(P->phi2);
+ sinphi = sin(Q->phi2);
+ cosphi = cos(Q->phi2);
m2 = pj_msfn(sinphi, cosphi, P->es);
ml2 = pj_qsfn(sinphi, P->e, P->one_es);
- P->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
+ Q->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
}
- P->ec = 1. - .5 * P->one_es * log((1. - P->e) /
+ Q->ec = 1. - .5 * P->one_es * log((1. - P->e) /
(1. + P->e)) / P->e;
- P->c = m1 * m1 + P->n * ml1;
- P->dd = 1. / P->n;
- P->rho0 = P->dd * sqrt(P->c - P->n * pj_qsfn(sin(P->phi0),
+ Q->c = m1 * m1 + Q->n * ml1;
+ Q->dd = 1. / Q->n;
+ Q->rho0 = Q->dd * sqrt(Q->c - Q->n * pj_qsfn(sin(P->phi0),
P->e, P->one_es));
} else {
- if (secant) P->n = .5 * (P->n + sin(P->phi2));
- P->n2 = P->n + P->n;
- P->c = cosphi * cosphi + P->n2 * sinphi;
- P->dd = 1. / P->n;
- P->rho0 = P->dd * sqrt(P->c - P->n2 * sin(P->phi0));
+ if (secant) Q->n = .5 * (Q->n + sin(Q->phi2));
+ Q->n2 = Q->n + Q->n;
+ Q->c = cosphi * cosphi + Q->n2 * sinphi;
+ Q->dd = 1. / Q->n;
+ Q->rho0 = Q->dd * sqrt(Q->c - Q->n2 * sin(P->phi0));
}
- P->inv = e_inverse; P->fwd = e_forward;
+
return P;
}
-ENTRY1(aea,en)
- P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
- P->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
-ENDENTRY(setup(P))
-ENTRY1(leac,en)
- P->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
- P->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? - HALFPI: HALFPI;
-ENDENTRY(setup(P))
+
+
+PJ *PROJECTION(aea) {
+ struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ if (0==Q)
+ return freeup_new (P);
+ P->opaque = Q;
+
+ P->pfree = freeup;
+ P->descr = des_aea;
+
+ Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+ Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
+ setup(P);
+ return P;
+}
+
+
+PJ *PROJECTION(leac) {
+ struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ if (0==Q)
+ return freeup_new (P);
+ P->opaque = Q;
+
+ P->pfree = freeup;
+ P->descr = des_leac;
+ Q->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
+ Q->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? - HALFPI: HALFPI;
+ setup(P);
+ return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_aea_selftest (void) {return 0;}
+#else
+
+int pj_aea_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char e_args[] = {"+proj=aea +ellps=GRS80 +lat_1=0 +lat_2=2"};
+ char s_args[] = {"+proj=aea +ellps=GRS80 +lat_1=0 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY e_fwd_expect[] = {
+ {222571.60875710563, 110653.32674302977},
+ {222706.30650839131, -110484.26714439997},
+ {-222571.60875710563, 110653.32674302977},
+ {-222706.30650839131, -110484.26714439997},
+ };
+
+ XY s_fwd_expect[] = {
+ {223334.08517088494, 111780.43188447191},
+ {223470.15499168713, -111610.33943099028},
+ {-223334.08517088494, 111780.43188447191},
+ {-223470.15499168713, -111610.33943099028},
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP e_inv_expect[] = {
+ {0.0017966310597749514, 0.00090436885862202158},
+ {0.0017966300767030448, -0.00090437009538581453},
+ {-0.0017966310597749514, 0.00090436885862202158},
+ {-0.0017966300767030448, -0.00090437009538581453},
+ };
+
+ LP s_inv_expect[] = {
+ {0.0017904935979658752, 0.00089524594491375306},
+ {0.0017904926216016812, -0.00089524716502493225},
+ {-0.0017904935979658752, 0.00089524594491375306},
+ {-0.0017904926216016812, -0.00089524716502493225},
+ };
+ return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
+}
+
+
+#endif
+
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_leac_selftest (void) {return 0;}
+#else
+
+int pj_leac_selftest (void) {
+ double tolerance_lp = 1e-10;
+ double tolerance_xy = 1e-7;
+
+ char e_args[] = {"+proj=leac +ellps=GRS80 +lat_1=0 +lat_2=2"};
+ char s_args[] = {"+proj=leac +ellps=GRS80 +lat_1=0 +lat_2=2"};
+
+ LP fwd_in[] = {
+ { 2, 1},
+ { 2,-1},
+ {-2, 1},
+ {-2,-1}
+ };
+
+ XY e_fwd_expect[] = {
+ {220685.14054297868, 112983.50088939646},
+ {224553.31227982609, -108128.63674487274},
+ {-220685.14054297868, 112983.50088939646},
+ {-224553.31227982609, -108128.63674487274},
+ };
+
+ XY s_fwd_expect[] = {
+ {221432.86859285168, 114119.45452653214},
+ {225331.72412711097, -109245.82943505641},
+ {-221432.86859285168, 114119.45452653214},
+ {-225331.72412711097, -109245.82943505641},
+ };
+
+ XY inv_in[] = {
+ { 200, 100},
+ { 200,-100},
+ {-200, 100},
+ {-200,-100}
+ };
+
+ LP e_inv_expect[] = {
+ {0.0017966446840328458, 0.00090435171340223211},
+ {0.0017966164523713021, -0.00090438724081843625},
+ {-0.0017966446840328458, 0.00090435171340223211},
+ {-0.0017966164523713021, -0.00090438724081843625},
+ };
+
+ LP s_inv_expect[] = {
+ {0.0017905070979748127, 0.00089522906964877795},
+ {0.001790479121519977, -0.00089526404022281043},
+ {-0.0017905070979748127, 0.00089522906964877795},
+ {-0.001790479121519977, -0.00089526404022281043},
+ };
+
+ return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
+}
+
+
+#endif
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-02 23:59:33 +0000