#define PJ_LIB__
#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, PJD_ERR_TOLERANCE_CONDITION);
        return xy;
    }
    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 * 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;
            p2 = p2_tmp;
        }
        xy.x = g - p2;                                 // G - P^2
        g = p2 + al2;                                  // P^2 + A^2
        // (29-1)
        xy.x = M_PI * (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, PJD_ERR_TOLERANCE_CONDITION);
            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, PJD_ERR_TOLERANCE_CONDITION);
        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)
        // (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, PJD_ERR_TOLERANCE_CONDITION);
        return lp;
    }

    return lp;
}


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

    return P;
}

