#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_imw_p_data {
    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_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(P->opaque);
    int err = 0;

    if (!pj_param(P->ctx, P->params, "tlat_1").i) {
        proj_log_error(P, _("Missing parameter: lat_1 should be specified"));
        err = PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    } else if (!pj_param(P->ctx, P->params, "tlat_2").i) {
        proj_log_error(P, _("Missing parameter: lat_2 should be specified"));
        err = PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE;
    } 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)
                  ? PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE
                  : 0;
        if (err) {
            proj_log_error(
                P, _("Illegal value for lat_1 and lat_2: |lat_1 - lat_2| "
                     "and |lat_1 + lat_2| should be > 0"));
        }
    }
    return err;
}

static PJ_XY loc_for(PJ_LP lp, PJ *P, double *yc) {
    struct pj_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(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_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(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, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                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, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        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_imw_p_data *>(P->opaque)->lam_1 * *sp;
    *y = *R * (1 - cos(F));
    *x = *R * sin(F);
}

static PJ *pj_imw_p_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_imw_p_data *>(P->opaque)->en)
        free(static_cast<struct pj_imw_p_data *>(P->opaque)->en);

    return pj_default_destructor(P, errlev);
}

PJ *PJ_PROJECTION(imw_p) {
    double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
    int err;
    struct pj_imw_p_data *Q = static_cast<struct pj_imw_p_data *>(
        calloc(1, sizeof(struct pj_imw_p_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    if (!(Q->en = pj_enfn(P->n)))
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    if ((err = phi12(P, &del, &sig)) != 0) {
        return pj_imw_p_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 = pj_imw_p_destructor;

    return P;
}

#undef TOL
#undef EPS