/* determine small t */
#include <math.h>
#include "proj.h"
#include "proj_internal.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 );
}