/* determine small t */

#include "proj.h"

#include "proj_internal.h"

#include <math.h>

double pj_tsfn(double phi, double sinphi, double e) {

    /****************************************************************************

     * Determine function ts(phi) defined in Snyder (1987), Eq. (7-10)

     * Inputs:

     *   phi = geographic latitude (radians)

     *   e = eccentricity of the ellipsoid (dimensionless)

     * Output:

     *   ts = exp(-psi) where psi is the isometric latitude (dimensionless)

     *      = 1 / (tan(chi) + sec(chi))

     * 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

     ***************************************************************************/

    double cosphi = cos(phi);

    // exp(-asinh(tan(phi))) = 1 / (tan(phi) + sec(phi))

    //                       = cos(phi) / (1 + sin(phi)) good for phi > 0

    //                       = (1 - sin(phi)) / cos(phi) good for phi < 0

    return exp(e * atanh(e * sinphi)) *

           (sinphi > 0 ? cosphi / (1 + sinphi) : (1 - sinphi) / cosphi);

}