ell_integral.tcc

Go to the documentation of this file.
00001 // Special functions -*- C++ -*-
00002 
00003 // Copyright (C) 2006, 2007, 2008, 2009
00004 // Free Software Foundation, Inc.
00005 //
00006 // This file is part of the GNU ISO C++ Library.  This library is free
00007 // software; you can redistribute it and/or modify it under the
00008 // terms of the GNU General Public License as published by the
00009 // Free Software Foundation; either version 3, or (at your option)
00010 // any later version.
00011 //
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 //
00017 // Under Section 7 of GPL version 3, you are granted additional
00018 // permissions described in the GCC Runtime Library Exception, version
00019 // 3.1, as published by the Free Software Foundation.
00020 
00021 // You should have received a copy of the GNU General Public License and
00022 // a copy of the GCC Runtime Library Exception along with this program;
00023 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
00024 // <http://www.gnu.org/licenses/>.
00025 
00026 /** @file tr1/ell_integral.tcc
00027  *  This is an internal header file, included by other library headers.
00028  *  You should not attempt to use it directly.
00029  */
00030 
00031 //
00032 // ISO C++ 14882 TR1: 5.2  Special functions
00033 //
00034 
00035 // Written by Edward Smith-Rowland based on:
00036 //   (1)  B. C. Carlson Numer. Math. 33, 1 (1979)
00037 //   (2)  B. C. Carlson, Special Functions of Applied Mathematics (1977)
00038 //   (3)  The Gnu Scientific Library, http://www.gnu.org/software/gsl
00039 //   (4)  Numerical Recipes in C, 2nd ed, by W. H. Press, S. A. Teukolsky,
00040 //        W. T. Vetterling, B. P. Flannery, Cambridge University Press
00041 //        (1992), pp. 261-269
00042 
00043 #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC
00044 #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1
00045 
00046 namespace std
00047 {
00048 namespace tr1
00049 {
00050 
00051   // [5.2] Special functions
00052 
00053   // Implementation-space details.
00054   namespace __detail
00055   {
00056 
00057     /**
00058      *   @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$
00059      *          of the first kind.
00060      * 
00061      *   The Carlson elliptic function of the first kind is defined by:
00062      *   @f[
00063      *       R_F(x,y,z) = \frac{1}{2} \int_0^\infty
00064      *                 \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}}
00065      *   @f]
00066      *
00067      *   @param  __x  The first of three symmetric arguments.
00068      *   @param  __y  The second of three symmetric arguments.
00069      *   @param  __z  The third of three symmetric arguments.
00070      *   @return  The Carlson elliptic function of the first kind.
00071      */
00072     template<typename _Tp>
00073     _Tp
00074     __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z)
00075     {
00076       const _Tp __min = std::numeric_limits<_Tp>::min();
00077       const _Tp __max = std::numeric_limits<_Tp>::max();
00078       const _Tp __lolim = _Tp(5) * __min;
00079       const _Tp __uplim = __max / _Tp(5);
00080 
00081       if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
00082         std::__throw_domain_error(__N("Argument less than zero "
00083                                       "in __ellint_rf."));
00084       else if (__x + __y < __lolim || __x + __z < __lolim
00085             || __y + __z < __lolim)
00086         std::__throw_domain_error(__N("Argument too small in __ellint_rf"));
00087       else
00088         {
00089           const _Tp __c0 = _Tp(1) / _Tp(4);
00090           const _Tp __c1 = _Tp(1) / _Tp(24);
00091           const _Tp __c2 = _Tp(1) / _Tp(10);
00092           const _Tp __c3 = _Tp(3) / _Tp(44);
00093           const _Tp __c4 = _Tp(1) / _Tp(14);
00094 
00095           _Tp __xn = __x;
00096           _Tp __yn = __y;
00097           _Tp __zn = __z;
00098 
00099           const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00100           const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6));
00101           _Tp __mu;
00102           _Tp __xndev, __yndev, __zndev;
00103 
00104           const unsigned int __max_iter = 100;
00105           for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00106             {
00107               __mu = (__xn + __yn + __zn) / _Tp(3);
00108               __xndev = 2 - (__mu + __xn) / __mu;
00109               __yndev = 2 - (__mu + __yn) / __mu;
00110               __zndev = 2 - (__mu + __zn) / __mu;
00111               _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
00112               __epsilon = std::max(__epsilon, std::abs(__zndev));
00113               if (__epsilon < __errtol)
00114                 break;
00115               const _Tp __xnroot = std::sqrt(__xn);
00116               const _Tp __ynroot = std::sqrt(__yn);
00117               const _Tp __znroot = std::sqrt(__zn);
00118               const _Tp __lambda = __xnroot * (__ynroot + __znroot)
00119                                  + __ynroot * __znroot;
00120               __xn = __c0 * (__xn + __lambda);
00121               __yn = __c0 * (__yn + __lambda);
00122               __zn = __c0 * (__zn + __lambda);
00123             }
00124 
00125           const _Tp __e2 = __xndev * __yndev - __zndev * __zndev;
00126           const _Tp __e3 = __xndev * __yndev * __zndev;
00127           const _Tp __s  = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2
00128                    + __c4 * __e3;
00129 
00130           return __s / std::sqrt(__mu);
00131         }
00132     }
00133 
00134 
00135     /**
00136      *   @brief Return the complete elliptic integral of the first kind
00137      *          @f$ K(k) @f$ by series expansion.
00138      * 
00139      *   The complete elliptic integral of the first kind is defined as
00140      *   @f[
00141      *     K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
00142      *                              {\sqrt{1 - k^2sin^2\theta}}
00143      *   @f]
00144      * 
00145      *   This routine is not bad as long as |k| is somewhat smaller than 1
00146      *   but is not is good as the Carlson elliptic integral formulation.
00147      * 
00148      *   @param  __k  The argument of the complete elliptic function.
00149      *   @return  The complete elliptic function of the first kind.
00150      */
00151     template<typename _Tp>
00152     _Tp
00153     __comp_ellint_1_series(const _Tp __k)
00154     {
00155 
00156       const _Tp __kk = __k * __k;
00157 
00158       _Tp __term = __kk / _Tp(4);
00159       _Tp __sum = _Tp(1) + __term;
00160 
00161       const unsigned int __max_iter = 1000;
00162       for (unsigned int __i = 2; __i < __max_iter; ++__i)
00163         {
00164           __term *= (2 * __i - 1) * __kk / (2 * __i);
00165           if (__term < std::numeric_limits<_Tp>::epsilon())
00166             break;
00167           __sum += __term;
00168         }
00169 
00170       return __numeric_constants<_Tp>::__pi_2() * __sum;
00171     }
00172 
00173 
00174     /**
00175      *   @brief  Return the complete elliptic integral of the first kind
00176      *           @f$ K(k) @f$ using the Carlson formulation.
00177      * 
00178      *   The complete elliptic integral of the first kind is defined as
00179      *   @f[
00180      *     K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
00181      *                                           {\sqrt{1 - k^2 sin^2\theta}}
00182      *   @f]
00183      *   where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the
00184      *   first kind.
00185      * 
00186      *   @param  __k  The argument of the complete elliptic function.
00187      *   @return  The complete elliptic function of the first kind.
00188      */
00189     template<typename _Tp>
00190     _Tp
00191     __comp_ellint_1(const _Tp __k)
00192     {
00193 
00194       if (__isnan(__k))
00195         return std::numeric_limits<_Tp>::quiet_NaN();
00196       else if (std::abs(__k) >= _Tp(1))
00197         return std::numeric_limits<_Tp>::quiet_NaN();
00198       else
00199         return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
00200     }
00201 
00202 
00203     /**
00204      *   @brief  Return the incomplete elliptic integral of the first kind
00205      *           @f$ F(k,\phi) @f$ using the Carlson formulation.
00206      * 
00207      *   The incomplete elliptic integral of the first kind is defined as
00208      *   @f[
00209      *     F(k,\phi) = \int_0^{\phi}\frac{d\theta}
00210      *                                   {\sqrt{1 - k^2 sin^2\theta}}
00211      *   @f]
00212      * 
00213      *   @param  __k  The argument of the elliptic function.
00214      *   @param  __phi  The integral limit argument of the elliptic function.
00215      *   @return  The elliptic function of the first kind.
00216      */
00217     template<typename _Tp>
00218     _Tp
00219     __ellint_1(const _Tp __k, const _Tp __phi)
00220     {
00221 
00222       if (__isnan(__k) || __isnan(__phi))
00223         return std::numeric_limits<_Tp>::quiet_NaN();
00224       else if (std::abs(__k) > _Tp(1))
00225         std::__throw_domain_error(__N("Bad argument in __ellint_1."));
00226       else
00227         {
00228           //  Reduce phi to -pi/2 < phi < +pi/2.
00229           const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
00230                                    + _Tp(0.5L));
00231           const _Tp __phi_red = __phi
00232                               - __n * __numeric_constants<_Tp>::__pi();
00233 
00234           const _Tp __s = std::sin(__phi_red);
00235           const _Tp __c = std::cos(__phi_red);
00236 
00237           const _Tp __F = __s
00238                         * __ellint_rf(__c * __c,
00239                                 _Tp(1) - __k * __k * __s * __s, _Tp(1));
00240 
00241           if (__n == 0)
00242             return __F;
00243           else
00244             return __F + _Tp(2) * __n * __comp_ellint_1(__k);
00245         }
00246     }
00247 
00248 
00249     /**
00250      *   @brief Return the complete elliptic integral of the second kind
00251      *          @f$ E(k) @f$ by series expansion.
00252      * 
00253      *   The complete elliptic integral of the second kind is defined as
00254      *   @f[
00255      *     E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
00256      *   @f]
00257      * 
00258      *   This routine is not bad as long as |k| is somewhat smaller than 1
00259      *   but is not is good as the Carlson elliptic integral formulation.
00260      * 
00261      *   @param  __k  The argument of the complete elliptic function.
00262      *   @return  The complete elliptic function of the second kind.
00263      */
00264     template<typename _Tp>
00265     _Tp
00266     __comp_ellint_2_series(const _Tp __k)
00267     {
00268 
00269       const _Tp __kk = __k * __k;
00270 
00271       _Tp __term = __kk;
00272       _Tp __sum = __term;
00273 
00274       const unsigned int __max_iter = 1000;
00275       for (unsigned int __i = 2; __i < __max_iter; ++__i)
00276         {
00277           const _Tp __i2m = 2 * __i - 1;
00278           const _Tp __i2 = 2 * __i;
00279           __term *= __i2m * __i2m * __kk / (__i2 * __i2);
00280           if (__term < std::numeric_limits<_Tp>::epsilon())
00281             break;
00282           __sum += __term / __i2m;
00283         }
00284 
00285       return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum);
00286     }
00287 
00288 
00289     /**
00290      *   @brief  Return the Carlson elliptic function of the second kind
00291      *           @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where
00292      *           @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function
00293      *           of the third kind.
00294      * 
00295      *   The Carlson elliptic function of the second kind is defined by:
00296      *   @f[
00297      *       R_D(x,y,z) = \frac{3}{2} \int_0^\infty
00298      *                 \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}}
00299      *   @f]
00300      *
00301      *   Based on Carlson's algorithms:
00302      *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
00303      *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
00304      *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
00305      *      by Press, Teukolsky, Vetterling, Flannery (1992)
00306      *
00307      *   @param  __x  The first of two symmetric arguments.
00308      *   @param  __y  The second of two symmetric arguments.
00309      *   @param  __z  The third argument.
00310      *   @return  The Carlson elliptic function of the second kind.
00311      */
00312     template<typename _Tp>
00313     _Tp
00314     __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
00315     {
00316       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00317       const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
00318       const _Tp __min = std::numeric_limits<_Tp>::min();
00319       const _Tp __max = std::numeric_limits<_Tp>::max();
00320       const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3));
00321       const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) / _Tp(3));
00322 
00323       if (__x < _Tp(0) || __y < _Tp(0))
00324         std::__throw_domain_error(__N("Argument less than zero "
00325                                       "in __ellint_rd."));
00326       else if (__x + __y < __lolim || __z < __lolim)
00327         std::__throw_domain_error(__N("Argument too small "
00328                                       "in __ellint_rd."));
00329       else
00330         {
00331           const _Tp __c0 = _Tp(1) / _Tp(4);
00332           const _Tp __c1 = _Tp(3) / _Tp(14);
00333           const _Tp __c2 = _Tp(1) / _Tp(6);
00334           const _Tp __c3 = _Tp(9) / _Tp(22);
00335           const _Tp __c4 = _Tp(3) / _Tp(26);
00336 
00337           _Tp __xn = __x;
00338           _Tp __yn = __y;
00339           _Tp __zn = __z;
00340           _Tp __sigma = _Tp(0);
00341           _Tp __power4 = _Tp(1);
00342 
00343           _Tp __mu;
00344           _Tp __xndev, __yndev, __zndev;
00345 
00346           const unsigned int __max_iter = 100;
00347           for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00348             {
00349               __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5);
00350               __xndev = (__mu - __xn) / __mu;
00351               __yndev = (__mu - __yn) / __mu;
00352               __zndev = (__mu - __zn) / __mu;
00353               _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
00354               __epsilon = std::max(__epsilon, std::abs(__zndev));
00355               if (__epsilon < __errtol)
00356                 break;
00357               _Tp __xnroot = std::sqrt(__xn);
00358               _Tp __ynroot = std::sqrt(__yn);
00359               _Tp __znroot = std::sqrt(__zn);
00360               _Tp __lambda = __xnroot * (__ynroot + __znroot)
00361                            + __ynroot * __znroot;
00362               __sigma += __power4 / (__znroot * (__zn + __lambda));
00363               __power4 *= __c0;
00364               __xn = __c0 * (__xn + __lambda);
00365               __yn = __c0 * (__yn + __lambda);
00366               __zn = __c0 * (__zn + __lambda);
00367             }
00368 
00369       // Note: __ea is an SPU badname.
00370           _Tp __eaa = __xndev * __yndev;
00371           _Tp __eb = __zndev * __zndev;
00372           _Tp __ec = __eaa - __eb;
00373           _Tp __ed = __eaa - _Tp(6) * __eb;
00374           _Tp __ef = __ed + __ec + __ec;
00375           _Tp __s1 = __ed * (-__c1 + __c3 * __ed
00376                                    / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef
00377                                    / _Tp(2));
00378           _Tp __s2 = __zndev
00379                    * (__c2 * __ef
00380                     + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa));
00381 
00382           return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
00383                                         / (__mu * std::sqrt(__mu));
00384         }
00385     }
00386 
00387 
00388     /**
00389      *   @brief  Return the complete elliptic integral of the second kind
00390      *           @f$ E(k) @f$ using the Carlson formulation.
00391      * 
00392      *   The complete elliptic integral of the second kind is defined as
00393      *   @f[
00394      *     E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
00395      *   @f]
00396      * 
00397      *   @param  __k  The argument of the complete elliptic function.
00398      *   @return  The complete elliptic function of the second kind.
00399      */
00400     template<typename _Tp>
00401     _Tp
00402     __comp_ellint_2(const _Tp __k)
00403     {
00404 
00405       if (__isnan(__k))
00406         return std::numeric_limits<_Tp>::quiet_NaN();
00407       else if (std::abs(__k) == 1)
00408         return _Tp(1);
00409       else if (std::abs(__k) > _Tp(1))
00410         std::__throw_domain_error(__N("Bad argument in __comp_ellint_2."));
00411       else
00412         {
00413           const _Tp __kk = __k * __k;
00414 
00415           return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
00416                - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3);
00417         }
00418     }
00419 
00420 
00421     /**
00422      *   @brief  Return the incomplete elliptic integral of the second kind
00423      *           @f$ E(k,\phi) @f$ using the Carlson formulation.
00424      * 
00425      *   The incomplete elliptic integral of the second kind is defined as
00426      *   @f[
00427      *     E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta}
00428      *   @f]
00429      * 
00430      *   @param  __k  The argument of the elliptic function.
00431      *   @param  __phi  The integral limit argument of the elliptic function.
00432      *   @return  The elliptic function of the second kind.
00433      */
00434     template<typename _Tp>
00435     _Tp
00436     __ellint_2(const _Tp __k, const _Tp __phi)
00437     {
00438 
00439       if (__isnan(__k) || __isnan(__phi))
00440         return std::numeric_limits<_Tp>::quiet_NaN();
00441       else if (std::abs(__k) > _Tp(1))
00442         std::__throw_domain_error(__N("Bad argument in __ellint_2."));
00443       else
00444         {
00445           //  Reduce phi to -pi/2 < phi < +pi/2.
00446           const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
00447                                    + _Tp(0.5L));
00448           const _Tp __phi_red = __phi
00449                               - __n * __numeric_constants<_Tp>::__pi();
00450 
00451           const _Tp __kk = __k * __k;
00452           const _Tp __s = std::sin(__phi_red);
00453           const _Tp __ss = __s * __s;
00454           const _Tp __sss = __ss * __s;
00455           const _Tp __c = std::cos(__phi_red);
00456           const _Tp __cc = __c * __c;
00457 
00458           const _Tp __E = __s
00459                         * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
00460                         - __kk * __sss
00461                         * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1))
00462                         / _Tp(3);
00463 
00464           if (__n == 0)
00465             return __E;
00466           else
00467             return __E + _Tp(2) * __n * __comp_ellint_2(__k);
00468         }
00469     }
00470 
00471 
00472     /**
00473      *   @brief  Return the Carlson elliptic function
00474      *           @f$ R_C(x,y) = R_F(x,y,y) @f$ where @f$ R_F(x,y,z) @f$
00475      *           is the Carlson elliptic function of the first kind.
00476      * 
00477      *   The Carlson elliptic function is defined by:
00478      *   @f[
00479      *       R_C(x,y) = \frac{1}{2} \int_0^\infty
00480      *                 \frac{dt}{(t + x)^{1/2}(t + y)}
00481      *   @f]
00482      *
00483      *   Based on Carlson's algorithms:
00484      *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
00485      *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
00486      *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
00487      *      by Press, Teukolsky, Vetterling, Flannery (1992)
00488      *
00489      *   @param  __x  The first argument.
00490      *   @param  __y  The second argument.
00491      *   @return  The Carlson elliptic function.
00492      */
00493     template<typename _Tp>
00494     _Tp
00495     __ellint_rc(const _Tp __x, const _Tp __y)
00496     {
00497       const _Tp __min = std::numeric_limits<_Tp>::min();
00498       const _Tp __max = std::numeric_limits<_Tp>::max();
00499       const _Tp __lolim = _Tp(5) * __min;
00500       const _Tp __uplim = __max / _Tp(5);
00501 
00502       if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
00503         std::__throw_domain_error(__N("Argument less than zero "
00504                                       "in __ellint_rc."));
00505       else
00506         {
00507           const _Tp __c0 = _Tp(1) / _Tp(4);
00508           const _Tp __c1 = _Tp(1) / _Tp(7);
00509           const _Tp __c2 = _Tp(9) / _Tp(22);
00510           const _Tp __c3 = _Tp(3) / _Tp(10);
00511           const _Tp __c4 = _Tp(3) / _Tp(8);
00512 
00513           _Tp __xn = __x;
00514           _Tp __yn = __y;
00515 
00516           const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00517           const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
00518           _Tp __mu;
00519           _Tp __sn;
00520 
00521           const unsigned int __max_iter = 100;
00522           for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00523             {
00524               __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
00525               __sn = (__yn + __mu) / __mu - _Tp(2);
00526               if (std::abs(__sn) < __errtol)
00527                 break;
00528               const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn)
00529                              + __yn;
00530               __xn = __c0 * (__xn + __lambda);
00531               __yn = __c0 * (__yn + __lambda);
00532             }
00533 
00534           _Tp __s = __sn * __sn
00535                   * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
00536 
00537           return (_Tp(1) + __s) / std::sqrt(__mu);
00538         }
00539     }
00540 
00541 
00542     /**
00543      *   @brief  Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$
00544      *           of the third kind.
00545      * 
00546      *   The Carlson elliptic function of the third kind is defined by:
00547      *   @f[
00548      *       R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty
00549      *       \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)}
00550      *   @f]
00551      *
00552      *   Based on Carlson's algorithms:
00553      *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
00554      *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
00555      *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
00556      *      by Press, Teukolsky, Vetterling, Flannery (1992)
00557      *
00558      *   @param  __x  The first of three symmetric arguments.
00559      *   @param  __y  The second of three symmetric arguments.
00560      *   @param  __z  The third of three symmetric arguments.
00561      *   @param  __p  The fourth argument.
00562      *   @return  The Carlson elliptic function of the fourth kind.
00563      */
00564     template<typename _Tp>
00565     _Tp
00566     __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p)
00567     {
00568       const _Tp __min = std::numeric_limits<_Tp>::min();
00569       const _Tp __max = std::numeric_limits<_Tp>::max();
00570       const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3));
00571       const _Tp __uplim = _Tp(0.3L)
00572                         * std::pow(_Tp(0.2L) * __max, _Tp(1)/_Tp(3));
00573 
00574       if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
00575         std::__throw_domain_error(__N("Argument less than zero "
00576                                       "in __ellint_rj."));
00577       else if (__x + __y < __lolim || __x + __z < __lolim
00578             || __y + __z < __lolim || __p < __lolim)
00579         std::__throw_domain_error(__N("Argument too small "
00580                                       "in __ellint_rj"));
00581       else
00582         {
00583           const _Tp __c0 = _Tp(1) / _Tp(4);
00584           const _Tp __c1 = _Tp(3) / _Tp(14);
00585           const _Tp __c2 = _Tp(1) / _Tp(3);
00586           const _Tp __c3 = _Tp(3) / _Tp(22);
00587           const _Tp __c4 = _Tp(3) / _Tp(26);
00588 
00589           _Tp __xn = __x;
00590           _Tp __yn = __y;
00591           _Tp __zn = __z;
00592           _Tp __pn = __p;
00593           _Tp __sigma = _Tp(0);
00594           _Tp __power4 = _Tp(1);
00595 
00596           const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00597           const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
00598 
00599           _Tp __lambda, __mu;
00600           _Tp __xndev, __yndev, __zndev, __pndev;
00601 
00602           const unsigned int __max_iter = 100;
00603           for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
00604             {
00605               __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5);
00606               __xndev = (__mu - __xn) / __mu;
00607               __yndev = (__mu - __yn) / __mu;
00608               __zndev = (__mu - __zn) / __mu;
00609               __pndev = (__mu - __pn) / __mu;
00610               _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
00611               __epsilon = std::max(__epsilon, std::abs(__zndev));
00612               __epsilon = std::max(__epsilon, std::abs(__pndev));
00613               if (__epsilon < __errtol)
00614                 break;
00615               const _Tp __xnroot = std::sqrt(__xn);
00616               const _Tp __ynroot = std::sqrt(__yn);
00617               const _Tp __znroot = std::sqrt(__zn);
00618               const _Tp __lambda = __xnroot * (__ynroot + __znroot)
00619                                  + __ynroot * __znroot;
00620               const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot)
00621                                 + __xnroot * __ynroot * __znroot;
00622               const _Tp __alpha2 = __alpha1 * __alpha1;
00623               const _Tp __beta = __pn * (__pn + __lambda)
00624                                       * (__pn + __lambda);
00625               __sigma += __power4 * __ellint_rc(__alpha2, __beta);
00626               __power4 *= __c0;
00627               __xn = __c0 * (__xn + __lambda);
00628               __yn = __c0 * (__yn + __lambda);
00629               __zn = __c0 * (__zn + __lambda);
00630               __pn = __c0 * (__pn + __lambda);
00631             }
00632 
00633       // Note: __ea is an SPU badname.
00634           _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev;
00635           _Tp __eb = __xndev * __yndev * __zndev;
00636           _Tp __ec = __pndev * __pndev;
00637           _Tp __e2 = __eaa - _Tp(3) * __ec;
00638           _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec);
00639           _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4)
00640                             - _Tp(3) * __c4 * __e3 / _Tp(2));
00641           _Tp __s2 = __eb * (__c2 / _Tp(2)
00642                    + __pndev * (-__c3 - __c3 + __pndev * __c4));
00643           _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3)
00644                    - __c2 * __pndev * __ec;
00645 
00646           return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
00647                                              / (__mu * std::sqrt(__mu));
00648         }
00649     }
00650 
00651 
00652     /**
00653      *   @brief Return the complete elliptic integral of the third kind
00654      *          @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the
00655      *          Carlson formulation.
00656      * 
00657      *   The complete elliptic integral of the third kind is defined as
00658      *   @f[
00659      *     \Pi(k,\nu) = \int_0^{\pi/2}
00660      *                   \frac{d\theta}
00661      *                 {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}}
00662      *   @f]
00663      * 
00664      *   @param  __k  The argument of the elliptic function.
00665      *   @param  __nu  The second argument of the elliptic function.
00666      *   @return  The complete elliptic function of the third kind.
00667      */
00668     template<typename _Tp>
00669     _Tp
00670     __comp_ellint_3(const _Tp __k, const _Tp __nu)
00671     {
00672 
00673       if (__isnan(__k) || __isnan(__nu))
00674         return std::numeric_limits<_Tp>::quiet_NaN();
00675       else if (__nu == _Tp(1))
00676         return std::numeric_limits<_Tp>::infinity();
00677       else if (std::abs(__k) > _Tp(1))
00678         std::__throw_domain_error(__N("Bad argument in __comp_ellint_3."));
00679       else
00680         {
00681           const _Tp __kk = __k * __k;
00682 
00683           return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
00684                - __nu
00685                * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu)
00686                / _Tp(3);
00687         }
00688     }
00689 
00690 
00691     /**
00692      *   @brief Return the incomplete elliptic integral of the third kind
00693      *          @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation.
00694      * 
00695      *   The incomplete elliptic integral of the third kind is defined as
00696      *   @f[
00697      *     \Pi(k,\nu,\phi) = \int_0^{\phi}
00698      *                       \frac{d\theta}
00699      *                            {(1 - \nu \sin^2\theta)
00700      *                             \sqrt{1 - k^2 \sin^2\theta}}
00701      *   @f]
00702      * 
00703      *   @param  __k  The argument of the elliptic function.
00704      *   @param  __nu  The second argument of the elliptic function.
00705      *   @param  __phi  The integral limit argument of the elliptic function.
00706      *   @return  The elliptic function of the third kind.
00707      */
00708     template<typename _Tp>
00709     _Tp
00710     __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi)
00711     {
00712 
00713       if (__isnan(__k) || __isnan(__nu) || __isnan(__phi))
00714         return std::numeric_limits<_Tp>::quiet_NaN();
00715       else if (std::abs(__k) > _Tp(1))
00716         std::__throw_domain_error(__N("Bad argument in __ellint_3."));
00717       else
00718         {
00719           //  Reduce phi to -pi/2 < phi < +pi/2.
00720           const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
00721                                    + _Tp(0.5L));
00722           const _Tp __phi_red = __phi
00723                               - __n * __numeric_constants<_Tp>::__pi();
00724 
00725           const _Tp __kk = __k * __k;
00726           const _Tp __s = std::sin(__phi_red);
00727           const _Tp __ss = __s * __s;
00728           const _Tp __sss = __ss * __s;
00729           const _Tp __c = std::cos(__phi_red);
00730           const _Tp __cc = __c * __c;
00731 
00732           const _Tp __Pi = __s
00733                          * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
00734                          - __nu * __sss
00735                          * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1),
00736                                        _Tp(1) + __nu * __ss) / _Tp(3);
00737 
00738           if (__n == 0)
00739             return __Pi;
00740           else
00741             return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu);
00742         }
00743     }
00744 
00745   } // namespace std::tr1::__detail
00746 }
00747 }
00748 
00749 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC
00750 

Generated on Thu Jul 23 21:16:04 2009 for libstdc++ by  doxygen 1.5.8