diff options
| author | Even Rouault <even.rouault@spatialys.com> | 2018-12-19 13:00:37 +0100 |
|---|---|---|
| committer | Even Rouault <even.rouault@spatialys.com> | 2018-12-26 10:08:55 +0100 |
| commit | df574ae332d57f556fd56314883b3354cab1d0ff (patch) | |
| tree | 63a68d40d7ed7932d6329d9c7baa340bb6294f7f /src/projections/lagrng.cpp | |
| parent | e6de172371ea203f6393d745641d66c82b5b13e2 (diff) | |
| download | PROJ-df574ae332d57f556fd56314883b3354cab1d0ff.tar.gz PROJ-df574ae332d57f556fd56314883b3354cab1d0ff.zip | |
cpp conversion: remove useless pj_, PJ_ and proj_ filename prefixes
Diffstat (limited to 'src/projections/lagrng.cpp')
| -rw-r--r-- | src/projections/lagrng.cpp | 98 |
1 files changed, 98 insertions, 0 deletions
diff --git a/src/projections/lagrng.cpp b/src/projections/lagrng.cpp new file mode 100644 index 00000000..8c0150aa --- /dev/null +++ b/src/projections/lagrng.cpp @@ -0,0 +1,98 @@ +#define PJ_LIB__ +#include <errno.h> +#include <math.h> + +#include "proj.h" +#include "projects.h" + +PROJ_HEAD(lagrng, "Lagrange") "\n\tMisc Sph\n\tW="; + +#define TOL 1e-10 + +namespace { // anonymous namespace +struct pj_opaque { + double a1; + double a2; + double hrw; + double hw; + double rw; + double w; +}; +} // anonymous namespace + + +static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ + XY xy = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + double v, c; + + if (fabs(fabs(lp.phi) - M_HALFPI) < TOL) { + xy.x = 0; + xy.y = lp.phi < 0 ? -2. : 2.; + } else { + lp.phi = sin(lp.phi); + v = Q->a1 * pow((1. + lp.phi)/(1. - lp.phi), Q->hrw); + lp.lam *= Q->rw; + c = 0.5 * (v + 1./v) + cos(lp.lam); + if (c < TOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + xy.x = 2. * sin(lp.lam) / c; + xy.y = (v - 1./v) / c; + } + return xy; +} + + +static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ + LP lp = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + double c, x2, y2p, y2m; + + if (fabs(fabs(xy.y) - 2.) < TOL) { + lp.phi = xy.y < 0 ? -M_HALFPI : M_HALFPI; + lp.lam = 0; + } else { + x2 = xy.x * xy.x; + y2p = 2. + xy.y; + y2m = 2. - xy.y; + c = y2p * y2m - x2; + if (fabs(c) < TOL) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + lp.phi = 2. * atan(pow((y2p * y2p + x2) / (Q->a2 * (y2m * y2m + x2)), Q->hw)) - M_HALFPI; + lp.lam = Q->w * atan2(4. * xy.x, c); + } + return lp; +} + + +PJ *PROJECTION(lagrng) { + double phi1; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + Q->w = pj_param(P->ctx, P->params, "dW").f; + if (Q->w <= 0) + return pj_default_destructor(P, PJD_ERR_W_OR_M_ZERO_OR_LESS); + Q->hw = 0.5 * Q->w; + Q->rw = 1. / Q->w; + Q->hrw = 0.5 * Q->rw; + phi1 = sin(pj_param(P->ctx, P->params, "rlat_1").f); + if (fabs(fabs(phi1) - 1.) < TOL) + return pj_default_destructor(P, PJD_ERR_LAT_LARGER_THAN_90); + + Q->a1 = pow((1. - phi1)/(1. + phi1), Q->hrw); + Q->a2 = Q->a1 * Q->a1; + + P->es = 0.; + P->inv = s_inverse; + P->fwd = s_forward; + + return P; +} + |
