00001 /** 00002 * \file EllipticFunction.hpp 00003 * \brief Header for GeographicLib::EllipticFunction class 00004 * 00005 * Copyright (c) Charles Karney (2008, 2009) <charles@karney.com> 00006 * and licensed under the LGPL. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP) 00011 #define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP "$Id: EllipticFunction.hpp 6720 2009-10-17 23:13:57Z ckarney $" 00012 00013 #include "GeographicLib/Constants.hpp" 00014 00015 namespace GeographicLib { 00016 00017 /** 00018 * \brief Elliptic functions needed for TransverseMercatorExact 00019 * 00020 * This provides the subset of elliptic functions needed for 00021 * TransverseMercatorExact. For a given ellipsoid, only parameters \e 00022 * e<sup>2</sup> and 1 - \e e<sup>2</sup> are needed. This class taken the 00023 * parameter as a constructor parameters and caches the values of the 00024 * required complete integrals. A method is provided for Jacobi elliptic 00025 * functions and for the incomplete elliptic integral of the second kind in 00026 * terms of the amplitude. 00027 * 00028 * The computation of the elliptic integrals uses the algorithms given in 00029 * - B. C. Carlson, 00030 * <a href="http://dx.doi.org/10.1007/BF02198293"> Computation of elliptic 00031 * integrals</a>, Numerical Algorithms 10, 13–26 (1995). 00032 * . 00033 * The computation of the Jacobi elliptic functions uses the algorithm given 00034 * in 00035 * - R. Bulirsch, 00036 * <a href="http://dx.doi.org/10.1007/BF01397975"> Numerical Calculation of 00037 * Elliptic Integrals and Elliptic Functions</a>, Numericshe Mathematik 7, 00038 * 78–90 (1965). 00039 * . 00040 * The notation follows Abramowitz and Stegun, Chapters 16 and 17. 00041 **********************************************************************/ 00042 class EllipticFunction { 00043 private: 00044 typedef Math::real real; 00045 static const real tol, tolRF, tolRD, tolRG0, tolJAC, tolJAC1; 00046 enum { num = 10 }; // Max depth required for sncndn. Probably 5 is enough. 00047 static real RF(real x, real y, real z) throw(); 00048 static real RD(real x, real y, real z) throw(); 00049 static real RG0(real x, real y) throw(); 00050 const real _m, _m1; 00051 mutable bool _init; 00052 mutable real _kc, _ec, _kec; 00053 bool Init() const throw(); 00054 public: 00055 00056 /** 00057 * Constructor with parameter \e m. 00058 **********************************************************************/ 00059 explicit EllipticFunction(real m) throw(); 00060 00061 /** 00062 * The parameter \e m. 00063 **********************************************************************/ 00064 Math::real m() const throw() { return _m; } 00065 00066 /** 00067 * The complementary parameter \e m' = (1 - \e m). 00068 **********************************************************************/ 00069 Math::real m1() const throw() { return _m1; } 00070 00071 /** 00072 * The complete integral of first kind, \e K(\e m). 00073 **********************************************************************/ 00074 Math::real K() const throw() { _init || Init(); return _kc; } 00075 00076 /** 00077 * The complete integral of second kind, \e E(\e m). 00078 **********************************************************************/ 00079 Math::real E() const throw() { _init || Init(); return _ec; } 00080 00081 /** 00082 * The difference \e K(\e m) - \e E(\e m) (which can be computed directly). 00083 **********************************************************************/ 00084 Math::real KE() const throw() { _init || Init(); return _kec; } 00085 00086 /** 00087 * The Jacobi elliptic functions sn(<i>x</i>|<i>m</i>), 00088 * cn(<i>x</i>|<i>m</i>), and dn(<i>x</i>|<i>m</i>) with argument \e x. 00089 * The results are returned in \e sn, \e cn, and \e dn. 00090 **********************************************************************/ 00091 void sncndn(real x, real& sn, real& cn, real& dn) const throw(); 00092 00093 /** 00094 * The incomplete integral of the second kind = int dn(\e w)<sup>2</sup> \e 00095 * dw (A+S 17.2.10). Instead of specifying the ampltiude \e phi, we 00096 * provide \e sn = sin(\e phi), \e cn = cos(\e phi), \e dn = sqrt(1 - \e m 00097 * sin<sup>2</sup>(\e phi)). 00098 **********************************************************************/ 00099 Math::real E(real sn, real cn, real dn) const throw(); 00100 }; 00101 00102 00103 } // namespace GeographicLib 00104 00105 #endif