// Changes to handle +over are: Copyright 2011-2014 Morelli Informatik

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

PROJ_HEAD(vandg, "van der Grinten (I)") "\n\tMisc Sph";

#define TOL 1.e-10
#define THIRD .33333333333333333333
#define C2_27 .07407407407407407407   // 2/27
#define PI4_3 4.18879020478639098458  // 4*pi/3
#define PISQ 9.86960440108935861869   // pi^2
#define TPISQ 19.73920880217871723738 // 2*pi^2
#define HPISQ 4.93480220054467930934  // pi^2/2

static PJ_XY vandg_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
    PJ_XY xy = {0.0, 0.0};
    double al, al2, g, g2, p2;
    // Comments tie this formulation to Snyder (1987), p. 241.
    p2 = fabs(lp.phi / M_HALFPI); // sin(theta) from (29-6)
    if ((p2 - TOL) > 1.) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return xy;
    }
    int sign = 1;
    if (P->over && fabs(lp.lam) > M_PI)
        sign = -1;
    if (p2 > 1.)
        p2 = 1.;
    if (fabs(lp.phi) <= TOL) {
        xy.x = lp.lam;
        xy.y = 0.;
    } else if (fabs(lp.lam) <= TOL || fabs(p2 - 1.) < TOL) {
        xy.x = 0.;
        xy.y = M_PI * tan(.5 * asin(p2));
        if (lp.phi < 0.)
            xy.y = -xy.y;
    } else {
        al = .5 * sign * fabs(M_PI / lp.lam - lp.lam / M_PI); // A from (29-3)
        al2 = al * al;                                        // A^2
        g = sqrt(1. - p2 * p2);                               // cos(theta)
        g = g / (p2 + g - 1.);                                // G from (29-4)
        g2 = g * g;                                           // G^2
        p2 = g * (2. / p2 - 1.);                              // P from (29-5)
        p2 = p2 * p2;                                         // P^2
        {
            // Force floating-point computations done on 80 bit on
            // i386 to go back to 64 bit. This is to avoid the test failures
            // of
            // https://github.com/OSGeo/PROJ/issues/1906#issuecomment-583168348
            // The choice of p2 is completely arbitrary.
            volatile double p2_tmp = p2;
            // cppcheck-suppress redundantAssignment
            p2 = p2_tmp;
        }
        xy.x = g - p2; // G - P^2
        g = p2 + al2;  // P^2 + A^2
        // (29-1)
        xy.x = M_PI *
               fabs(al * xy.x + sqrt(al2 * xy.x * xy.x - g * (g2 - p2))) / g;
        if (lp.lam < 0.)
            xy.x = -xy.x;
        xy.y = fabs(xy.x / M_PI);
        // y from (29-2) has been expressed in terms of x here
        xy.y = 1. - xy.y * (xy.y + 2. * al);
        if (xy.y < -TOL) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
        if (xy.y < 0.)
            xy.y = 0.;
        else
            xy.y = sqrt(xy.y) * (lp.phi < 0. ? -M_PI : M_PI);
    }

    return xy;
}

static PJ_LP vandg_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    double t, c0, c1, c2, c3, al, r2, r, m, d, ay, x2, y2;
    // Comments tie this formulation to Snyder (1987), p. 242.
    x2 = xy.x * xy.x; // pi^2 * X^2
    if ((ay = fabs(xy.y)) < TOL) {
        lp.phi = 0.;
        t = x2 * x2 + TPISQ * (x2 + HPISQ);
        lp.lam = fabs(xy.x) <= TOL ? 0. : .5 * (x2 - PISQ + sqrt(t)) / xy.x;
        return (lp);
    }
    y2 = xy.y * xy.y;             // pi^2 * Y^2
    r = x2 + y2;                  // pi^2 * (X^2+Y^2)
    r2 = r * r;                   // pi^4 * (X^2+Y^2)^2
    c1 = -M_PI * ay * (r + PISQ); // pi^4 * c1 (29-11)
    // pi^4 * c3 (29-13)
    c3 = r2 + M_TWOPI * (ay * r + M_PI * (y2 + M_PI * (ay + M_HALFPI)));
    c2 = c1 + PISQ * (r - 3. * y2); // pi^4 * c2 (29-12)
    c0 = M_PI * ay;                 // pi^2 * Y
    c2 /= c3;                       // c2/c3
    al = c1 / c3 - THIRD * c2 * c2; // a1 (29-15)
    m = 2. * sqrt(-THIRD * al);     // m1 (29-16)
    d = C2_27 * c2 * c2 * c2 + (c0 * c0 - THIRD * c2 * c1) / c3; // d (29-14)
    const double al_mul_m = al * m;                              // a1*m1
    if (fabs(al_mul_m) < 1e-16) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return proj_coord_error().lp;
    }
    d = 3. * d / al_mul_m; // cos(3*theta1) (29-17)
    t = fabs(d);
    if ((t - TOL) <= 1.) {
        d = t > 1. ? (d > 0. ? 0. : M_PI) : acos(d); // 3*theta1 (29-17)
        if (r > PISQ) {
            // This code path is triggered for coordinates generated in the
            // forward path when |long|>180deg and +over
            d = M_TWOPI - d;
        }
        // (29-18) but change pi/3 to 4*pi/3 to flip sign of cos
        lp.phi = M_PI * (m * cos(d * THIRD + PI4_3) - THIRD * c2);
        if (xy.y < 0.)
            lp.phi = -lp.phi;
        t = r2 + TPISQ * (x2 - y2 + HPISQ);
        lp.lam = fabs(xy.x) <= TOL
                     ? 0.
                     : .5 * (r - PISQ + (t <= 0. ? 0. : sqrt(t))) / xy.x;
    } else {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return lp;
    }

    return lp;
}

PJ *PJ_PROJECTION(vandg) {
    P->es = 0.;
    P->inv = vandg_s_inverse;
    P->fwd = vandg_s_forward;

    return P;
}

#undef TOL
#undef THIRD
#undef C2_27
#undef PI4_3
#undef PISQ
#undef TPISQ
#undef HPISQ