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

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

PROJ_HEAD(gn_sinu, "General Sinusoidal Series") "\n\tPCyl, Sph\n\tm= n=";
PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)") "\n\tPCyl, Sph&Ell";
PROJ_HEAD(eck6, "Eckert VI") "\n\tPCyl, Sph";
PROJ_HEAD(mbtfps, "McBryde-Thomas Flat-Polar Sinusoidal") "\n\tPCyl, Sph";

#define EPS10 1e-10
#define MAX_ITER 8
#define LOOP_TOL 1e-7

namespace { // anonymous namespace
struct pj_gn_sinu_data {
    double *en;
    double m, n, C_x, C_y;
};
} // anonymous namespace

static PJ_XY gn_sinu_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};

    const double s = sin(lp.phi);
    const double c = cos(lp.phi);
    xy.y = pj_mlfn(lp.phi, s, c,
                   static_cast<struct pj_gn_sinu_data *>(P->opaque)->en);
    xy.x = lp.lam * c / sqrt(1. - P->es * s * s);
    return xy;
}

static PJ_LP gn_sinu_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    double s;

    lp.phi =
        pj_inv_mlfn(xy.y, static_cast<struct pj_gn_sinu_data *>(P->opaque)->en);
    s = fabs(lp.phi);
    if (s < M_HALFPI) {
        s = sin(lp.phi);
        lp.lam = xy.x * sqrt(1. - P->es * s * s) / cos(lp.phi);
    } else if ((s - EPS10) < M_HALFPI) {
        lp.lam = 0.;
    } else {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
    }

    return lp;
}

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

    if (Q->m == 0.0)
        lp.phi = Q->n != 1. ? aasin(P->ctx, Q->n * sin(lp.phi)) : lp.phi;
    else {
        int i;

        const double k = Q->n * sin(lp.phi);
        for (i = MAX_ITER; i; --i) {
            const double V =
                (Q->m * lp.phi + sin(lp.phi) - k) / (Q->m + cos(lp.phi));
            lp.phi -= V;
            if (fabs(V) < LOOP_TOL)
                break;
        }
        if (!i) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
    }
    xy.x = Q->C_x * lp.lam * (Q->m + cos(lp.phi));
    xy.y = Q->C_y * lp.phi;

    return xy;
}

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

    xy.y /= Q->C_y;
    lp.phi = (Q->m != 0.0)
                 ? aasin(P->ctx, (Q->m * xy.y + sin(xy.y)) / Q->n)
                 : (Q->n != 1. ? aasin(P->ctx, sin(xy.y) / Q->n) : xy.y);
    lp.lam = xy.x / (Q->C_x * (Q->m + cos(xy.y)));
    return lp;
}

static PJ *pj_gn_sinu_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 pj_gn_sinu_data *>(P->opaque)->en);
    return pj_default_destructor(P, errlev);
}

/* for spheres, only */
static void pj_gn_sinu_setup(PJ *P) {
    struct pj_gn_sinu_data *Q =
        static_cast<struct pj_gn_sinu_data *>(P->opaque);
    P->es = 0;
    P->inv = gn_sinu_s_inverse;
    P->fwd = gn_sinu_s_forward;

    Q->C_y = sqrt((Q->m + 1.) / Q->n);
    Q->C_x = Q->C_y / (Q->m + 1.);
}

PJ *PJ_PROJECTION(sinu) {
    struct pj_gn_sinu_data *Q = static_cast<struct pj_gn_sinu_data *>(
        calloc(1, sizeof(struct pj_gn_sinu_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_gn_sinu_destructor;

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

    if (P->es != 0.0) {
        P->inv = gn_sinu_e_inverse;
        P->fwd = gn_sinu_e_forward;
    } else {
        Q->n = 1.;
        Q->m = 0.;
        pj_gn_sinu_setup(P);
    }
    return P;
}

PJ *PJ_PROJECTION(eck6) {
    struct pj_gn_sinu_data *Q = static_cast<struct pj_gn_sinu_data *>(
        calloc(1, sizeof(struct pj_gn_sinu_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_gn_sinu_destructor;

    Q->m = 1.;
    Q->n = 2.570796326794896619231321691;
    pj_gn_sinu_setup(P);

    return P;
}

PJ *PJ_PROJECTION(mbtfps) {
    struct pj_gn_sinu_data *Q = static_cast<struct pj_gn_sinu_data *>(
        calloc(1, sizeof(struct pj_gn_sinu_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_gn_sinu_destructor;

    Q->m = 0.5;
    Q->n = 1.785398163397448309615660845;
    pj_gn_sinu_setup(P);

    return P;
}

PJ *PJ_PROJECTION(gn_sinu) {
    struct pj_gn_sinu_data *Q = static_cast<struct pj_gn_sinu_data *>(
        calloc(1, sizeof(struct pj_gn_sinu_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_gn_sinu_destructor;

    if (!pj_param(P->ctx, P->params, "tn").i) {
        proj_log_error(P, _("Missing parameter n."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }
    if (!pj_param(P->ctx, P->params, "tm").i) {
        proj_log_error(P, _("Missing parameter m."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    Q->n = pj_param(P->ctx, P->params, "dn").f;
    Q->m = pj_param(P->ctx, P->params, "dm").f;
    if (Q->n <= 0) {
        proj_log_error(P, _("Invalid value for n: it should be > 0."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    if (Q->m < 0) {
        proj_log_error(P, _("Invalid value for m: it should be >= 0."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }

    pj_gn_sinu_setup(P);

    return P;
}