#include <math.h>

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

PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph";

// C_x = 2 / sqrt(4 * M_PI + M_PI * M_PI)
constexpr double C_x = .42223820031577120149;

// C_y = 2 * sqrt(M_PI / (4 + M_PI))
constexpr double C_y = 1.32650042817700232218;

// RC_y = 1. / C_y
constexpr double RC_y = .75386330736002178205;

// C_p = 2 + M_PI / 2
constexpr double C_p = 3.57079632679489661922;

// RC_p = 1. / C_p
constexpr double RC_p = .28004957675577868795;

static PJ_XY eck4_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
    PJ_XY xy = {0.0, 0.0};
    int i;
    (void)P;

    // Newton's iterative method to find theta such as
    // theta + sin(theta) * cos(theta) + 2 * sin(theta) == C_p * sin(phi)
    const double p = C_p * sin(lp.phi);
    double V = lp.phi * lp.phi;
    double theta = lp.phi * (0.895168 + V * (0.0218849 + V * 0.00826809));
    constexpr int NITER = 6;
    for (i = NITER; i; --i) {
        const double c = cos(theta);
        const double s = sin(theta);
        V = (theta + s * (c + 2.) - p) / (1. + c * (c + 2.) - s * s);
        theta -= V;
        constexpr double EPS = 1e-7;
        if (fabs(V) < EPS)
            break;
    }
    if (!i) {
        xy.x = C_x * lp.lam;
        xy.y = theta < 0. ? -C_y : C_y;
    } else {
        xy.x = C_x * lp.lam * (1. + cos(theta));
        xy.y = C_y * sin(theta);
    }
    return xy;
}

static PJ_LP eck4_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    PJ_LP lp = {0.0, 0.0};

    const double sin_theta = xy.y * RC_y;
    const double one_minus_abs_sin_theta = 1.0 - fabs(sin_theta);
    if (one_minus_abs_sin_theta >= 0.0 && one_minus_abs_sin_theta <= 1e-12) {
        lp.lam = xy.x / C_x;
        lp.phi = sin_theta > 0 ? M_PI / 2 : -M_PI / 2;
    } else {
        const double theta = aasin(P->ctx, sin_theta);
        const double cos_theta = cos(theta);
        lp.lam = xy.x / (C_x * (1. + cos_theta));
        const double sin_phi = (theta + sin_theta * (cos_theta + 2.)) * RC_p;
        lp.phi = aasin(P->ctx, sin_phi);
    }
    if (!P->over) {
        const double fabs_lam_minus_pi = fabs(lp.lam) - M_PI;
        if (fabs_lam_minus_pi > 0.0) {
            if (fabs_lam_minus_pi > 1e-10) {
                proj_errno_set(
                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                return lp;
            } else {
                lp.lam = lp.lam > 0 ? M_PI : -M_PI;
            }
        }
    }
    return lp;
}

PJ *PJ_PROJECTION(eck4) {
    P->es = 0.0;
    P->inv = eck4_s_inverse;
    P->fwd = eck4_s_forward;

    return P;
}