/*****************************************************************************

               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 "proj_internal.h"

PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative")
    "\n\tConic, Sph&Ell\n\tlat_0=";

#define MAX_ITER 10
#define DEL_TOL 1e-12

namespace { // anonymous namespace
struct pj_opaque {
    double *en;
    double r0, l, M0;
    double C;
};
} // anonymous namespace


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 PJ_XY lcca_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
    PJ_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;
    const double lam_mul_l = lp.lam * Q->l;
    xy.x = P->k0 * (r * sin(lam_mul_l));
    xy.y = P->k0 * (Q->r0 - r * cos(lam_mul_l) );
    return xy;
}


static PJ_LP lcca_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
    PJ_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 (nullptr==P)
        return nullptr;

    if (nullptr==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 (nullptr==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 = lcca_e_inverse;
    P->fwd = lcca_e_forward;
    P->destructor = destructor;

    return P;
}