#define PJ_LIB__

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

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

PROJ_HEAD(imw_p, "International Map of the World Polyconic")
    "\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";

#define TOL 1e-10
#define EPS 1e-10

namespace { // anonymous namespace
enum Mode {
    NONE_IS_ZERO  =  0, /* phi_1 and phi_2 != 0 */
    PHI_1_IS_ZERO =  1, /* phi_1 = 0 */
    PHI_2_IS_ZERO = -1  /* phi_2 = 0 */
};
} // anonymous namespace

namespace { // anonymous namespace
struct pj_opaque {
    double  P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
    double  phi_1, phi_2, lam_1;
    double  *en;
    enum Mode mode;
};
} // anonymous namespace


static int phi12(PJ *P, double *del, double *sig) {
    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
    int err = 0;

    if (!pj_param(P->ctx, P->params, "tlat_1").i ||
        !pj_param(P->ctx, P->params, "tlat_2").i) {
        err = PJD_ERR_LAT_1_2_UNSPECIFIED;
    } else {
        Q->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
        Q->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
        *del = 0.5 * (Q->phi_2 - Q->phi_1);
        *sig = 0.5 * (Q->phi_2 + Q->phi_1);
        err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? PJD_ERR_ABS_LAT1_EQ_ABS_LAT2 : 0;
    }
    return err;
}


static PJ_XY loc_for(PJ_LP lp, PJ *P, double *yc) {
    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
    PJ_XY xy;

    if (lp.phi == 0.0) {
        xy.x = lp.lam;
        xy.y = 0.;
    } else {
        double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;

        sp = sin(lp.phi);
        m = pj_mlfn(lp.phi, sp, cos(lp.phi), Q->en);
        xa = Q->Pp + Q->Qp * m;
        ya = Q->P + Q->Q * m;
        R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
        C = sqrt(R * R - xa * xa);
        if (lp.phi < 0.) C = - C;
        C += ya - R;
        if (Q->mode == PHI_2_IS_ZERO) {
            xb = lp.lam;
            yb = Q->C2;
        } else {
            t = lp.lam * Q->sphi_2;
            xb = Q->R_2 * sin(t);
            yb = Q->C2 + Q->R_2 * (1. - cos(t));
        }
        if (Q->mode == PHI_1_IS_ZERO) {
            xc = lp.lam;
            *yc = 0.;
        } else {
            t = lp.lam * Q->sphi_1;
            xc = Q->R_1 * sin(t);
            *yc = Q->R_1 * (1. - cos(t));
        }
        D = (xb - xc)/(yb - *yc);
        B = xc + D * (C + R - *yc);
        xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
        if (lp.phi > 0)
            xy.x = - xy.x;
        xy.x = (B + xy.x) / (1. + D * D);
        xy.y = sqrt(R * R - xy.x * xy.x);
        if (lp.phi > 0)
            xy.y = - xy.y;
        xy.y += C + R;
    }
    return xy;
}


static PJ_XY imw_p_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
    double yc;
    PJ_XY xy = loc_for(lp, P, &yc);
    return (xy);
}


static PJ_LP imw_p_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);
    PJ_XY t;
    double yc = 0.0;
    int i = 0;
    const int N_MAX_ITER = 1000; /* Arbitrarily chosen number... */

    lp.phi = Q->phi_2;
    lp.lam = xy.x / cos(lp.phi);
    do {
        t = loc_for(lp, P, &yc);
        const double denom = t.y - yc;
        if( denom != 0 || fabs(t.y - xy.y) > TOL )
        {
            if( denom == 0 ) {
                proj_errno_set(P, PJD_ERR_NON_CONVERGENT);
                return proj_coord_error().lp;
            }
            lp.phi = ((lp.phi - Q->phi_1) * (xy.y - yc) / denom) + Q->phi_1;
        }
        if( t.x != 0 && fabs(t.x - xy.x) > TOL )
            lp.lam = lp.lam * xy.x / t.x;
        i ++;
    } while (i < N_MAX_ITER &&
             (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL));

    if( i == N_MAX_ITER )
    {
        proj_errno_set(P, PJD_ERR_NON_CONVERGENT);
        return proj_coord_error().lp;
    }

    return lp;
}


static void xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
    double F;

    *sp = sin(phi);
    *R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
    F = static_cast<struct pj_opaque*>(P->opaque)->lam_1 * *sp;
    *y = *R * (1 - cos(F));
    *x = *R * sin(F);
}


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

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

    if( static_cast<struct pj_opaque*>(P->opaque)->en )
        pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en);

    return pj_default_destructor(P, errlev);
}


PJ *PROJECTION(imw_p) {
    double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
    int err;
    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;

    if (!(Q->en = pj_enfn(P->es))) return pj_default_destructor (P, ENOMEM);
    if( (err = phi12(P, &del, &sig)) != 0) {
        return destructor(P, err);
    }
    if (Q->phi_2 < Q->phi_1) { /* make sure P->phi_1 most southerly */
        del = Q->phi_1;
        Q->phi_1 = Q->phi_2;
        Q->phi_2 = del;
    }
    if (pj_param(P->ctx, P->params, "tlon_1").i)
        Q->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
    else { /* use predefined based upon latitude */
        sig = fabs(sig * RAD_TO_DEG);
        if (sig <= 60)      sig = 2.;
        else if (sig <= 76) sig = 4.;
        else                sig = 8.;
        Q->lam_1 = sig * DEG_TO_RAD;
    }
    Q->mode = NONE_IS_ZERO;
    if (Q->phi_1 != 0.0)
        xy(P, Q->phi_1, &x1, &y1, &Q->sphi_1, &Q->R_1);
    else {
        Q->mode = PHI_1_IS_ZERO;
        y1 = 0.;
        x1 = Q->lam_1;
    }
    if (Q->phi_2 != 0.0)
        xy(P, Q->phi_2, &x2, &T2, &Q->sphi_2, &Q->R_2);
    else {
        Q->mode = PHI_2_IS_ZERO;
        T2 = 0.;
        x2 = Q->lam_1;
    }
    m1 = pj_mlfn(Q->phi_1, Q->sphi_1, cos(Q->phi_1), Q->en);
    m2 = pj_mlfn(Q->phi_2, Q->sphi_2, cos(Q->phi_2), Q->en);
    t = m2 - m1;
    s = x2 - x1;
    y2 = sqrt(t * t - s * s) + y1;
    Q->C2 = y2 - T2;
    t = 1. / t;
    Q->P = (m2 * y1 - m1 * y2) * t;
    Q->Q = (y2 - y1) * t;
    Q->Pp = (m2 * x1 - m1 * x2) * t;
    Q->Qp = (x2 - x1) * t;

    P->fwd = imw_p_e_forward;
    P->inv = imw_p_e_inverse;
    P->destructor = destructor;

    return P;
}