diff options
Diffstat (limited to 'src/PJ_lcca.cpp')
| -rw-r--r-- | src/PJ_lcca.cpp | 162 |
1 files changed, 162 insertions, 0 deletions
diff --git a/src/PJ_lcca.cpp b/src/PJ_lcca.cpp new file mode 100644 index 00000000..cbb18709 --- /dev/null +++ b/src/PJ_lcca.cpp @@ -0,0 +1,162 @@ +/***************************************************************************** + + Lambert Conformal Conic Alternative + ----------------------------------- + + This is Gerald Evenden's 2003 implementation of an alternative + "almost" LCC, which has been in use historically, but which + should NOT be used for new projects - i.e: use this implementation + if you need interoperability with old data represented in this + projection, but not in any other case. + + The code was originally discussed on the PROJ.4 mailing list in + a thread archived over at + + http://lists.maptools.org/pipermail/proj/2003-March/000644.html + + It was discussed again in the thread starting at + + http://lists.maptools.org/pipermail/proj/2017-October/007828.html + and continuing at + http://lists.maptools.org/pipermail/proj/2017-November/007831.html + + which prompted Clifford J. Mugnier to add these clarifying notes: + + The French Army Truncated Cubic Lambert (partially conformal) Conic + projection is the Legal system for the projection in France between + the late 1800s and 1948 when the French Legislature changed the law + to recognize the fully conformal version. + + It was (might still be in one or two North African prior French + Colonies) used in North Africa in Algeria, Tunisia, & Morocco, as + well as in Syria during the Levant. + + Last time I have seen it used was about 30+ years ago in + Algeria when it was used to define Lease Block boundaries for + Petroleum Exploration & Production. + + (signed) + + Clifford J. Mugnier, c.p., c.m.s. + Chief of Geodesy + LSU Center for GeoInformatics + Dept. of Civil Engineering + LOUISIANA STATE UNIVERSITY + +*****************************************************************************/ + +#define PJ_LIB__ + +#include <errno.h> +#include <math.h> + +#include "proj.h" +#include "projects.h" + +PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative") + "\n\tConic, Sph&Ell\n\tlat_0="; + +#define MAX_ITER 10 +#define DEL_TOL 1e-12 + +struct pj_opaque { + double *en; + double r0, l, M0; + double C; +}; + + +static double fS(double S, double C) { /* func to compute dr */ + + return S * ( 1. + S * S * C); +} + + +static double fSp(double S, double C) { /* deriv of fs */ + + return 1. + 3.* S * S * C; +} + + +static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ + XY xy = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + double S, r, dr; + + S = pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), Q->en) - Q->M0; + dr = fS(S, Q->C); + r = Q->r0 - dr; + xy.x = P->k0 * (r * sin( lp.lam *= Q->l ) ); + xy.y = P->k0 * (Q->r0 - r * cos(lp.lam) ); + return xy; +} + + +static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */ + LP lp = {0.0,0.0}; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + double theta, dr, S, dif; + int i; + + xy.x /= P->k0; + xy.y /= P->k0; + theta = atan2(xy.x , Q->r0 - xy.y); + dr = xy.y - xy.x * tan(0.5 * theta); + lp.lam = theta / Q->l; + S = dr; + for (i = MAX_ITER; i ; --i) { + S -= (dif = (fS(S, Q->C) - dr) / fSp(S, Q->C)); + if (fabs(dif) < DEL_TOL) break; + } + if (!i) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + lp.phi = pj_inv_mlfn(P->ctx, S + Q->M0, P->es, Q->en); + + return lp; +} + + +static PJ *destructor (PJ *P, int errlev) { + if (0==P) + return 0; + + if (0==P->opaque) + return pj_default_destructor (P, errlev); + + pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en); + return pj_default_destructor (P, errlev); +} + + +PJ *PROJECTION(lcca) { + double s2p0, N0, R0, tan0; + struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + if (0==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + (Q->en = pj_enfn(P->es)); + if (!Q->en) + return pj_default_destructor (P, ENOMEM); + + if (P->phi0 == 0.) { + return destructor(P, PJD_ERR_LAT_0_IS_ZERO); + } + Q->l = sin(P->phi0); + Q->M0 = pj_mlfn(P->phi0, Q->l, cos(P->phi0), Q->en); + s2p0 = Q->l * Q->l; + R0 = 1. / (1. - P->es * s2p0); + N0 = sqrt(R0); + R0 *= P->one_es * N0; + tan0 = tan(P->phi0); + Q->r0 = N0 / tan0; + Q->C = 1. / (6. * R0 * N0); + + P->inv = e_inverse; + P->fwd = e_forward; + P->destructor = destructor; + + return P; +} |