#define PJ_LIB__

#include <errno.h>
#include <math.h>

#include "proj.h"
#include "proj_internal.h"

PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";


# define C1 .16666666666666666666
# define C2 .00833333333333333333
# define C3 .04166666666666666666
# define C4 .33333333333333333333
# define C5 .06666666666666666666


namespace { // anonymous namespace
struct cass_data {
    double *en;
    double m0;
    bool hyperbolic;
};
} // anonymous namespace



static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);

    const double sinphi = sin (lp.phi);
    const double cosphi = cos (lp.phi);
    const double M = pj_mlfn (lp.phi, sinphi, cosphi, Q->en);

    const double nu_square = 1./(1. - P->es * sinphi*sinphi);
    const double nu  = sqrt(nu_square);
    const double tanphi = tan(lp.phi);
    const double T = tanphi * tanphi;
    const double A = lp.lam * cosphi;
    const double C = P->es * (cosphi * cosphi) / (1 - P->es);
    const double A2 = A * A;

    xy.x = nu * A * (1. - A2 * T *
        (C1 - (8. - T + 8. * C) * A2 * C2));
    xy.y = M - Q->m0 + nu * tanphi * A2 *
        (.5 + (5. - T + 6. * C) * A2 * C3);
    if( Q->hyperbolic )
    {
        const double rho = nu_square * (1. - P->es) * nu;
        xy.y -= xy.y * xy.y * xy.y / (6 * rho * nu);
    }

    return xy;
}


static PJ_XY cass_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
    PJ_XY xy = {0.0, 0.0};
    xy.x  =  asin (cos (lp.phi) * sin (lp.lam));
    xy.y  =  atan2 (tan (lp.phi), cos (lp.lam)) - P->phi0;
    return xy;
}


static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);

    const double phi1 = pj_inv_mlfn (P->ctx, Q->m0 + xy.y, P->es, Q->en);
    const double tanphi1 = tan (phi1);
    const double T1 = tanphi1*tanphi1;
    const double sinphi1 = sin (phi1);
    const double nu1_square = 1. / (1. - P->es * sinphi1 * sinphi1);
    const double nu1 = sqrt (nu1_square);
    const double rho1 = nu1_square * (1. - P->es) * nu1;
    const double D  = xy.x / nu1;
    const double D2 = D * D;
    lp.phi = phi1 - (nu1 * tanphi1 / rho1) * D2 *
        (.5 - (1. + 3. * T1) * D2 * C3);
    lp.lam = D * (1. + T1 * D2 *
        (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1);

    if( Q->hyperbolic )
    {
        // EPSG guidance note 7-2 suggests a custom approximation for the
        // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
        // generic inversion method
        lp = pj_generic_inverse_2d(xy, P, lp);
    }

    return lp;
}


static PJ_LP cass_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
    PJ_LP lp = {0.0,0.0};
    double dd;
    lp.phi = asin(sin(dd = xy.y + P->phi0) * cos(xy.x));
    lp.lam = atan2(tan(xy.x), cos(dd));
    return lp;
}

static PJ *destructor (PJ *P, int errlev) {                        /* Destructor */
    if (nullptr==P)
        return nullptr;

    if (nullptr==P->opaque)
        return pj_default_destructor (P, errlev);

    free (static_cast<struct cass_data*>(P->opaque)->en);
    return pj_default_destructor (P, errlev);
}



PJ *PROJECTION(cass) {

    /* Spheroidal? */
    if (0==P->es) {
        P->inv = cass_s_inverse;
        P->fwd = cass_s_forward;
        return P;
    }

    /* otherwise it's ellipsoidal */
    auto Q = static_cast<struct cass_data*>(calloc (1, sizeof (struct cass_data)));
    P->opaque = Q;
    if (nullptr==P->opaque)
        return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->destructor = destructor;

    Q->en = pj_enfn (P->es);
    if (nullptr==Q->en)
        return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);

    Q->m0 = pj_mlfn (P->phi0,  sin (P->phi0),  cos (P->phi0), Q->en);
    if (pj_param_exists(P->params, "hyperbolic"))
        Q->hyperbolic = true;
    P->inv = cass_e_inverse;
    P->fwd = cass_e_forward;

    return P;
}