/*
 * Implementation of the Guyou, Pierce Quincuncial, Adams Hemisphere in a
 * Square, Adams World in a Square I & II projections.
 *
 * Based on original code from libproj4 written by Gerald Evenden. Adapted to
 * modern PROJ by Kristian Evers. Original code found in file src/proj_guyou.c,
 * see
 * https://github.com/rouault/libproj4/blob/master/libproject-1.01/src/proj_guyou.c
 * for reference.
 * Fix for Peirce Quincuncial projection to diamond or square by Toby C.
 * Wilkinson to correctly flip out southern hemisphere into the four triangles
 * of Peirce's quincunx. The fix inspired by a similar rotate and translate
 * solution applied by Jonathan Feinberg for cartopy, see
 * https://github.com/jonathf/cartopy/blob/8172cac7fc45cafc86573d408ce85b74258a9c28/lib/cartopy/peircequincuncial.py
 * Added original code for horizontal and vertical arrangement of hemispheres by
 * Toby C. Wilkinson to allow creation of lateral quincuncial projections, such
 * as Grieger's Triptychial, see, e.g.:
 * - Grieger, B. (2020). “Optimized global map projections for specific
 * applications: the triptychial projection and the Spilhaus projection”.
 * EGU2020-9885. https://doi.org/10.5194/egusphere-egu2020-9885
 *
 * Copyright (c) 2005, 2006, 2009 Gerald I. Evenden
 * Copyright (c) 2020 Kristian Evers
 * Copyright (c) 2021 Toby C Wilkinson
 *
 * Related material
 * ----------------
 *
 *  CONFORMAL PROJECTION OF THE SPHERE WITHIN A SQUARE, 1929, OSCAR S. ADAMS,
 *  U.S. COAST AND GEODETIC SURVEY, Special Publication No.153,
 *  ftp://ftp.library.noaa.gov/docs.lib/htdocs/rescue/cgs_specpubs/QB275U35no1531929.pdf
 *
 *  https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection
 *  https://en.wikipedia.org/wiki/Adams_hemisphere-in-a-square_projection
 *  https://en.wikipedia.org/wiki/Peirce_quincuncial_projection
 */

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

#include <algorithm>

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

PROJ_HEAD(guyou, "Guyou") "\n\tMisc Sph No inv";
PROJ_HEAD(peirce_q, "Peirce Quincuncial") "\n\tMisc Sph No inv";
PROJ_HEAD(adams_hemi, "Adams Hemisphere in a Square") "\n\tMisc Sph No inv";
PROJ_HEAD(adams_ws1, "Adams World in a Square I") "\n\tMisc Sph No inv";
PROJ_HEAD(adams_ws2, "Adams World in a Square II") "\n\tMisc Sph No inv";

namespace { // anonymous namespace

enum projection_type {
    GUYOU,
    PEIRCE_Q,
    ADAMS_HEMI,
    ADAMS_WS1,
    ADAMS_WS2,
};

enum peirce_shape {
    PEIRCE_Q_SQUARE,
    PEIRCE_Q_DIAMOND,
    PEIRCE_Q_NHEMISPHERE,
    PEIRCE_Q_SHEMISPHERE,
    PEIRCE_Q_HORIZONTAL,
    PEIRCE_Q_VERTICAL,
};

struct pj_adams_data {
    projection_type mode;
    peirce_shape pqshape;
    double scrollx = 0.0;
    double scrolly = 0.0;
};

} // anonymous namespace

#define TOL 1e-9
#define RSQRT2 0.7071067811865475244008443620

static double ell_int_5(double phi) {
    /* Procedure to compute elliptic integral of the first kind
     * where k^2=0.5.  Precision good to better than 1e-7
     * The approximation is performed with an even Chebyshev
     * series, thus the coefficients below are the even values
     * and where series evaluation  must be multiplied by the argument. */
    constexpr double C0 = 2.19174570831038;
    static const double C[] = {
        -8.58691003636495e-07, 2.02692115653689e-07, 3.12960480765314e-05,
        5.30394739921063e-05,  -0.0012804644680613,  -0.00575574836830288,
        0.0914203033408211,
    };

    double y = phi * M_2_PI;
    y = 2. * y * y - 1.;
    double y2 = 2. * y;
    double d1 = 0.0;
    double d2 = 0.0;
    for (double c : C) {
        double temp = d1;
        d1 = y2 * d1 - d2 + c;
        d2 = temp;
    }

    return phi * (y * d1 - d2 + 0.5 * C0);
}

static PJ_XY adams_forward(PJ_LP lp, PJ *P) {
    double a = 0., b = 0.;
    bool sm = false, sn = false;
    PJ_XY xy;
    const struct pj_adams_data *Q =
        static_cast<const struct pj_adams_data *>(P->opaque);

    switch (Q->mode) {
    case GUYOU:
        if ((fabs(lp.lam) - TOL) > M_PI_2) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return proj_coord_error().xy;
        }

        if (fabs(fabs(lp.phi) - M_PI_2) < TOL) {
            xy.x = 0;
            xy.y = lp.phi < 0 ? -1.85407 : 1.85407;
            return xy;
        } else {
            const double sl = sin(lp.lam);
            const double sp = sin(lp.phi);
            const double cp = cos(lp.phi);
            a = aacos(P->ctx, (cp * sl - sp) * RSQRT2);
            b = aacos(P->ctx, (cp * sl + sp) * RSQRT2);
            sm = lp.lam < 0.;
            sn = lp.phi < 0.;
        }
        break;
    case PEIRCE_Q: {
        /* lam0 - note that the original Peirce model used a central meridian of
         * around -70deg, but the default within proj is +lon0=0 */
        if (Q->pqshape == PEIRCE_Q_NHEMISPHERE) {
            if (lp.phi < -TOL) {
                proj_errno_set(
                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                return proj_coord_error().xy;
            }
        }
        if (Q->pqshape == PEIRCE_Q_SHEMISPHERE) {
            if (lp.phi > -TOL) {
                proj_errno_set(
                    P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                return proj_coord_error().xy;
            }
        }
        const double sl = sin(lp.lam);
        const double cl = cos(lp.lam);
        const double cp = cos(lp.phi);
        a = aacos(P->ctx, cp * (sl + cl) * RSQRT2);
        b = aacos(P->ctx, cp * (sl - cl) * RSQRT2);
        sm = sl < 0.;
        sn = cl > 0.;
    } break;
    case ADAMS_HEMI: {
        const double sp = sin(lp.phi);
        if ((fabs(lp.lam) - TOL) > M_PI_2) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return proj_coord_error().xy;
        }
        a = cos(lp.phi) * sin(lp.lam);
        sm = (sp + a) < 0.;
        sn = (sp - a) < 0.;
        a = aacos(P->ctx, a);
        b = M_PI_2 - lp.phi;
    } break;
    case ADAMS_WS1: {
        const double sp = tan(0.5 * lp.phi);
        b = cos(aasin(P->ctx, sp)) * sin(0.5 * lp.lam);
        a = aacos(P->ctx, (b - sp) * RSQRT2);
        b = aacos(P->ctx, (b + sp) * RSQRT2);
        sm = lp.lam < 0.;
        sn = lp.phi < 0.;
    } break;
    case ADAMS_WS2: {
        const double spp = tan(0.5 * lp.phi);
        a = cos(aasin(P->ctx, spp)) * sin(0.5 * lp.lam);
        sm = (spp + a) < 0.;
        sn = (spp - a) < 0.;
        b = aacos(P->ctx, spp);
        a = aacos(P->ctx, a);
    } break;
    }

    double m = aasin(P->ctx, sqrt((1. + std::min(0.0, cos(a + b)))));
    if (sm)
        m = -m;

    double n = aasin(P->ctx, sqrt(fabs(1. - std::max(0.0, cos(a - b)))));
    if (sn)
        n = -n;

    xy.x = ell_int_5(m);
    xy.y = ell_int_5(n);

    if (Q->mode == PEIRCE_Q) {
        /* Constant complete elliptic integral of the first kind with m=0.5,
         * calculated using
         * https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipk.html
         * . Used as basic as scaled shift distance */
        constexpr double shd = 1.8540746773013719 * 2;

        /* For square and diamond Quincuncial projections, spin out southern
         * hemisphere to triangular segments of quincunx (before rotation for
         * square)*/
        if (Q->pqshape == PEIRCE_Q_SQUARE || (Q->pqshape == PEIRCE_Q_DIAMOND)) {
            if (lp.phi < 0.) { /* fold out segments */
                if (lp.lam < (-0.75 * M_PI))
                    xy.y = shd -
                           xy.y; /* top left segment, shift up and reflect y */
                if ((lp.lam < (-0.25 * M_PI)) && (lp.lam >= (-0.75 * M_PI)))
                    xy.x = -shd -
                           xy.x; /* left segment, shift left and reflect x */
                if ((lp.lam < (0.25 * M_PI)) && (lp.lam >= (-0.25 * M_PI)))
                    xy.y = -shd -
                           xy.y; /* bottom segment, shift down and reflect y */
                if ((lp.lam < (0.75 * M_PI)) && (lp.lam >= (0.25 * M_PI)))
                    xy.x = shd -
                           xy.x; /* right segment, shift right and reflect x */
                if (lp.lam >= (0.75 * M_PI))
                    xy.y = shd -
                           xy.y; /* top right segment, shift up and reflect y */
            }
        }

        /* For square types rotate xy by 45 deg */
        if (Q->pqshape == PEIRCE_Q_SQUARE) {
            const double temp = xy.x;
            xy.x = RSQRT2 * (xy.x - xy.y);
            xy.y = RSQRT2 * (temp + xy.y);
        }

        /* For rectangle Quincuncial projs, spin out southern hemisphere to east
         * (horizontal) or north (vertical) after rotation */
        if (Q->pqshape == PEIRCE_Q_HORIZONTAL) {
            if (lp.phi < 0.) {
                xy.x = shd - xy.x; /* reflect x to east */
            }
            xy.x = xy.x - (shd / 2); /* shift everything so origin is in middle
                                        of two hemispheres */
        }
        if (Q->pqshape == PEIRCE_Q_VERTICAL) {
            if (lp.phi < 0.) {
                xy.y = shd - xy.y; /* reflect y to north */
            }
            xy.y = xy.y - (shd / 2); /* shift everything so origin is in middle
                                        of two hemispheres */
        }

        // if o_scrollx param present, scroll x
        if (!(Q->scrollx == 0.0) && (Q->pqshape == PEIRCE_Q_HORIZONTAL)) {
            double xscale = 2.0;
            double xthresh = shd / 2;
            xy.x = xy.x +
                   (Q->scrollx *
                    (xthresh * 2 * xscale)); /*shift relative to proj width*/
            if (xy.x >= (xthresh * xscale)) {
                xy.x = xy.x - (shd * xscale);
            } else if (xy.x < -(xthresh * xscale)) {
                xy.x = xy.x + (shd * xscale);
            }
        }

        // if o_scrolly param present, scroll y
        if (!(Q->scrolly == 0.0) && (Q->pqshape == PEIRCE_Q_VERTICAL)) {
            double yscale = 2.0;
            double ythresh = shd / 2;
            xy.y = xy.y +
                   (Q->scrolly *
                    (ythresh * 2 * yscale)); /*shift relative to proj height*/
            if (xy.y >= (ythresh * yscale)) {
                xy.y = xy.y - (shd * yscale);
            } else if (xy.y < -(ythresh * yscale)) {
                xy.y = xy.y + (shd * yscale);
            }
        }
    }

    if (Q->mode == ADAMS_HEMI || Q->mode == ADAMS_WS2) { /* rotate by 45deg. */
        const double temp = xy.x;
        xy.x = RSQRT2 * (xy.x - xy.y);
        xy.y = RSQRT2 * (temp + xy.y);
    }

    return xy;
}

static PJ_LP adams_inverse(PJ_XY xy, PJ *P) {
    PJ_LP lp;

    // Only implemented for ADAMS_WS2
    // Uses Newton-Raphson method on the following pair of functions:
    //      f_x(lam,phi) = adams_forward(lam, phi).x - xy.x
    //      f_y(lam,phi) = adams_forward(lam, phi).y - xy.y

    // Initial guess (very rough, especially at high northings)
    // The magic values are got with:
    //  echo 0   90 | src/proj -f "%.8f" +proj=adams_ws2 +R=1
    //  echo 180 0  | src/proj -f "%.8f" +proj=adams_ws2 +R=1
    lp.phi = std::max(std::min(xy.y / 2.62181347, 1.0), -1.0) * M_HALFPI;
    lp.lam =
        fabs(lp.phi) >= M_HALFPI
            ? 0
            : std::max(std::min(xy.x / 2.62205760 / cos(lp.phi), 1.0), -1.0) *
                  M_PI;

    constexpr double deltaXYTolerance = 1e-10;
    return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
}

static PJ_LP peirce_q_square_inverse(PJ_XY xy, PJ *P) {
    /* Heuristics based on trial and repeat process */
    PJ_LP lp;
    lp.phi = 0;
    if (xy.x == 0 && xy.y < 0) {
        lp.lam = -M_PI / 4;
        if (fabs(xy.y) < 2.622057580396)
            lp.phi = M_PI / 4;
    } else if (xy.x > 0 && fabs(xy.y) < 1e-7)
        lp.lam = M_PI / 4;
    else if (xy.x < 0 && fabs(xy.y) < 1e-7) {
        lp.lam = -3 * M_PI / 4;
        lp.phi = M_PI / 2 / 2.622057574224 * xy.x + M_PI / 2;
    } else if (fabs(xy.x) < 1e-7 && xy.y > 0)
        lp.lam = 3 * M_PI / 4;
    else if (xy.x >= 0 && xy.y <= 0) {
        lp.lam = 0;
        if (xy.x == 0 && xy.y == 0) {
            lp.phi = M_PI / 2;
            return lp;
        }
    } else if (xy.x >= 0 && xy.y >= 0)
        lp.lam = M_PI / 2;
    else if (xy.x <= 0 && xy.y >= 0) {
        if (fabs(xy.x) < fabs(xy.y))
            lp.lam = M_PI * 0.9;
        else
            lp.lam = -M_PI * 0.9;
    } else /* if( xy.x <= 0 && xy.y <= 0 ) */
        lp.lam = -M_PI / 2;

    constexpr double deltaXYTolerance = 1e-10;
    return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
}

static PJ_LP peirce_q_diamond_inverse(PJ_XY xy, PJ *P) {
    /* Heuristics based on a trial and repeat process */
    PJ_LP lp;
    lp.phi = 0;
    if (xy.x >= 0 && xy.y <= 0) {
        lp.lam = M_PI / 4;
        if (xy.x > 0 && xy.y == 0) {
            lp.lam = M_PI / 2;
            lp.phi = 0;
        } else if (xy.x == 0 && xy.y == 0) {
            lp.lam = 0;
            lp.phi = M_PI / 2;
            return lp;
        } else if (xy.x == 0 && xy.y < 0) {
            lp.lam = 0;
            lp.phi = M_PI / 4;
        }
    } else if (xy.x >= 0 && xy.y >= 0)
        lp.lam = 3 * M_PI / 4;
    else if (xy.x <= 0 && xy.y >= 0) {
        lp.lam = -3 * M_PI / 4;
    } else /* if( xy.x <= 0 && xy.y <= 0 ) */
        lp.lam = -M_PI / 4;

    if (fabs(xy.x) > 1.8540746773013719 + 1e-3 ||
        fabs(xy.y) > 1.8540746773013719 + 1e-3) {
        lp.phi = -M_PI / 4;
    }

    constexpr double deltaXYTolerance = 1e-10;
    return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
}

static PJ *pj_adams_setup(PJ *P, projection_type mode) {
    struct pj_adams_data *Q = static_cast<struct pj_adams_data *>(
        calloc(1, sizeof(struct pj_adams_data)));

    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    P->es = 0;
    P->fwd = adams_forward;

    Q->mode = mode;
    if (mode == ADAMS_WS2)
        P->inv = adams_inverse;

    if (mode == PEIRCE_Q) {
        // Quincuncial projections shape options: square, diamond, hemisphere,
        // horizontal (rectangle) or vertical (rectangle)
        const char *pqshape = pj_param(P->ctx, P->params, "sshape").s;

        if (!pqshape)
            pqshape = "diamond"; /* default if shape value not supplied */

        if (strcmp(pqshape, "square") == 0) {
            Q->pqshape = PEIRCE_Q_SQUARE;
            P->inv = peirce_q_square_inverse;
        } else if (strcmp(pqshape, "diamond") == 0) {
            Q->pqshape = PEIRCE_Q_DIAMOND;
            P->inv = peirce_q_diamond_inverse;
        } else if (strcmp(pqshape, "nhemisphere") == 0) {
            Q->pqshape = PEIRCE_Q_NHEMISPHERE;
        } else if (strcmp(pqshape, "shemisphere") == 0) {
            Q->pqshape = PEIRCE_Q_SHEMISPHERE;
        } else if (strcmp(pqshape, "horizontal") == 0) {
            Q->pqshape = PEIRCE_Q_HORIZONTAL;
            if (pj_param(P->ctx, P->params, "tscrollx").i) {
                double scrollx;
                scrollx = pj_param(P->ctx, P->params, "dscrollx").f;
                if (scrollx > 1 || scrollx < -1) {
                    proj_log_error(
                        P, _("Invalid value for scrollx: |scrollx| should "
                             "between -1 and 1"));
                    return pj_default_destructor(
                        P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
                }
                Q->scrollx = scrollx;
            }
        } else if (strcmp(pqshape, "vertical") == 0) {
            Q->pqshape = PEIRCE_Q_VERTICAL;
            if (pj_param(P->ctx, P->params, "tscrolly").i) {
                double scrolly;
                scrolly = pj_param(P->ctx, P->params, "dscrolly").f;
                if (scrolly > 1 || scrolly < -1) {
                    proj_log_error(
                        P, _("Invalid value for scrolly: |scrolly| should "
                             "between -1 and 1"));
                    return pj_default_destructor(
                        P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
                }
                Q->scrolly = scrolly;
            }
        } else {
            proj_log_error(P,
                           _("peirce_q: invalid value for 'shape' parameter"));
            return pj_default_destructor(P,
                                         PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
    }

    return P;
}

PJ *PJ_PROJECTION(guyou) { return pj_adams_setup(P, GUYOU); }

PJ *PJ_PROJECTION(peirce_q) { return pj_adams_setup(P, PEIRCE_Q); }

PJ *PJ_PROJECTION(adams_hemi) { return pj_adams_setup(P, ADAMS_HEMI); }

PJ *PJ_PROJECTION(adams_ws1) { return pj_adams_setup(P, ADAMS_WS1); }

PJ *PJ_PROJECTION(adams_ws2) { return pj_adams_setup(P, ADAMS_WS2); }

#undef TOL
#undef RSQRT2