1522 lines
58 KiB
C++
1522 lines
58 KiB
C++
/**
|
|
* \file AuxLatitude.cpp
|
|
* \brief Implementation for the GeographicLib::AuxLatitude and
|
|
* GeographicLib::AuxAngle classes.
|
|
*
|
|
* \note This is just sample code. It is not part of GeographicLib itself.
|
|
*
|
|
* This file is an implementation of the methods described in
|
|
* - C. F. F. Karney,
|
|
* On auxiliary latitudes,
|
|
* Technical Report, SRI International, December 2022.
|
|
* https://arxiv.org/abs/2212.05818
|
|
* .
|
|
* Copyright (c) Charles Karney (2022) <charles@karney.com> and licensed under
|
|
* the MIT/X11 License. For more information, see
|
|
* https://geographiclib.sourceforge.io/
|
|
**********************************************************************/
|
|
|
|
#include "AuxLatitude.hpp"
|
|
|
|
/// \cond SKIP
|
|
|
|
#if !defined(AUXLATITUDE_UNOPT)
|
|
// If 1 use the unoptimized formulas in the paper. Otherwise use the
|
|
// optimizations to improve roundoff for high eccentricity.
|
|
#define AUXLATITUDE_UNOPT 0
|
|
#endif
|
|
|
|
#if !defined(AUXLATITUDE_NATIVE_ELLIPTIC)
|
|
// If 1 use the templated implementation of RD and RF in this file. Otherwise
|
|
// use the non-templated versions in the EllipticFunction class.
|
|
#define AUXLATITUDE_NATIVE_ELLIPTIC 1
|
|
#endif
|
|
|
|
/// \endcond
|
|
|
|
#if !AUXLATITUDE_NATIVE_ELLIPTIC
|
|
// Use GeographicLib::EllipticFunction class
|
|
#include <GeographicLib/EllipticFunction.hpp>
|
|
#endif
|
|
|
|
#if defined(_MSC_VER)
|
|
// Squelch warnings about constant conditional expressions
|
|
# pragma warning (disable: 4127)
|
|
#endif
|
|
|
|
namespace GeographicLib {
|
|
|
|
using namespace std;
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxAngle<T>::NaN() {
|
|
return AuxAngle(std::numeric_limits<real>::quiet_NaN(),
|
|
std::numeric_limits<real>::quiet_NaN());
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxAngle<T>::normalized() const {
|
|
if ( isnan( tan() ) ||
|
|
(fabs(_y) > std::numeric_limits<real>::max()/2 &&
|
|
fabs(_x) > std::numeric_limits<real>::max()/2) )
|
|
// deal with
|
|
// (0,0), (inf,inf), (nan,nan), (nan,x), (y,nan), (toobig,toobig)
|
|
return NaN();
|
|
real r = hypot(_y, _x),
|
|
y = _y/r, x = _x/r;
|
|
// deal with r = inf, then one of y,x
|
|
if (isnan(y)) y = copysign(real(1), _y);
|
|
if (isnan(x)) x = copysign(real(1), _x);
|
|
return AuxAngle(y, x);
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxAngle<T>::copyquadrant(const AuxAngle& p) const {
|
|
return AuxAngle(copysign(y(), p.y()), copysign(x(), p.x()));
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T>& AuxAngle<T>::operator+=(const AuxAngle& p) {
|
|
// Do nothing if p.tan() == 0 to preserve signs of y() and x()
|
|
if (p.tan() != 0) {
|
|
real x = _x * p._x - _y * p._y;
|
|
_y = _y * p._x + _x * p._y;
|
|
_x = x;
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
template<typename T>
|
|
AuxLatitude<T>::AuxLatitude(real f)
|
|
: tol_( sqrt(numeric_limits<real>::epsilon()) )
|
|
#if AUXLATITUDE_UNOPT
|
|
, bmin_( numeric_limits<real>::min() )
|
|
, bmax_( numeric_limits<real>::max() )
|
|
#else
|
|
, bmin_( log2(numeric_limits<real>::min()) )
|
|
, bmax_( log2(numeric_limits<real>::max()) )
|
|
#endif
|
|
, _f( f )
|
|
, _fm1( 1 - _f )
|
|
, _e2( _f * (2 - _f) )
|
|
, _e2m1( _fm1 * _fm1 )
|
|
, _e12( _e2/(1 - _e2) )
|
|
, _e12p1( 1 / _e2m1 )
|
|
, _n( _f/(2 - _f) )
|
|
, _e( sqrt(fabs(_e2)) )
|
|
, _e1( sqrt(fabs(_e12)) )
|
|
, _n2( _n * _n )
|
|
#if AUXLATITUDE_UNOPT
|
|
, _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? asinh(_e1) : atan(_e)) / _e) )
|
|
#else
|
|
, _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? atanh(_e) : atan(_e)) / _e) )
|
|
#endif
|
|
{
|
|
if (!(isfinite(_f) && _f < 1 &&
|
|
isfinite(_n) && fabs(_n) < 1))
|
|
throw GeographicErr("Bad ellipsoid parameters");
|
|
fill(_c, _c + Lmax * AUXNUMBER * AUXNUMBER,
|
|
numeric_limits<real>::quiet_NaN());
|
|
}
|
|
|
|
template<typename T>
|
|
AuxLatitude<T>::AuxLatitude(real a, real b)
|
|
: tol_( sqrt(numeric_limits<real>::epsilon()) )
|
|
#if AUXLATITUDE_UNOPT
|
|
, bmin_( numeric_limits<real>::min() )
|
|
, bmax_( numeric_limits<real>::max() )
|
|
#else
|
|
, bmin_( log2(numeric_limits<real>::min()) )
|
|
, bmax_( log2(numeric_limits<real>::max()) )
|
|
#endif
|
|
, _f( (a - b) / a )
|
|
, _fm1( b / a )
|
|
, _e2( ((a - b) * (a + b)) / (a * a) )
|
|
, _e2m1( (b * b) / (a * a) )
|
|
, _e12( ((a - b) * (a + b)) / (b * b) )
|
|
, _e12p1( (a * a) / (b * b) )
|
|
, _n( (a - b) / (a + b) )
|
|
, _e( sqrt(fabs(a - b) * (a + b)) / a )
|
|
, _e1( sqrt(fabs(a - b) * (a + b)) / b )
|
|
, _n2( _n * _n )
|
|
#if AUXLATITUDE_UNOPT
|
|
, _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? asinh(_e1) : atan(_e)) / _e) )
|
|
#else
|
|
, _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? atanh(_e) : atan(_e)) / _e) )
|
|
#endif
|
|
{
|
|
if (!(isfinite(_f) && _f < 1 &&
|
|
isfinite(_n) && fabs(_n) < 1))
|
|
throw GeographicErr("Bad ellipsoid parameters");
|
|
fill(_c, _c + Lmax * AUXNUMBER * AUXNUMBER,
|
|
numeric_limits<real>::quiet_NaN());
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::Parametric(const angle& phi, real* diff) const {
|
|
if (diff) *diff = _fm1;
|
|
return angle(phi.y() * _fm1, phi.x());
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::Geocentric(const angle& phi, real* diff) const {
|
|
if (diff) *diff = _e2m1;
|
|
return angle(phi.y() * _e2m1, phi.x());
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::Rectifying(const angle& phi, real* diff) const {
|
|
angle beta(Parametric(phi).normalized());
|
|
real sbeta = fabs(beta.y()), cbeta = fabs(beta.x());
|
|
real a = 1, b = _fm1, ka = _e2, kb = -_e12, ka1 = _e2m1, kb1 = _e12p1,
|
|
smu, cmu, mr;
|
|
if (_f < 0) {
|
|
swap(a, b); swap(ka, kb); swap(ka1, kb1); swap(sbeta, cbeta);
|
|
}
|
|
// now a,b = larger/smaller semiaxis
|
|
// beta now measured from larger semiaxis
|
|
// kb,ka = modulus-squared for distance from beta = 0,pi/2
|
|
// NB kb <= 0; 0 <= ka <= 1
|
|
// sa = b*E(beta,sqrt(kb)), sb = a*E(beta',sqrt(ka))
|
|
// 1 - ka * (1 - sb2) = 1 -ka + ka*sb2
|
|
real
|
|
sb2 = sbeta * sbeta,
|
|
cb2 = cbeta * cbeta,
|
|
db2 = 1 - kb * sb2,
|
|
da2 = ka1 + ka * sb2,
|
|
// DLMF Eq. 19.25.9
|
|
sa = b * sbeta * ( RF(cb2, db2, 1) - kb * sb2 * RD(cb2, db2, 1) / 3 ),
|
|
// DLMF Eq. 19.25.10 with complementary angles
|
|
sb = a * cbeta * ( ka1 * RF(sb2, da2, 1)
|
|
+ ka * ka1 * cb2 * RD(sb2, 1, da2) / 3
|
|
+ ka * sbeta / sqrt(da2) );
|
|
// sa + sb = 2*EllipticFunction::RG(a*a, b*b) = a*E(e) = b*E(i*e')
|
|
// mr = a*E(e)*(2/pi) = b*E(i*e')*(2/pi)
|
|
mr = (2 * (sa + sb)) / Math::pi<real>();
|
|
smu = sin(sa / mr);
|
|
cmu = sin(sb / mr);
|
|
if (_f < 0) { swap(smu, cmu); swap(a, b); }
|
|
// mu is normalized
|
|
angle mu(angle(smu, cmu).copyquadrant(phi));
|
|
if (diff) {
|
|
real cphi = phi.normalized().x(), tphi = phi.tan();
|
|
if (isfinite(tphi) || isnan(tphi)) {
|
|
cmu = mu.x(); cbeta = beta.x();
|
|
*diff = _fm1 * b/mr * Math::sq(cbeta / cmu) * (cbeta / cphi);
|
|
} else
|
|
*diff = _fm1 * mr/a;
|
|
}
|
|
return mu;
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::Conformal(const angle& phi, real* diff) const {
|
|
real tphi = fabs(phi.tan()), tchi = tphi;
|
|
if ( !( !isfinite(tphi) || tphi == 0 || _f == 0 ) ) {
|
|
real scphi = sc(tphi),
|
|
sig = sinh(_e2 * atanhee(tphi) ),
|
|
scsig = sc(sig);
|
|
#if AUXLATITUDE_UNOPT
|
|
tchi = tphi * scsig - sig * scphi;
|
|
#else
|
|
if (_f <= 0) {
|
|
tchi = tphi * scsig - sig * scphi;
|
|
} else {
|
|
// The general expression for tchi is
|
|
// tphi * scsig - sig * scphi
|
|
// This involves cancellation if f > 0, so change to
|
|
// (tphi - sig) * (tphi + sig) / (tphi * scsig + sig * scphi)
|
|
// To control overflow, write as (sigtphi = sig / tphi)
|
|
// (tphi - sig) * (1 + sigtphi) / (scsig + sigtphi * scphi)
|
|
real sigtphi = sig / tphi, tphimsig;
|
|
if (sig < tphi / 2)
|
|
tphimsig = tphi - sig;
|
|
else {
|
|
// Still have possibly dangerous cancellation in tphi - sig.
|
|
//
|
|
// Write tphi - sig = (1 - e) * Dg(1, e)
|
|
// Dg(x, y) = (g(x) - g(y)) / (x - y)
|
|
// g(x) = sinh(x * atanh(sphi * x))
|
|
// Note sinh(atanh(sphi)) = tphi
|
|
// Turn the crank on divided differences, substitute
|
|
// sphi = tphi/sc(tphi)
|
|
// atanh(x) = asinh(x/sqrt(1-x^2))
|
|
real
|
|
em1 = _e2m1 / (1 + _e), // 1 - e
|
|
atanhs = asinh(tphi), // atanh(sphi)
|
|
scbeta = sc(_fm1 * tphi), // sec(beta)
|
|
scphibeta = sc(tphi) / scbeta, // sec(phi)/sec(beta)
|
|
atanhes = asinh(_e * tphi / scbeta), // atanh(e * sphi)
|
|
t1 = (atanhs - _e * atanhes)/2,
|
|
t2 = asinh(em1 * (tphi * scphibeta)) / em1,
|
|
Dg = cosh((atanhs + _e * atanhes)/2) * (sinh(t1) / t1)
|
|
* ((atanhs + atanhes)/2 + (1 + _e)/2 * t2);
|
|
tphimsig = em1 * Dg; // tphi - sig
|
|
}
|
|
tchi = tphimsig * (1 + sigtphi) / (scsig + sigtphi * scphi);
|
|
}
|
|
#endif
|
|
}
|
|
angle chi(angle(tchi).copyquadrant(phi));
|
|
if (diff) {
|
|
if (isfinite(tphi) || isnan(tphi)) {
|
|
real cchi = chi.normalized().x(),
|
|
cphi = phi.normalized().x(),
|
|
cbeta = Parametric(phi).normalized().x();
|
|
*diff = _e2m1 * (cbeta / cchi) * (cbeta / cphi);
|
|
} else {
|
|
#if AUXLATITUDE_UNOPT
|
|
real ss = _f > 0 ? sinh(_e * asinh(_e1)) : sinh(-_e * atan(_e));
|
|
*diff = _f > 0 ? 1/( sc(ss) + ss ) : sc(ss) - ss;
|
|
#else
|
|
real ss = _f > 0 ? sinh(_e * atanh(_e)) : sinh(-_e * atan(_e));
|
|
*diff = sc(ss) - ss;
|
|
#endif
|
|
}
|
|
}
|
|
return chi;
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::Authalic(const angle& phi, real* diff) const {
|
|
real tphi = fabs(phi.tan());
|
|
angle xi(phi), phin(phi.normalized());
|
|
if ( !( !isfinite(tphi) || tphi == 0 || _f == 0 ) ) {
|
|
real qv = q(tphi),
|
|
Dqp = Dq(tphi),
|
|
Dqm = (_q + qv) / (1 + fabs(phin.y())); // Dq(-tphi)
|
|
xi = angle( copysign(qv, phi.y()), phin.x() * sqrt(Dqp * Dqm) );
|
|
}
|
|
if (diff) {
|
|
if (isfinite(tphi) || isnan(tphi)) {
|
|
real
|
|
cbeta = Parametric(phi).normalized().x(),
|
|
cxi = xi.normalized().x();
|
|
*diff =
|
|
(2/_q) * Math::sq(cbeta / cxi) * (cbeta / cxi) * (cbeta / phin.x());
|
|
} else
|
|
*diff = _e2m1 * sqrt(_q/2);
|
|
}
|
|
return xi;
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::ToAuxiliary(int auxout, const angle& phi,
|
|
real* diff) const {
|
|
switch (auxout) {
|
|
case GEOGRAPHIC: if (diff) *diff = 1; return phi; break;
|
|
case PARAMETRIC: return Parametric(phi, diff); break;
|
|
case GEOCENTRIC: return Geocentric(phi, diff); break;
|
|
case RECTIFYING: return Rectifying(phi, diff); break;
|
|
case CONFORMAL : return Conformal (phi, diff); break;
|
|
case AUTHALIC : return Authalic (phi, diff); break;
|
|
default:
|
|
if (diff) *diff = numeric_limits<real>::quiet_NaN();
|
|
return angle::NaN();
|
|
break;
|
|
}
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::FromAuxiliary(int auxin, const angle& zeta,
|
|
int* niter) const {
|
|
int n = 0; if (niter) *niter = n;
|
|
real tphi = _fm1;
|
|
switch (auxin) {
|
|
case GEOGRAPHIC: return zeta; break;
|
|
// case PARAMETRIC: break;
|
|
case PARAMETRIC: return angle(zeta.y() / _fm1, zeta.x()); break;
|
|
// case GEOCENTRIC: tphi *= _fm1 ; break;
|
|
case GEOCENTRIC: return angle(zeta.y() / _e2m1, zeta.x()); break;
|
|
case RECTIFYING: tphi *= sqrt(_fm1); break;
|
|
case CONFORMAL : tphi *= _fm1 ; break;
|
|
case AUTHALIC : tphi *= cbrt(_fm1); break;
|
|
default: return angle::NaN(); break;
|
|
}
|
|
|
|
// Drop through to solution by Newton's method
|
|
#if AUXLATITUDE_UNOPT
|
|
real tzeta = fabs(zeta.tan());
|
|
if (!isfinite(tzeta)) return zeta;
|
|
tphi = tzeta / tphi;
|
|
real bmin = fmin(tphi, bmin_), bmax = fmax(tphi, bmax_);
|
|
for (int sign = 0, osign = 0, ntrip = 0; n < numit_;) {
|
|
++n;
|
|
real diff;
|
|
angle zeta1(ToAuxiliary(auxin, angle(tphi), &diff));
|
|
real tzeta1 = zeta1.tan();
|
|
osign = sign;
|
|
if (tzeta1 == tzeta)
|
|
break;
|
|
else if (tzeta1 > tzeta) {
|
|
sign = 1;
|
|
bmax = tphi;
|
|
} else {
|
|
sign = -1;
|
|
bmin = tphi;
|
|
}
|
|
real dtphi = -(tzeta1 - tzeta) / diff;
|
|
tphi += dtphi;
|
|
if (!(fabs(dtphi) >= tol_ * tphi))
|
|
break;
|
|
// Skip the Newton guards if UNOPT
|
|
if (false && ((sign * osign < 0 && n - ntrip > 2) ||
|
|
tphi >= bmax || tphi <= bmin)) {
|
|
sign = 0; ntrip = n;
|
|
tphi = sqrt(bmin * bmax);
|
|
}
|
|
}
|
|
#else
|
|
real tzeta = fabs(zeta.tan()), ltzeta = log2(tzeta);
|
|
if (!isfinite(ltzeta)) return zeta;
|
|
tphi = tzeta / tphi;
|
|
real ltphi = log2(tphi),
|
|
bmin = fmin(ltphi, bmin_), bmax = fmax(ltphi, bmax_);
|
|
for (int sign = 0, osign = 0, ntrip = 0; n < numit_;) {
|
|
++n;
|
|
real diff;
|
|
angle zeta1(ToAuxiliary(auxin, angle(tphi), &diff));
|
|
real tzeta1 = zeta1.tan(), ltzeta1 = log2(tzeta1);
|
|
// Convert derivative from dtan(zeta)/dtan(phi) to
|
|
// dlog(tan(zeta))/dlog(tan(phi))
|
|
diff *= tphi/tzeta1;
|
|
osign = sign;
|
|
if (tzeta1 == tzeta)
|
|
break;
|
|
else if (tzeta1 > tzeta) {
|
|
sign = 1;
|
|
bmax = ltphi;
|
|
} else {
|
|
sign = -1;
|
|
bmin = ltphi;
|
|
}
|
|
real dltphi = -(ltzeta1 - ltzeta) / diff;
|
|
ltphi += dltphi;
|
|
tphi = exp2(ltphi);
|
|
if (!(fabs(dltphi) >= tol_)) {
|
|
++n;
|
|
// Final Newton iteration without the logs
|
|
zeta1 = ToAuxiliary(auxin, angle(tphi), &diff);
|
|
tphi -= (zeta1.tan() - tzeta) / diff;
|
|
break;
|
|
}
|
|
if ((sign * osign < 0 && n - ntrip > 2) ||
|
|
ltphi >= bmax || ltphi <= bmin) {
|
|
sign = 0; ntrip = n;
|
|
ltphi = (bmin + bmax) / 2;
|
|
tphi = exp2(ltphi);
|
|
}
|
|
}
|
|
#endif
|
|
if (niter) *niter = n;
|
|
return angle(tphi).copyquadrant(zeta);
|
|
}
|
|
|
|
template<typename T>
|
|
AuxAngle<T> AuxLatitude<T>::Convert(int auxin, int auxout, const angle& zeta,
|
|
bool series) const {
|
|
int k = ind(auxout, auxin);
|
|
if (k < 0) return angle::NaN();
|
|
if (auxin == auxout) return zeta;
|
|
if (series) {
|
|
if ( isnan(_c[Lmax * k]) ) fillcoeff(auxin, auxout, k);
|
|
angle zetan(zeta.normalized());
|
|
real d = Clenshaw(zetan.y(), zetan.x(), _c + Lmax * k, Lmax);
|
|
zetan += angle::radians(d);
|
|
return zetan;
|
|
} else {
|
|
if (auxin < 3 && auxout < 3)
|
|
return angle(zeta.y() * pow(_fm1, auxout - auxin), zeta.x());
|
|
else
|
|
return ToAuxiliary(auxout, FromAuxiliary(auxin, zeta));
|
|
}
|
|
}
|
|
|
|
template<typename T>
|
|
typename AuxLatitude<T>::real AuxLatitude<T>::RD(real x, real y, real z) {
|
|
#if AUXLATITUDE_NATIVE_ELLIPTIC
|
|
// Carlson, eqs 2.28 - 2.34
|
|
static const real
|
|
tolRD = pow(real(0.2) * (numeric_limits<real>::epsilon() * real(0.01)),
|
|
1/real(8));
|
|
real
|
|
A0 = (x + y + 3*z)/5,
|
|
An = A0,
|
|
Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRD,
|
|
x0 = x,
|
|
y0 = y,
|
|
z0 = z,
|
|
mul = 1,
|
|
s = 0;
|
|
while (Q >= mul * fabs(An)) {
|
|
// Max 7 trips
|
|
real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
|
|
s += 1/(mul * sqrt(z0) * (z0 + lam));
|
|
An = (An + lam)/4;
|
|
x0 = (x0 + lam)/4;
|
|
y0 = (y0 + lam)/4;
|
|
z0 = (z0 + lam)/4;
|
|
mul *= 4;
|
|
}
|
|
real
|
|
X = (A0 - x) / (mul * An),
|
|
Y = (A0 - y) / (mul * An),
|
|
Z = -(X + Y) / 3,
|
|
E2 = X*Y - 6*Z*Z,
|
|
E3 = (3*X*Y - 8*Z*Z)*Z,
|
|
E4 = 3 * (X*Y - Z*Z) * Z*Z,
|
|
E5 = X*Y*Z*Z*Z;
|
|
// https://dlmf.nist.gov/19.36.E2
|
|
// Polynomial is
|
|
// (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
|
|
// - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
|
|
// - 9*(E3*E4+E2*E5)/68)
|
|
return ((471240 - 540540 * E2) * E5 +
|
|
(612612 * E2 - 540540 * E3 - 556920) * E4 +
|
|
E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
|
|
E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
|
|
(4084080 * mul * An * sqrt(An)) + 3 * s;
|
|
#else
|
|
return real(EllipticFunction::RD(Math::real(x), Math::real(y),
|
|
Math::real(z)));
|
|
#endif
|
|
}
|
|
|
|
template<typename T>
|
|
typename AuxLatitude<T>::real AuxLatitude<T>::RF(real x, real y, real z) {
|
|
#if AUXLATITUDE_NATIVE_ELLIPTIC
|
|
// Carlson, eqs 2.2 - 2.7
|
|
static const real tolRF =
|
|
pow(3 * numeric_limits<real>::epsilon() * real(0.01), 1/real(8));
|
|
real
|
|
A0 = (x + y + z)/3,
|
|
An = A0,
|
|
Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRF,
|
|
x0 = x,
|
|
y0 = y,
|
|
z0 = z,
|
|
mul = 1;
|
|
while (Q >= mul * fabs(An)) {
|
|
// Max 6 trips
|
|
real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
|
|
An = (An + lam)/4;
|
|
x0 = (x0 + lam)/4;
|
|
y0 = (y0 + lam)/4;
|
|
z0 = (z0 + lam)/4;
|
|
mul *= 4;
|
|
}
|
|
real
|
|
X = (A0 - x) / (mul * An),
|
|
Y = (A0 - y) / (mul * An),
|
|
Z = - (X + Y),
|
|
E2 = X*Y - Z*Z,
|
|
E3 = X*Y*Z;
|
|
// https://dlmf.nist.gov/19.36.E1
|
|
// Polynomial is
|
|
// (1 - E2/10 + E3/14 + E2^2/24 - 3*E2*E3/44
|
|
// - 5*E2^3/208 + 3*E3^2/104 + E2^2*E3/16)
|
|
// convert to Horner form...
|
|
return (E3 * (6930 * E3 + E2 * (15015 * E2 - 16380) + 17160) +
|
|
E2 * ((10010 - 5775 * E2) * E2 - 24024) + 240240) /
|
|
(240240 * sqrt(An));
|
|
#else
|
|
return real(EllipticFunction::RF(Math::real(x), Math::real(y),
|
|
Math::real(z)));
|
|
#endif
|
|
}
|
|
|
|
template<typename T>
|
|
typename AuxLatitude<T>::real AuxLatitude<T>::atanhee(real tphi) const {
|
|
#if AUXLATITUDE_UNOPT
|
|
real sphi = sn(tphi);
|
|
return _f == 0 ? sphi :
|
|
(_f < 0 ? atan( _e * sphi ) : atanh( _e * sphi )) / _e;
|
|
#else
|
|
real s = _f <= 0 ? sn(tphi) : sn(_fm1 * tphi);
|
|
return _f == 0 ? s :
|
|
// atanh(e * sphi) = asinh(e' * sbeta)
|
|
(_f < 0 ? atan( _e * s ) : asinh( _e1 * s )) / _e;
|
|
#endif
|
|
}
|
|
|
|
template<typename T>
|
|
typename AuxLatitude<T>::real AuxLatitude<T>::q(real tphi) const {
|
|
#if AUXLATITUDE_UNOPT
|
|
real sphi = sn(tphi);
|
|
return atanhee(tphi) + sphi / (1 - _e2 * sphi*sphi);
|
|
#else
|
|
real scbeta = sc(_fm1 * tphi);
|
|
return atanhee(tphi) + (tphi / scbeta) * (sc(tphi) / scbeta);
|
|
#endif
|
|
}
|
|
|
|
template<typename T>
|
|
typename AuxLatitude<T>::real AuxLatitude<T>::Dq(real tphi) const {
|
|
#if AUXLATITUDE_UNOPT
|
|
real sphi = sn(tphi), d = 1 - sphi;
|
|
if (tphi <= 0)
|
|
// This branch is not reached; this case is open-coded in Authalic.
|
|
return (_q - q(tphi)) / d;
|
|
else if (d == 0)
|
|
return 2 / Math::sq(_e2m1);
|
|
else
|
|
return (_f == 0 ? 1 / (1 - _e2 * sphi) :
|
|
(_f < 0 ? atan(_e * d / (1 - _e2 * sphi)) / (_e * d) :
|
|
atanh(_e * d / (1 - _e2 * sphi)) / (_e * d) )) +
|
|
(1 + _e2 * sphi) / ((1 - _e2 * sphi*sphi) * _e2m1);
|
|
#else
|
|
real scphi = sc(tphi), sphi = sn(tphi),
|
|
// d = (1 - sphi) can underflow to zero for large tphi
|
|
d = tphi > 0 ? 1 / (scphi * scphi * (1 + sphi)) : 1 - sphi;
|
|
if (tphi <= 0)
|
|
// This branch is not reached; this case is open-coded in Authalic.
|
|
return (_q - q(tphi)) / d;
|
|
else if (d == 0)
|
|
return 2 / Math::sq(_e2m1);
|
|
else {
|
|
// General expression for Dq(1, sphi) is
|
|
// atanh(e * d / (1 - e2 * sphi)) / (e * d) +
|
|
// (1 + e2 * sphi) / ((1 - e2 * sphi * sphi) * e2m1);
|
|
// atanh( e * d / (1 - e2 * sphi))
|
|
// = atanh( e * d * scphi/(scphi - e2 * tphi))
|
|
// =
|
|
real scbeta = sc(_fm1 * tphi);
|
|
return (_f == 0 ? 1 :
|
|
(_f > 0 ? asinh(_e1 * d * scphi / scbeta) :
|
|
atan(_e * d / (1 - _e2 * sphi))) / (_e * d) ) +
|
|
(_f > 0 ?
|
|
((scphi + _e2 * tphi) / (_e2m1 * scbeta)) * (scphi / scbeta) :
|
|
(1 + _e2 * sphi) / ((1 - _e2 * sphi*sphi) * _e2m1) );
|
|
}
|
|
#endif
|
|
}
|
|
|
|
template<typename T>
|
|
void AuxLatitude<T>::fillcoeff(int auxin, int auxout, int k) const {
|
|
#if GEOGRAPHICLIB_AUXLATITUDE_ORDER == 4
|
|
static const real coeffs[] = {
|
|
// C[phi,phi] skipped
|
|
// C[phi,beta]; even coeffs only
|
|
0, 1,
|
|
0, 1/real(2),
|
|
1/real(3),
|
|
1/real(4),
|
|
// C[phi,theta]; even coeffs only
|
|
-2, 2,
|
|
-4, 2,
|
|
8/real(3),
|
|
4,
|
|
// C[phi,mu]; even coeffs only
|
|
-27/real(32), 3/real(2),
|
|
-55/real(32), 21/real(16),
|
|
151/real(96),
|
|
1097/real(512),
|
|
// C[phi,chi]
|
|
116/real(45), -2, -2/real(3), 2,
|
|
-227/real(45), -8/real(5), 7/real(3),
|
|
-136/real(35), 56/real(15),
|
|
4279/real(630),
|
|
// C[phi,xi]
|
|
-2582/real(14175), -16/real(35), 4/real(45), 4/real(3),
|
|
-11966/real(14175), 152/real(945), 46/real(45),
|
|
3802/real(14175), 3044/real(2835),
|
|
6059/real(4725),
|
|
// C[beta,phi]; even coeffs only
|
|
0, -1,
|
|
0, 1/real(2),
|
|
-1/real(3),
|
|
1/real(4),
|
|
// C[beta,beta] skipped
|
|
// C[beta,theta]; even coeffs only
|
|
0, 1,
|
|
0, 1/real(2),
|
|
1/real(3),
|
|
1/real(4),
|
|
// C[beta,mu]; even coeffs only
|
|
-9/real(32), 1/real(2),
|
|
-37/real(96), 5/real(16),
|
|
29/real(96),
|
|
539/real(1536),
|
|
// C[beta,chi]
|
|
38/real(45), -1/real(3), -2/real(3), 1,
|
|
-7/real(9), -14/real(15), 5/real(6),
|
|
-34/real(21), 16/real(15),
|
|
2069/real(1260),
|
|
// C[beta,xi]
|
|
-1082/real(14175), -46/real(315), 4/real(45), 1/real(3),
|
|
-338/real(2025), 68/real(945), 17/real(90),
|
|
1102/real(14175), 461/real(2835),
|
|
3161/real(18900),
|
|
// C[theta,phi]; even coeffs only
|
|
2, -2,
|
|
-4, 2,
|
|
-8/real(3),
|
|
4,
|
|
// C[theta,beta]; even coeffs only
|
|
0, -1,
|
|
0, 1/real(2),
|
|
-1/real(3),
|
|
1/real(4),
|
|
// C[theta,theta] skipped
|
|
// C[theta,mu]; even coeffs only
|
|
-23/real(32), -1/real(2),
|
|
-5/real(96), 5/real(16),
|
|
1/real(32),
|
|
283/real(1536),
|
|
// C[theta,chi]
|
|
4/real(9), -2/real(3), -2/real(3), 0,
|
|
-23/real(45), -4/real(15), 1/real(3),
|
|
-24/real(35), 2/real(5),
|
|
83/real(126),
|
|
// C[theta,xi]
|
|
-2102/real(14175), -158/real(315), 4/real(45), -2/real(3),
|
|
934/real(14175), -16/real(945), 16/real(45),
|
|
922/real(14175), -232/real(2835),
|
|
719/real(4725),
|
|
// C[mu,phi]; even coeffs only
|
|
9/real(16), -3/real(2),
|
|
-15/real(32), 15/real(16),
|
|
-35/real(48),
|
|
315/real(512),
|
|
// C[mu,beta]; even coeffs only
|
|
3/real(16), -1/real(2),
|
|
1/real(32), -1/real(16),
|
|
-1/real(48),
|
|
-5/real(512),
|
|
// C[mu,theta]; even coeffs only
|
|
13/real(16), 1/real(2),
|
|
33/real(32), -1/real(16),
|
|
-5/real(16),
|
|
-261/real(512),
|
|
// C[mu,mu] skipped
|
|
// C[mu,chi]
|
|
41/real(180), 5/real(16), -2/real(3), 1/real(2),
|
|
557/real(1440), -3/real(5), 13/real(48),
|
|
-103/real(140), 61/real(240),
|
|
49561/real(161280),
|
|
// C[mu,xi]
|
|
-1609/real(28350), 121/real(1680), 4/real(45), -1/real(6),
|
|
16463/real(453600), 26/real(945), -29/real(720),
|
|
449/real(28350), -1003/real(45360),
|
|
-40457/real(2419200),
|
|
// C[chi,phi]
|
|
-82/real(45), 4/real(3), 2/real(3), -2,
|
|
-13/real(9), -16/real(15), 5/real(3),
|
|
34/real(21), -26/real(15),
|
|
1237/real(630),
|
|
// C[chi,beta]
|
|
-16/real(45), 0, 2/real(3), -1,
|
|
19/real(45), -2/real(5), 1/real(6),
|
|
16/real(105), -1/real(15),
|
|
17/real(1260),
|
|
// C[chi,theta]
|
|
-2/real(9), 2/real(3), 2/real(3), 0,
|
|
43/real(45), 4/real(15), -1/real(3),
|
|
2/real(105), -2/real(5),
|
|
-55/real(126),
|
|
// C[chi,mu]
|
|
1/real(360), -37/real(96), 2/real(3), -1/real(2),
|
|
437/real(1440), -1/real(15), -1/real(48),
|
|
37/real(840), -17/real(480),
|
|
-4397/real(161280),
|
|
// C[chi,chi] skipped
|
|
// C[chi,xi]
|
|
-2312/real(14175), -88/real(315), 34/real(45), -2/real(3),
|
|
6079/real(14175), -184/real(945), 1/real(45),
|
|
772/real(14175), -106/real(2835),
|
|
-167/real(9450),
|
|
// C[xi,phi]
|
|
538/real(4725), 88/real(315), -4/real(45), -4/real(3),
|
|
-2482/real(14175), 8/real(105), 34/real(45),
|
|
-898/real(14175), -1532/real(2835),
|
|
6007/real(14175),
|
|
// C[xi,beta]
|
|
34/real(675), 32/real(315), -4/real(45), -1/real(3),
|
|
74/real(2025), -4/real(315), -7/real(90),
|
|
2/real(14175), -83/real(2835),
|
|
-797/real(56700),
|
|
// C[xi,theta]
|
|
778/real(4725), 62/real(105), -4/real(45), 2/real(3),
|
|
12338/real(14175), -32/real(315), 4/real(45),
|
|
-1618/real(14175), -524/real(2835),
|
|
-5933/real(14175),
|
|
// C[xi,mu]
|
|
1297/real(18900), -817/real(10080), -4/real(45), 1/real(6),
|
|
-29609/real(453600), -2/real(35), 49/real(720),
|
|
-2917/real(56700), 4463/real(90720),
|
|
331799/real(7257600),
|
|
// C[xi,chi]
|
|
2458/real(4725), 46/real(315), -34/real(45), 2/real(3),
|
|
3413/real(14175), -256/real(315), 19/real(45),
|
|
-15958/real(14175), 248/real(567),
|
|
16049/real(28350),
|
|
// C[xi,xi] skipped
|
|
};
|
|
static const int ptrs[] = {
|
|
0, 0, 6, 12, 18, 28, 38, 44, 44, 50, 56, 66, 76, 82, 88, 88, 94, 104,
|
|
114, 120, 126, 132, 132, 142, 152, 162, 172, 182, 192, 192, 202, 212,
|
|
222, 232, 242, 252, 252,
|
|
};
|
|
#elif GEOGRAPHICLIB_AUXLATITUDE_ORDER == 6
|
|
static const real coeffs[] = {
|
|
// C[phi,phi] skipped
|
|
// C[phi,beta]; even coeffs only
|
|
0, 0, 1,
|
|
0, 0, 1/real(2),
|
|
0, 1/real(3),
|
|
0, 1/real(4),
|
|
1/real(5),
|
|
1/real(6),
|
|
// C[phi,theta]; even coeffs only
|
|
2, -2, 2,
|
|
6, -4, 2,
|
|
-8, 8/real(3),
|
|
-16, 4,
|
|
32/real(5),
|
|
32/real(3),
|
|
// C[phi,mu]; even coeffs only
|
|
269/real(512), -27/real(32), 3/real(2),
|
|
6759/real(4096), -55/real(32), 21/real(16),
|
|
-417/real(128), 151/real(96),
|
|
-15543/real(2560), 1097/real(512),
|
|
8011/real(2560),
|
|
293393/real(61440),
|
|
// C[phi,chi]
|
|
-2854/real(675), 26/real(45), 116/real(45), -2, -2/real(3), 2,
|
|
2323/real(945), 2704/real(315), -227/real(45), -8/real(5), 7/real(3),
|
|
73814/real(2835), -1262/real(105), -136/real(35), 56/real(15),
|
|
-399572/real(14175), -332/real(35), 4279/real(630),
|
|
-144838/real(6237), 4174/real(315),
|
|
601676/real(22275),
|
|
// C[phi,xi]
|
|
28112932/real(212837625), 60136/real(467775), -2582/real(14175),
|
|
-16/real(35), 4/real(45), 4/real(3),
|
|
251310128/real(638512875), -21016/real(51975), -11966/real(14175),
|
|
152/real(945), 46/real(45),
|
|
-8797648/real(10945935), -94388/real(66825), 3802/real(14175),
|
|
3044/real(2835),
|
|
-1472637812/real(638512875), 41072/real(93555), 6059/real(4725),
|
|
455935736/real(638512875), 768272/real(467775),
|
|
4210684958LL/real(1915538625),
|
|
// C[beta,phi]; even coeffs only
|
|
0, 0, -1,
|
|
0, 0, 1/real(2),
|
|
0, -1/real(3),
|
|
0, 1/real(4),
|
|
-1/real(5),
|
|
1/real(6),
|
|
// C[beta,beta] skipped
|
|
// C[beta,theta]; even coeffs only
|
|
0, 0, 1,
|
|
0, 0, 1/real(2),
|
|
0, 1/real(3),
|
|
0, 1/real(4),
|
|
1/real(5),
|
|
1/real(6),
|
|
// C[beta,mu]; even coeffs only
|
|
205/real(1536), -9/real(32), 1/real(2),
|
|
1335/real(4096), -37/real(96), 5/real(16),
|
|
-75/real(128), 29/real(96),
|
|
-2391/real(2560), 539/real(1536),
|
|
3467/real(7680),
|
|
38081/real(61440),
|
|
// C[beta,chi]
|
|
-3118/real(4725), -1/real(3), 38/real(45), -1/real(3), -2/real(3), 1,
|
|
-247/real(270), 50/real(21), -7/real(9), -14/real(15), 5/real(6),
|
|
17564/real(2835), -5/real(3), -34/real(21), 16/real(15),
|
|
-49877/real(14175), -28/real(9), 2069/real(1260),
|
|
-28244/real(4455), 883/real(315),
|
|
797222/real(155925),
|
|
// C[beta,xi]
|
|
7947332/real(212837625), 11824/real(467775), -1082/real(14175),
|
|
-46/real(315), 4/real(45), 1/real(3),
|
|
39946703/real(638512875), -16672/real(155925), -338/real(2025),
|
|
68/real(945), 17/real(90),
|
|
-255454/real(1563705), -101069/real(467775), 1102/real(14175),
|
|
461/real(2835),
|
|
-189032762/real(638512875), 1786/real(18711), 3161/real(18900),
|
|
80274086/real(638512875), 88868/real(467775),
|
|
880980241/real(3831077250LL),
|
|
// C[theta,phi]; even coeffs only
|
|
-2, 2, -2,
|
|
6, -4, 2,
|
|
8, -8/real(3),
|
|
-16, 4,
|
|
-32/real(5),
|
|
32/real(3),
|
|
// C[theta,beta]; even coeffs only
|
|
0, 0, -1,
|
|
0, 0, 1/real(2),
|
|
0, -1/real(3),
|
|
0, 1/real(4),
|
|
-1/real(5),
|
|
1/real(6),
|
|
// C[theta,theta] skipped
|
|
// C[theta,mu]; even coeffs only
|
|
499/real(1536), -23/real(32), -1/real(2),
|
|
6565/real(12288), -5/real(96), 5/real(16),
|
|
-77/real(128), 1/real(32),
|
|
-4037/real(7680), 283/real(1536),
|
|
1301/real(7680),
|
|
17089/real(61440),
|
|
// C[theta,chi]
|
|
-3658/real(4725), 2/real(9), 4/real(9), -2/real(3), -2/real(3), 0,
|
|
61/real(135), 68/real(45), -23/real(45), -4/real(15), 1/real(3),
|
|
9446/real(2835), -46/real(35), -24/real(35), 2/real(5),
|
|
-34712/real(14175), -80/real(63), 83/real(126),
|
|
-2362/real(891), 52/real(45),
|
|
335882/real(155925),
|
|
// C[theta,xi]
|
|
216932/real(2627625), 109042/real(467775), -2102/real(14175),
|
|
-158/real(315), 4/real(45), -2/real(3),
|
|
117952358/real(638512875), -7256/real(155925), 934/real(14175),
|
|
-16/real(945), 16/real(45),
|
|
-7391576/real(54729675), -25286/real(66825), 922/real(14175),
|
|
-232/real(2835),
|
|
-67048172/real(638512875), 268/real(18711), 719/real(4725),
|
|
46774256/real(638512875), 14354/real(467775),
|
|
253129538/real(1915538625),
|
|
// C[mu,phi]; even coeffs only
|
|
-3/real(32), 9/real(16), -3/real(2),
|
|
135/real(2048), -15/real(32), 15/real(16),
|
|
105/real(256), -35/real(48),
|
|
-189/real(512), 315/real(512),
|
|
-693/real(1280),
|
|
1001/real(2048),
|
|
// C[mu,beta]; even coeffs only
|
|
-1/real(32), 3/real(16), -1/real(2),
|
|
-9/real(2048), 1/real(32), -1/real(16),
|
|
3/real(256), -1/real(48),
|
|
3/real(512), -5/real(512),
|
|
-7/real(1280),
|
|
-7/real(2048),
|
|
// C[mu,theta]; even coeffs only
|
|
-15/real(32), 13/real(16), 1/real(2),
|
|
-1673/real(2048), 33/real(32), -1/real(16),
|
|
349/real(256), -5/real(16),
|
|
963/real(512), -261/real(512),
|
|
-921/real(1280),
|
|
-6037/real(6144),
|
|
// C[mu,mu] skipped
|
|
// C[mu,chi]
|
|
7891/real(37800), -127/real(288), 41/real(180), 5/real(16), -2/real(3),
|
|
1/real(2),
|
|
-1983433/real(1935360), 281/real(630), 557/real(1440), -3/real(5),
|
|
13/real(48),
|
|
167603/real(181440), 15061/real(26880), -103/real(140), 61/real(240),
|
|
6601661/real(7257600), -179/real(168), 49561/real(161280),
|
|
-3418889/real(1995840), 34729/real(80640),
|
|
212378941/real(319334400),
|
|
// C[mu,xi]
|
|
12674323/real(851350500), -384229/real(14968800), -1609/real(28350),
|
|
121/real(1680), 4/real(45), -1/real(6),
|
|
-31621753811LL/real(1307674368000LL), -431/real(17325),
|
|
16463/real(453600), 26/real(945), -29/real(720),
|
|
-32844781/real(1751349600), 3746047/real(119750400), 449/real(28350),
|
|
-1003/real(45360),
|
|
10650637121LL/real(326918592000LL), 629/real(53460),
|
|
-40457/real(2419200),
|
|
205072597/real(20432412000LL), -1800439/real(119750400),
|
|
-59109051671LL/real(3923023104000LL),
|
|
// C[chi,phi]
|
|
4642/real(4725), 32/real(45), -82/real(45), 4/real(3), 2/real(3), -2,
|
|
-1522/real(945), 904/real(315), -13/real(9), -16/real(15), 5/real(3),
|
|
-12686/real(2835), 8/real(5), 34/real(21), -26/real(15),
|
|
-24832/real(14175), -12/real(5), 1237/real(630),
|
|
109598/real(31185), -734/real(315),
|
|
444337/real(155925),
|
|
// C[chi,beta]
|
|
-998/real(4725), 2/real(5), -16/real(45), 0, 2/real(3), -1,
|
|
-2/real(27), -22/real(105), 19/real(45), -2/real(5), 1/real(6),
|
|
116/real(567), -22/real(105), 16/real(105), -1/real(15),
|
|
2123/real(14175), -8/real(105), 17/real(1260),
|
|
128/real(4455), -1/real(105),
|
|
149/real(311850),
|
|
// C[chi,theta]
|
|
1042/real(4725), -14/real(45), -2/real(9), 2/real(3), 2/real(3), 0,
|
|
-712/real(945), -4/real(45), 43/real(45), 4/real(15), -1/real(3),
|
|
274/real(2835), 124/real(105), 2/real(105), -2/real(5),
|
|
21068/real(14175), -16/real(105), -55/real(126),
|
|
-9202/real(31185), -22/real(45),
|
|
-90263/real(155925),
|
|
// C[chi,mu]
|
|
-96199/real(604800), 81/real(512), 1/real(360), -37/real(96), 2/real(3),
|
|
-1/real(2),
|
|
1118711/real(3870720), -46/real(105), 437/real(1440), -1/real(15),
|
|
-1/real(48),
|
|
-5569/real(90720), 209/real(4480), 37/real(840), -17/real(480),
|
|
830251/real(7257600), 11/real(504), -4397/real(161280),
|
|
108847/real(3991680), -4583/real(161280),
|
|
-20648693/real(638668800),
|
|
// C[chi,chi] skipped
|
|
// C[chi,xi]
|
|
-55271278/real(212837625), 27128/real(93555), -2312/real(14175),
|
|
-88/real(315), 34/real(45), -2/real(3),
|
|
106691108/real(638512875), -65864/real(155925), 6079/real(14175),
|
|
-184/real(945), 1/real(45),
|
|
5921152/real(54729675), -14246/real(467775), 772/real(14175),
|
|
-106/real(2835),
|
|
75594328/real(638512875), -5312/real(467775), -167/real(9450),
|
|
2837636/real(638512875), -248/real(13365),
|
|
-34761247/real(1915538625),
|
|
// C[xi,phi]
|
|
-44732/real(2837835), 20824/real(467775), 538/real(4725), 88/real(315),
|
|
-4/real(45), -4/real(3),
|
|
-12467764/real(212837625), -37192/real(467775), -2482/real(14175),
|
|
8/real(105), 34/real(45),
|
|
100320856/real(1915538625), 54968/real(467775), -898/real(14175),
|
|
-1532/real(2835),
|
|
-5884124/real(70945875), 24496/real(467775), 6007/real(14175),
|
|
-839792/real(19348875), -23356/real(66825),
|
|
570284222/real(1915538625),
|
|
// C[xi,beta]
|
|
-70496/real(8513505), 2476/real(467775), 34/real(675), 32/real(315),
|
|
-4/real(45), -1/real(3),
|
|
53836/real(212837625), 3992/real(467775), 74/real(2025), -4/real(315),
|
|
-7/real(90),
|
|
-661844/real(1915538625), 7052/real(467775), 2/real(14175),
|
|
-83/real(2835),
|
|
1425778/real(212837625), 934/real(467775), -797/real(56700),
|
|
390088/real(212837625), -3673/real(467775),
|
|
-18623681/real(3831077250LL),
|
|
// C[xi,theta]
|
|
-4286228/real(42567525), -193082/real(467775), 778/real(4725),
|
|
62/real(105), -4/real(45), 2/real(3),
|
|
-61623938/real(70945875), 92696/real(467775), 12338/real(14175),
|
|
-32/real(315), 4/real(45),
|
|
427003576/real(1915538625), 612536/real(467775), -1618/real(14175),
|
|
-524/real(2835),
|
|
427770788/real(212837625), -8324/real(66825), -5933/real(14175),
|
|
-9153184/real(70945875), -320044/real(467775),
|
|
-1978771378/real(1915538625),
|
|
// C[xi,mu]
|
|
-9292991/real(302702400), 7764059/real(239500800), 1297/real(18900),
|
|
-817/real(10080), -4/real(45), 1/real(6),
|
|
36019108271LL/real(871782912000LL), 35474/real(467775),
|
|
-29609/real(453600), -2/real(35), 49/real(720),
|
|
3026004511LL/real(30648618000LL), -4306823/real(59875200),
|
|
-2917/real(56700), 4463/real(90720),
|
|
-368661577/real(4036032000LL), -102293/real(1871100),
|
|
331799/real(7257600),
|
|
-875457073/real(13621608000LL), 11744233/real(239500800),
|
|
453002260127LL/real(7846046208000LL),
|
|
// C[xi,chi]
|
|
2706758/real(42567525), -55222/real(93555), 2458/real(4725),
|
|
46/real(315), -34/real(45), 2/real(3),
|
|
-340492279/real(212837625), 516944/real(467775), 3413/real(14175),
|
|
-256/real(315), 19/real(45),
|
|
4430783356LL/real(1915538625), 206834/real(467775), -15958/real(14175),
|
|
248/real(567),
|
|
62016436/real(70945875), -832976/real(467775), 16049/real(28350),
|
|
-651151712/real(212837625), 15602/real(18711),
|
|
2561772812LL/real(1915538625),
|
|
// C[xi,xi] skipped
|
|
};
|
|
static const int ptrs[] = {
|
|
0, 0, 12, 24, 36, 57, 78, 90, 90, 102, 114, 135, 156, 168, 180, 180, 192,
|
|
213, 234, 246, 258, 270, 270, 291, 312, 333, 354, 375, 396, 396, 417,
|
|
438, 459, 480, 501, 522, 522,
|
|
};
|
|
#elif GEOGRAPHICLIB_AUXLATITUDE_ORDER == 8
|
|
static const real coeffs[] = {
|
|
// C[phi,phi] skipped
|
|
// C[phi,beta]; even coeffs only
|
|
0, 0, 0, 1,
|
|
0, 0, 0, 1/real(2),
|
|
0, 0, 1/real(3),
|
|
0, 0, 1/real(4),
|
|
0, 1/real(5),
|
|
0, 1/real(6),
|
|
1/real(7),
|
|
1/real(8),
|
|
// C[phi,theta]; even coeffs only
|
|
-2, 2, -2, 2,
|
|
-8, 6, -4, 2,
|
|
16, -8, 8/real(3),
|
|
40, -16, 4,
|
|
-32, 32/real(5),
|
|
-64, 32/real(3),
|
|
128/real(7),
|
|
32,
|
|
// C[phi,mu]; even coeffs only
|
|
-6607/real(24576), 269/real(512), -27/real(32), 3/real(2),
|
|
-155113/real(122880), 6759/real(4096), -55/real(32), 21/real(16),
|
|
87963/real(20480), -417/real(128), 151/real(96),
|
|
2514467/real(245760), -15543/real(2560), 1097/real(512),
|
|
-69119/real(6144), 8011/real(2560),
|
|
-5962461/real(286720), 293393/real(61440),
|
|
6459601/real(860160),
|
|
332287993/real(27525120),
|
|
// C[phi,chi]
|
|
189416/real(99225), 16822/real(4725), -2854/real(675), 26/real(45),
|
|
116/real(45), -2, -2/real(3), 2,
|
|
141514/real(8505), -31256/real(1575), 2323/real(945), 2704/real(315),
|
|
-227/real(45), -8/real(5), 7/real(3),
|
|
-2363828/real(31185), 98738/real(14175), 73814/real(2835),
|
|
-1262/real(105), -136/real(35), 56/real(15),
|
|
14416399/real(935550), 11763988/real(155925), -399572/real(14175),
|
|
-332/real(35), 4279/real(630),
|
|
258316372/real(1216215), -2046082/real(31185), -144838/real(6237),
|
|
4174/real(315),
|
|
-2155215124LL/real(14189175), -115444544/real(2027025),
|
|
601676/real(22275),
|
|
-170079376/real(1216215), 38341552/real(675675),
|
|
1383243703/real(11351340),
|
|
// C[phi,xi]
|
|
-1683291094/real(37574026875LL), 22947844/real(1915538625),
|
|
28112932/real(212837625), 60136/real(467775), -2582/real(14175),
|
|
-16/real(35), 4/real(45), 4/real(3),
|
|
-14351220203LL/real(488462349375LL), 1228352/real(3007125),
|
|
251310128/real(638512875), -21016/real(51975), -11966/real(14175),
|
|
152/real(945), 46/real(45),
|
|
505559334506LL/real(488462349375LL), 138128272/real(147349125),
|
|
-8797648/real(10945935), -94388/real(66825), 3802/real(14175),
|
|
3044/real(2835),
|
|
973080708361LL/real(488462349375LL), -45079184/real(29469825),
|
|
-1472637812/real(638512875), 41072/real(93555), 6059/real(4725),
|
|
-1385645336626LL/real(488462349375LL), -550000184/real(147349125),
|
|
455935736/real(638512875), 768272/real(467775),
|
|
-2939205114427LL/real(488462349375LL), 443810768/real(383107725),
|
|
4210684958LL/real(1915538625),
|
|
101885255158LL/real(54273594375LL), 387227992/real(127702575),
|
|
1392441148867LL/real(325641566250LL),
|
|
// C[beta,phi]; even coeffs only
|
|
0, 0, 0, -1,
|
|
0, 0, 0, 1/real(2),
|
|
0, 0, -1/real(3),
|
|
0, 0, 1/real(4),
|
|
0, -1/real(5),
|
|
0, 1/real(6),
|
|
-1/real(7),
|
|
1/real(8),
|
|
// C[beta,beta] skipped
|
|
// C[beta,theta]; even coeffs only
|
|
0, 0, 0, 1,
|
|
0, 0, 0, 1/real(2),
|
|
0, 0, 1/real(3),
|
|
0, 0, 1/real(4),
|
|
0, 1/real(5),
|
|
0, 1/real(6),
|
|
1/real(7),
|
|
1/real(8),
|
|
// C[beta,mu]; even coeffs only
|
|
-4879/real(73728), 205/real(1536), -9/real(32), 1/real(2),
|
|
-86171/real(368640), 1335/real(4096), -37/real(96), 5/real(16),
|
|
2901/real(4096), -75/real(128), 29/real(96),
|
|
1082857/real(737280), -2391/real(2560), 539/real(1536),
|
|
-28223/real(18432), 3467/real(7680),
|
|
-733437/real(286720), 38081/real(61440),
|
|
459485/real(516096),
|
|
109167851/real(82575360),
|
|
// C[beta,chi]
|
|
-25666/real(99225), 4769/real(4725), -3118/real(4725), -1/real(3),
|
|
38/real(45), -1/real(3), -2/real(3), 1,
|
|
193931/real(42525), -14404/real(4725), -247/real(270), 50/real(21),
|
|
-7/real(9), -14/real(15), 5/real(6),
|
|
-1709614/real(155925), -36521/real(14175), 17564/real(2835), -5/real(3),
|
|
-34/real(21), 16/real(15),
|
|
-637699/real(85050), 2454416/real(155925), -49877/real(14175),
|
|
-28/real(9), 2069/real(1260),
|
|
48124558/real(1216215), -20989/real(2835), -28244/real(4455),
|
|
883/real(315),
|
|
-16969807/real(1091475), -2471888/real(184275), 797222/real(155925),
|
|
-1238578/real(42525), 2199332/real(225225),
|
|
87600385/real(4540536),
|
|
// C[beta,xi]
|
|
-5946082372LL/real(488462349375LL), 9708931/real(1915538625),
|
|
7947332/real(212837625), 11824/real(467775), -1082/real(14175),
|
|
-46/real(315), 4/real(45), 1/real(3),
|
|
190673521/real(69780335625LL), 164328266/real(1915538625),
|
|
39946703/real(638512875), -16672/real(155925), -338/real(2025),
|
|
68/real(945), 17/real(90),
|
|
86402898356LL/real(488462349375LL), 236067184/real(1915538625),
|
|
-255454/real(1563705), -101069/real(467775), 1102/real(14175),
|
|
461/real(2835),
|
|
110123070361LL/real(488462349375LL), -98401826/real(383107725),
|
|
-189032762/real(638512875), 1786/real(18711), 3161/real(18900),
|
|
-200020620676LL/real(488462349375LL), -802887278/real(1915538625),
|
|
80274086/real(638512875), 88868/real(467775),
|
|
-296107325077LL/real(488462349375LL), 66263486/real(383107725),
|
|
880980241/real(3831077250LL),
|
|
4433064236LL/real(18091198125LL), 37151038/real(127702575),
|
|
495248998393LL/real(1302566265000LL),
|
|
// C[theta,phi]; even coeffs only
|
|
2, -2, 2, -2,
|
|
-8, 6, -4, 2,
|
|
-16, 8, -8/real(3),
|
|
40, -16, 4,
|
|
32, -32/real(5),
|
|
-64, 32/real(3),
|
|
-128/real(7),
|
|
32,
|
|
// C[theta,beta]; even coeffs only
|
|
0, 0, 0, -1,
|
|
0, 0, 0, 1/real(2),
|
|
0, 0, -1/real(3),
|
|
0, 0, 1/real(4),
|
|
0, -1/real(5),
|
|
0, 1/real(6),
|
|
-1/real(7),
|
|
1/real(8),
|
|
// C[theta,theta] skipped
|
|
// C[theta,mu]; even coeffs only
|
|
-14321/real(73728), 499/real(1536), -23/real(32), -1/real(2),
|
|
-201467/real(368640), 6565/real(12288), -5/real(96), 5/real(16),
|
|
2939/real(4096), -77/real(128), 1/real(32),
|
|
1155049/real(737280), -4037/real(7680), 283/real(1536),
|
|
-19465/real(18432), 1301/real(7680),
|
|
-442269/real(286720), 17089/real(61440),
|
|
198115/real(516096),
|
|
48689387/real(82575360),
|
|
// C[theta,chi]
|
|
64424/real(99225), 76/real(225), -3658/real(4725), 2/real(9), 4/real(9),
|
|
-2/real(3), -2/real(3), 0,
|
|
2146/real(1215), -2728/real(945), 61/real(135), 68/real(45),
|
|
-23/real(45), -4/real(15), 1/real(3),
|
|
-95948/real(10395), 428/real(945), 9446/real(2835), -46/real(35),
|
|
-24/real(35), 2/real(5),
|
|
29741/real(85050), 4472/real(525), -34712/real(14175), -80/real(63),
|
|
83/real(126),
|
|
280108/real(13365), -17432/real(3465), -2362/real(891), 52/real(45),
|
|
-48965632/real(4729725), -548752/real(96525), 335882/real(155925),
|
|
-197456/real(15795), 51368/real(12285),
|
|
1461335/real(174636),
|
|
// C[theta,xi]
|
|
-230886326/real(6343666875LL), -189115382/real(1915538625),
|
|
216932/real(2627625), 109042/real(467775), -2102/real(14175),
|
|
-158/real(315), 4/real(45), -2/real(3),
|
|
-11696145869LL/real(69780335625LL), 288456008/real(1915538625),
|
|
117952358/real(638512875), -7256/real(155925), 934/real(14175),
|
|
-16/real(945), 16/real(45),
|
|
91546732346LL/real(488462349375LL), 478700902/real(1915538625),
|
|
-7391576/real(54729675), -25286/real(66825), 922/real(14175),
|
|
-232/real(2835),
|
|
218929662961LL/real(488462349375LL), -67330724/real(383107725),
|
|
-67048172/real(638512875), 268/real(18711), 719/real(4725),
|
|
-129039188386LL/real(488462349375LL), -117954842/real(273648375),
|
|
46774256/real(638512875), 14354/real(467775),
|
|
-178084928947LL/real(488462349375LL), 2114368/real(34827975),
|
|
253129538/real(1915538625),
|
|
6489189398LL/real(54273594375LL), 13805944/real(127702575),
|
|
59983985827LL/real(325641566250LL),
|
|
// C[mu,phi]; even coeffs only
|
|
57/real(2048), -3/real(32), 9/real(16), -3/real(2),
|
|
-105/real(4096), 135/real(2048), -15/real(32), 15/real(16),
|
|
-105/real(2048), 105/real(256), -35/real(48),
|
|
693/real(16384), -189/real(512), 315/real(512),
|
|
693/real(2048), -693/real(1280),
|
|
-1287/real(4096), 1001/real(2048),
|
|
-6435/real(14336),
|
|
109395/real(262144),
|
|
// C[mu,beta]; even coeffs only
|
|
19/real(2048), -1/real(32), 3/real(16), -1/real(2),
|
|
7/real(4096), -9/real(2048), 1/real(32), -1/real(16),
|
|
-3/real(2048), 3/real(256), -1/real(48),
|
|
-11/real(16384), 3/real(512), -5/real(512),
|
|
7/real(2048), -7/real(1280),
|
|
9/real(4096), -7/real(2048),
|
|
-33/real(14336),
|
|
-429/real(262144),
|
|
// C[mu,theta]; even coeffs only
|
|
509/real(2048), -15/real(32), 13/real(16), 1/real(2),
|
|
2599/real(4096), -1673/real(2048), 33/real(32), -1/real(16),
|
|
-2989/real(2048), 349/real(256), -5/real(16),
|
|
-43531/real(16384), 963/real(512), -261/real(512),
|
|
5545/real(2048), -921/real(1280),
|
|
16617/real(4096), -6037/real(6144),
|
|
-19279/real(14336),
|
|
-490925/real(262144),
|
|
// C[mu,mu] skipped
|
|
// C[mu,chi]
|
|
-18975107/real(50803200), 72161/real(387072), 7891/real(37800),
|
|
-127/real(288), 41/real(180), 5/real(16), -2/real(3), 1/real(2),
|
|
148003883/real(174182400), 13769/real(28800), -1983433/real(1935360),
|
|
281/real(630), 557/real(1440), -3/real(5), 13/real(48),
|
|
79682431/real(79833600), -67102379/real(29030400), 167603/real(181440),
|
|
15061/real(26880), -103/real(140), 61/real(240),
|
|
-40176129013LL/real(7664025600LL), 97445/real(49896),
|
|
6601661/real(7257600), -179/real(168), 49561/real(161280),
|
|
2605413599LL/real(622702080), 14644087/real(9123840),
|
|
-3418889/real(1995840), 34729/real(80640),
|
|
175214326799LL/real(58118860800LL), -30705481/real(10378368),
|
|
212378941/real(319334400),
|
|
-16759934899LL/real(3113510400LL), 1522256789/real(1383782400),
|
|
1424729850961LL/real(743921418240LL),
|
|
// C[mu,xi]
|
|
-375027460897LL/real(125046361440000LL),
|
|
7183403063LL/real(560431872000LL), 12674323/real(851350500),
|
|
-384229/real(14968800), -1609/real(28350), 121/real(1680), 4/real(45),
|
|
-1/real(6),
|
|
30410873385097LL/real(2000741783040000LL),
|
|
1117820213/real(122594472000LL), -31621753811LL/real(1307674368000LL),
|
|
-431/real(17325), 16463/real(453600), 26/real(945), -29/real(720),
|
|
151567502183LL/real(17863765920000LL),
|
|
-116359346641LL/real(3923023104000LL), -32844781/real(1751349600),
|
|
3746047/real(119750400), 449/real(28350), -1003/real(45360),
|
|
-317251099510901LL/real(8002967132160000LL), -13060303/real(766215450),
|
|
10650637121LL/real(326918592000LL), 629/real(53460),
|
|
-40457/real(2419200),
|
|
-2105440822861LL/real(125046361440000LL),
|
|
146875240637LL/real(3923023104000LL), 205072597/real(20432412000LL),
|
|
-1800439/real(119750400),
|
|
91496147778023LL/real(2000741783040000LL), 228253559/real(24518894400LL),
|
|
-59109051671LL/real(3923023104000LL),
|
|
126430355893LL/real(13894040160000LL),
|
|
-4255034947LL/real(261534873600LL),
|
|
-791820407649841LL/real(42682491371520000LL),
|
|
// C[chi,phi]
|
|
1514/real(1323), -8384/real(4725), 4642/real(4725), 32/real(45),
|
|
-82/real(45), 4/real(3), 2/real(3), -2,
|
|
142607/real(42525), -2288/real(1575), -1522/real(945), 904/real(315),
|
|
-13/real(9), -16/real(15), 5/real(3),
|
|
120202/real(51975), 44644/real(14175), -12686/real(2835), 8/real(5),
|
|
34/real(21), -26/real(15),
|
|
-1097407/real(187110), 1077964/real(155925), -24832/real(14175),
|
|
-12/real(5), 1237/real(630),
|
|
-12870194/real(1216215), 1040/real(567), 109598/real(31185),
|
|
-734/real(315),
|
|
-126463/real(72765), -941912/real(184275), 444337/real(155925),
|
|
3463678/real(467775), -2405834/real(675675),
|
|
256663081/real(56756700),
|
|
// C[chi,beta]
|
|
1384/real(11025), -34/real(4725), -998/real(4725), 2/real(5),
|
|
-16/real(45), 0, 2/real(3), -1,
|
|
-12616/real(42525), 1268/real(4725), -2/real(27), -22/real(105),
|
|
19/real(45), -2/real(5), 1/real(6),
|
|
1724/real(51975), -1858/real(14175), 116/real(567), -22/real(105),
|
|
16/real(105), -1/real(15),
|
|
115249/real(935550), -26836/real(155925), 2123/real(14175), -8/real(105),
|
|
17/real(1260),
|
|
140836/real(1216215), -424/real(6237), 128/real(4455), -1/real(105),
|
|
210152/real(4729725), -31232/real(2027025), 149/real(311850),
|
|
30208/real(6081075), -499/real(225225),
|
|
-68251/real(113513400),
|
|
// C[chi,theta]
|
|
-1738/real(11025), 18/real(175), 1042/real(4725), -14/real(45),
|
|
-2/real(9), 2/real(3), 2/real(3), 0,
|
|
23159/real(42525), 332/real(945), -712/real(945), -4/real(45),
|
|
43/real(45), 4/real(15), -1/real(3),
|
|
13102/real(31185), -1352/real(945), 274/real(2835), 124/real(105),
|
|
2/real(105), -2/real(5),
|
|
-2414843/real(935550), 1528/real(4725), 21068/real(14175), -16/real(105),
|
|
-55/real(126),
|
|
60334/real(93555), 20704/real(10395), -9202/real(31185), -22/real(45),
|
|
40458083/real(14189175), -299444/real(675675), -90263/real(155925),
|
|
-3818498/real(6081075), -8962/real(12285),
|
|
-4259027/real(4365900),
|
|
// C[chi,mu]
|
|
-7944359/real(67737600), 5406467/real(38707200), -96199/real(604800),
|
|
81/real(512), 1/real(360), -37/real(96), 2/real(3), -1/real(2),
|
|
-24749483/real(348364800), -51841/real(1209600), 1118711/real(3870720),
|
|
-46/real(105), 437/real(1440), -1/real(15), -1/real(48),
|
|
6457463/real(17740800), -9261899/real(58060800), -5569/real(90720),
|
|
209/real(4480), 37/real(840), -17/real(480),
|
|
-324154477/real(7664025600LL), -466511/real(2494800),
|
|
830251/real(7257600), 11/real(504), -4397/real(161280),
|
|
-22894433/real(124540416), 8005831/real(63866880), 108847/real(3991680),
|
|
-4583/real(161280),
|
|
2204645983LL/real(12915302400LL), 16363163/real(518918400),
|
|
-20648693/real(638668800),
|
|
497323811/real(12454041600LL), -219941297/real(5535129600LL),
|
|
-191773887257LL/real(3719607091200LL),
|
|
// C[chi,chi] skipped
|
|
// C[chi,xi]
|
|
-17451293242LL/real(488462349375LL), 308365186/real(1915538625),
|
|
-55271278/real(212837625), 27128/real(93555), -2312/real(14175),
|
|
-88/real(315), 34/real(45), -2/real(3),
|
|
-101520127208LL/real(488462349375LL), 149984636/real(1915538625),
|
|
106691108/real(638512875), -65864/real(155925), 6079/real(14175),
|
|
-184/real(945), 1/real(45),
|
|
10010741462LL/real(37574026875LL), -99534832/real(383107725),
|
|
5921152/real(54729675), -14246/real(467775), 772/real(14175),
|
|
-106/real(2835),
|
|
1615002539/real(75148053750LL), -35573728/real(273648375),
|
|
75594328/real(638512875), -5312/real(467775), -167/real(9450),
|
|
-3358119706LL/real(488462349375LL), 130601488/real(1915538625),
|
|
2837636/real(638512875), -248/real(13365),
|
|
46771947158LL/real(488462349375LL), -3196/real(3553875),
|
|
-34761247/real(1915538625),
|
|
-18696014/real(18091198125LL), -2530364/real(127702575),
|
|
-14744861191LL/real(651283132500LL),
|
|
// C[xi,phi]
|
|
-88002076/real(13956067125LL), -86728/real(16372125),
|
|
-44732/real(2837835), 20824/real(467775), 538/real(4725), 88/real(315),
|
|
-4/real(45), -4/real(3),
|
|
-2641983469LL/real(488462349375LL), -895712/real(147349125),
|
|
-12467764/real(212837625), -37192/real(467775), -2482/real(14175),
|
|
8/real(105), 34/real(45),
|
|
8457703444LL/real(488462349375LL), 240616/real(4209975),
|
|
100320856/real(1915538625), 54968/real(467775), -898/real(14175),
|
|
-1532/real(2835),
|
|
-4910552477LL/real(97692469875LL), -4832848/real(147349125),
|
|
-5884124/real(70945875), 24496/real(467775), 6007/real(14175),
|
|
9393713176LL/real(488462349375LL), 816824/real(13395375),
|
|
-839792/real(19348875), -23356/real(66825),
|
|
-4532926649LL/real(97692469875LL), 1980656/real(54729675),
|
|
570284222/real(1915538625),
|
|
-14848113968LL/real(488462349375LL), -496894276/real(1915538625),
|
|
224557742191LL/real(976924698750LL),
|
|
// C[xi,beta]
|
|
29232878/real(97692469875LL), -18484/real(4343625), -70496/real(8513505),
|
|
2476/real(467775), 34/real(675), 32/real(315), -4/real(45), -1/real(3),
|
|
-324943819/real(488462349375LL), -4160804/real(1915538625),
|
|
53836/real(212837625), 3992/real(467775), 74/real(2025), -4/real(315),
|
|
-7/real(90),
|
|
-168643106/real(488462349375LL), 237052/real(383107725),
|
|
-661844/real(1915538625), 7052/real(467775), 2/real(14175),
|
|
-83/real(2835),
|
|
113042383/real(97692469875LL), -2915326/real(1915538625),
|
|
1425778/real(212837625), 934/real(467775), -797/real(56700),
|
|
-558526274/real(488462349375LL), 6064888/real(1915538625),
|
|
390088/real(212837625), -3673/real(467775),
|
|
155665021/real(97692469875LL), 41288/real(29469825),
|
|
-18623681/real(3831077250LL),
|
|
504234982/real(488462349375LL), -6205669/real(1915538625),
|
|
-8913001661LL/real(3907698795000LL),
|
|
// C[xi,theta]
|
|
182466964/real(8881133625LL), 53702182/real(212837625),
|
|
-4286228/real(42567525), -193082/real(467775), 778/real(4725),
|
|
62/real(105), -4/real(45), 2/real(3),
|
|
367082779691LL/real(488462349375LL), -32500616/real(273648375),
|
|
-61623938/real(70945875), 92696/real(467775), 12338/real(14175),
|
|
-32/real(315), 4/real(45),
|
|
-42668482796LL/real(488462349375LL), -663111728/real(383107725),
|
|
427003576/real(1915538625), 612536/real(467775), -1618/real(14175),
|
|
-524/real(2835),
|
|
-327791986997LL/real(97692469875LL), 421877252/real(1915538625),
|
|
427770788/real(212837625), -8324/real(66825), -5933/real(14175),
|
|
74612072536LL/real(488462349375LL), 6024982024LL/real(1915538625),
|
|
-9153184/real(70945875), -320044/real(467775),
|
|
489898512247LL/real(97692469875LL), -46140784/real(383107725),
|
|
-1978771378/real(1915538625),
|
|
-42056042768LL/real(488462349375LL), -2926201612LL/real(1915538625),
|
|
-2209250801969LL/real(976924698750LL),
|
|
// C[xi,mu]
|
|
39534358147LL/real(2858202547200LL),
|
|
-25359310709LL/real(1743565824000LL), -9292991/real(302702400),
|
|
7764059/real(239500800), 1297/real(18900), -817/real(10080), -4/real(45),
|
|
1/real(6),
|
|
-13216941177599LL/real(571640509440000LL),
|
|
-14814966289LL/real(245188944000LL), 36019108271LL/real(871782912000LL),
|
|
35474/real(467775), -29609/real(453600), -2/real(35), 49/real(720),
|
|
-27782109847927LL/real(250092722880000LL),
|
|
99871724539LL/real(1569209241600LL), 3026004511LL/real(30648618000LL),
|
|
-4306823/real(59875200), -2917/real(56700), 4463/real(90720),
|
|
168979300892599LL/real(1600593426432000LL),
|
|
2123926699/real(15324309000LL), -368661577/real(4036032000LL),
|
|
-102293/real(1871100), 331799/real(7257600),
|
|
1959350112697LL/real(9618950880000LL),
|
|
-493031379277LL/real(3923023104000LL), -875457073/real(13621608000LL),
|
|
11744233/real(239500800),
|
|
-145659994071373LL/real(800296713216000LL),
|
|
-793693009/real(9807557760LL), 453002260127LL/real(7846046208000LL),
|
|
-53583096419057LL/real(500185445760000LL),
|
|
103558761539LL/real(1426553856000LL),
|
|
12272105438887727LL/real(128047474114560000LL),
|
|
// C[xi,chi]
|
|
-64724382148LL/real(97692469875LL), 16676974/real(30405375),
|
|
2706758/real(42567525), -55222/real(93555), 2458/real(4725),
|
|
46/real(315), -34/real(45), 2/real(3),
|
|
85904355287LL/real(37574026875LL), 158999572/real(1915538625),
|
|
-340492279/real(212837625), 516944/real(467775), 3413/real(14175),
|
|
-256/real(315), 19/real(45),
|
|
2986003168LL/real(37574026875LL), -7597644214LL/real(1915538625),
|
|
4430783356LL/real(1915538625), 206834/real(467775), -15958/real(14175),
|
|
248/real(567),
|
|
-375566203/real(39037950), 851209552/real(174139875),
|
|
62016436/real(70945875), -832976/real(467775), 16049/real(28350),
|
|
5106181018156LL/real(488462349375LL), 3475643362LL/real(1915538625),
|
|
-651151712/real(212837625), 15602/real(18711),
|
|
34581190223LL/real(8881133625LL), -10656173804LL/real(1915538625),
|
|
2561772812LL/real(1915538625),
|
|
-5150169424688LL/real(488462349375LL), 873037408/real(383107725),
|
|
7939103697617LL/real(1953849397500LL),
|
|
// C[xi,xi] skipped
|
|
};
|
|
static const int ptrs[] = {
|
|
0, 0, 20, 40, 60, 96, 132, 152, 152, 172, 192, 228, 264, 284, 304, 304,
|
|
324, 360, 396, 416, 436, 456, 456, 492, 528, 564, 600, 636, 672, 672,
|
|
708, 744, 780, 816, 852, 888, 888,
|
|
};
|
|
#else
|
|
#error "Unsupported value for GEOGRAPHICLIB_AUXLATITUDE_ORDER"
|
|
#endif
|
|
|
|
static_assert(sizeof(ptrs) / sizeof(int) == AUXNUMBER*AUXNUMBER+1,
|
|
"Mismatch in size of ptrs array");
|
|
static_assert(sizeof(coeffs) / sizeof(real) ==
|
|
(RECTIFYING+1)*RECTIFYING *
|
|
(Lmax * (Lmax + 3) - 2*(Lmax/2))/4 // Even only arrays
|
|
+ (AUXNUMBER*(AUXNUMBER-1) - (RECTIFYING+1)*RECTIFYING) *
|
|
(Lmax * (Lmax + 1))/2,
|
|
"Mismatch in size of coeffs array");
|
|
|
|
if (k < 0) return; // auxout or auxin out of range
|
|
if (auxout == auxin)
|
|
fill(_c + Lmax * k, _c + Lmax * (k + 1), 0);
|
|
else {
|
|
int o = ptrs[k];
|
|
real d = _n;
|
|
if (auxin <= RECTIFYING && auxout <= RECTIFYING) {
|
|
for (int l = 0; l < Lmax; ++l) {
|
|
int m = (Lmax - l - 1) / 2; // order of polynomial in n^2
|
|
_c[Lmax * k + l] = d * Math::polyval(m, coeffs + o, _n2);
|
|
o += m + 1;
|
|
d *= _n;
|
|
}
|
|
} else {
|
|
for (int l = 0; l < Lmax; ++l) {
|
|
int m = (Lmax - l - 1); // order of polynomial in n
|
|
_c[Lmax * k + l] = d * Math::polyval(m, coeffs + o, _n);
|
|
o += m + 1;
|
|
d *= _n;
|
|
}
|
|
}
|
|
// assert (o == ptrs[AUXNUMBER * auxout + auxin + 1])
|
|
}
|
|
}
|
|
|
|
template<typename T>
|
|
typename AuxLatitude<T>::real
|
|
AuxLatitude<T>::Clenshaw(real szeta, real czeta,
|
|
const real c[], int K, bool alt) {
|
|
// Evaluate
|
|
// y = sum(c[k] * sin( (2*k+2) * zeta), i, 0, K-1)
|
|
// Approx operation count = (K + 5) mult and (2 * K + 2) add
|
|
// If alt, use the Reinsch optimizations for small szeta and small czeta.
|
|
int k = K;
|
|
real u0 = 0; // accumulator for sum
|
|
if (alt && 2*fabs(szeta) < 1) {
|
|
real X = -4 * szeta * szeta,
|
|
d0 = 0; // other accumulator for sum
|
|
for (; k > 0;) {
|
|
d0 = X * u0 + d0 + c[--k];
|
|
u0 = d0 + u0;
|
|
}
|
|
} else if (alt && 2*fabs(czeta) < 1) {
|
|
real X = 4 * czeta * czeta,
|
|
d0 = 0; // other accumulator for sum
|
|
for (; k > 0;) {
|
|
d0 = X * u0 - d0 + c[--k];
|
|
u0 = d0 - u0;
|
|
}
|
|
} else {
|
|
real X = 2 * (czeta - szeta) * (czeta + szeta), // 2 * cos(2*zeta)
|
|
u1 = 0; // other accumulator for sum
|
|
for (; k > 0;) {
|
|
real t = X * u0 - u1 + c[--k];
|
|
u1 = u0; u0 = t;
|
|
}
|
|
}
|
|
return 2 * szeta * czeta * u0; // sin(2*zeta) * u0
|
|
}
|
|
|
|
template class AuxAngle<Math::real>;
|
|
template class AuxLatitude<Math::real>;
|
|
#if GEOGRAPHICLIB_PRECISION != 2
|
|
template class AuxAngle<double>;
|
|
template class AuxLatitude<double>;
|
|
#endif
|
|
|
|
} // namespace GeographicLib
|