/* Determine latitude angle phi-2. */

#include <algorithm>
#include <limits>
#include <math.h>

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

double pj_sinhpsi2tanphi(PJ_CONTEXT *ctx, const double taup, const double e) {
    /****************************************************************************
     * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi).  The code is taken
     * from GeographicLib::Math::tauf(taup, e).
     *
     * Here
     *   phi = geographic latitude (radians)
     * psi is the isometric latitude
     *   psi = asinh(tan(phi)) - e * atanh(e * sin(phi))
     *       = asinh(tan(chi))
     * chi is the conformal latitude
     *
     * The representation of latitudes via their tangents, tan(phi) and
     *tan(chi), maintains full *relative* accuracy close to latitude = 0 and +/-
     *pi/2. This is sometimes important, e.g., to compute the scale of the
     *transverse Mercator projection which involves cos(phi)/cos(chi) tan(phi)
     *
     * From Karney (2011), Eq. 7,
     *
     *   tau' = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi)))
     *        = tan(phi) * cosh(e * atanh(e * sin(phi))) -
     *          sec(phi) * sinh(e * atanh(e * sin(phi)))
     *        = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma
     * where
     *   sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) ))
     *
     * For e small,
     *
     *    tau' = (1 - e^2) * tau
     *
     * The relation tau'(tau) can therefore by reliably inverted by Newton's
     * method with
     *
     *    tau = tau' / (1 - e^2)
     *
     * as an initial guess.  Newton's method requires dtau'/dtau.  Noting that
     *
     *   dsigma/dtau = e^2 * sqrt(1 + sigma^2) /
     *                 (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2))
     *   d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2)
     *
     * we have
     *
     *   dtau'/dtau = (1 - e^2) * sqrt(1 + tau'^2) * sqrt(1 + tau^2) /
     *                (1 + (1 - e^2) * tau^2)
     *
     * This works fine unless tau^2 and tau'^2 overflows.  This may be partially
     * cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau).  However, nan
     * will still be generated with tau' = inf, since (inf - inf) will appear in
     * the Newton iteration.
     *
     * If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps),
     * sqrt(1 + tau^2) = |tau| and
     *
     *   tau' = exp(- e * atanh(e)) * tau
     *
     * So
     *
     *   tau = exp(e * atanh(e)) * tau'
     *
     * can be returned unless |tau| >= 2/sqrt(eps); this then avoids overflow
     * problems for large tau' and returns the correct result for tau' = +/-inf
     * and nan.
     *
     * Newton's method usually take 2 iterations to converge to double precision
     * accuracy (for WGS84 flattening).  However only 1 iteration is needed for
     * |chi| < 3.35 deg.  In addition, only 1 iteration is needed for |chi| >
     * 89.18 deg (tau' > 70), if tau = exp(e * atanh(e)) * tau' is used as the
     * starting guess.
     ****************************************************************************/

    constexpr int numit = 5;
    // min iterations = 1, max iterations = 2; mean = 1.954
    static const double rooteps = sqrt(std::numeric_limits<double>::epsilon());
    static const double tol = rooteps / 10; // the criterion for Newton's method
    static const double tmax =
        2 / rooteps; // threshold for large arg limit exact
    const double e2m = 1 - e * e;
    const double stol = tol * std::max(1.0, fabs(taup));
    // The initial guess.  70 corresponds to chi = 89.18 deg (see above)
    double tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m;
    if (!(fabs(tau) < tmax)) // handles +/-inf and nan and e = 1
        return tau;
    // If we need to deal with e > 1, then we could include:
    // if (e2m < 0) return std::numeric_limits<double>::quiet_NaN();
    int i = numit;
    for (; i; --i) {
        double tau1 = sqrt(1 + tau * tau);
        double sig = sinh(e * atanh(e * tau / tau1));
        double taupa = sqrt(1 + sig * sig) * tau - sig * tau1;
        double dtau = ((taup - taupa) * (1 + e2m * (tau * tau)) /
                       (e2m * tau1 * sqrt(1 + taupa * taupa)));
        tau += dtau;
        if (!(fabs(dtau) >= stol)) // backwards test to allow nans to succeed.
            break;
    }
    if (i == 0)
        proj_context_errno_set(ctx, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    return tau;
}

/*****************************************************************************/
double pj_phi2(PJ_CONTEXT *ctx, const double ts0, const double e) {
    /****************************************************************************
     * Determine latitude angle phi-2.
     * Inputs:
     *   ts = exp(-psi) where psi is the isometric latitude (dimensionless)
     *        this variable is defined in Snyder (1987), Eq. (7-10)
     *   e = eccentricity of the ellipsoid (dimensionless)
     * Output:
     *   phi = geographic latitude (radians)
     * Here isometric latitude is defined by
     *   psi = log( tan(pi/4 + phi/2) *
     *              ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
     *       = asinh(tan(phi)) - e * atanh(e * sin(phi))
     *       = asinh(tan(chi))
     *   chi = conformal latitude
     *
     * This routine converts t = exp(-psi) to
     *
     *   tau' = tan(chi) = sinh(psi) = (1/t - t)/2
     *
     * returns atan(sinpsi2tanphi(tau'))
     ***************************************************************************/
    return atan(pj_sinhpsi2tanphi(ctx, (1 / ts0 - ts0) / 2, e));
}