CLHEP/GenericFunctions/EllipticIntegral.icc

0001 // -*- C++ -*-
0002 // $Id:
0003 #include "gsl/gsl_sf_ellint.h"
0004 #include <cmath>
0005 #include <signal.h>
0006 #include <assert.h>
0007
0008
0009 namespace Genfun {
0010 namespace EllipticIntegral {
0011
0012 //-----------------------------------------------------------------------------//
0013 //                           FIRST KIND                                        //
0014 //-----------------------------------------------------------------------------//
0015
0016 FUNCTION_OBJECT_IMP(FirstKind)
0017
0018 inline
0019 FirstKind::FirstKind():
0020   _k("K", 1.0,0.0,1.0)
0021 {
0022 }
0023
0024 inline
0025 FirstKind::~FirstKind() {
0026 }
0027
0028 inline
0029 FirstKind::FirstKind(const FirstKind & right):
0030   _k(right._k)
0031 {
0032 }
0033
0034
0035 inline
0036 Parameter & FirstKind::k() {
0037   return _k;
0038 }
0039
0040 inline
0041 const Parameter & FirstKind::k() const {
0042   return _k;
0043 }
0044
0045
0046 inline
0047 double FirstKind::operator() (double x) const {
0048   gsl_sf_result result;
0049   int status = gsl_sf_ellint_F_e(x,_k.getValue(), GSL_PREC_DOUBLE,  &result);
0050   if (status!=0) {
0051     std::cerr << "Warning, GSL function gsl_sf_ellint_F_impl"
0052           << " return code" << status << std::endl;
0053     raise(SIGFPE);
0054   }
0055   return result.val;
0056 }
0057 //-----------------------------------------------------------------------------//
0058 //                           SECOND KIND                                       //
0059 //-----------------------------------------------------------------------------//
0060
0061 FUNCTION_OBJECT_IMP(SecondKind)
0062
0063 inline
0064 SecondKind::SecondKind():
0065   _k("K", 1.0,0.0,1.0)
0066 {
0067 }
0068
0069 inline
0070 SecondKind::~SecondKind() {
0071 }
0072
0073 inline
0074 SecondKind::SecondKind(const SecondKind & right):
0075   _k(right._k)
0076 {
0077 }
0078
0079
0080 inline
0081 Parameter & SecondKind::k() {
0082   return _k;
0083 }
0084
0085 inline
0086 const Parameter & SecondKind::k() const {
0087   return _k;
0088 }
0089
0090
0091 inline
0092 double SecondKind::operator() (double x) const {
0093   gsl_sf_result result;
0094   int status = gsl_sf_ellint_E_e(x,_k.getValue(), GSL_PREC_DOUBLE,  &result);
0095   if (status!=0) {
0096     std::cerr << "Warning, GSL function gsl_sf_ellint_E_impl"
0097           << " return code" << status << std::endl;
0098     raise(SIGFPE);
0099   }
0100   return result.val;
0101 }
0102 //-----------------------------------------------------------------------------//
0103 //                           THIRD KIND                                        //
0104 //-----------------------------------------------------------------------------//
0105 FUNCTION_OBJECT_IMP(ThirdKind)
0106
0107 inline
0108 ThirdKind::ThirdKind():
0109   _k("K", 1.0,0.0, 1.0),
0110   _n("N", 1.0,0.0,10.0)
0111 {
0112 }
0113
0114 inline
0115 ThirdKind::~ThirdKind() {
0116 }
0117
0118 inline
0119 ThirdKind::ThirdKind(const ThirdKind & right):
0120   _k(right._k),
0121   _n(right._n)
0122 {
0123 }
0124
0125
0126 inline
0127 Parameter & ThirdKind::k() {
0128   return _k;
0129 }
0130
0131 inline
0132 const Parameter & ThirdKind::k() const {
0133   return _k;
0134 }
0135
0136
0137 inline
0138 Parameter & ThirdKind::n() {
0139   return _n;
0140 }
0141
0142 inline
0143 const Parameter & ThirdKind::n() const {
0144   return _n;
0145 }
0146
0147
0148 inline
0149 double ThirdKind::operator() (double x) const {
0150   gsl_sf_result result;
0151   int status = gsl_sf_ellint_P_e(x,_k.getValue(),_n.getValue(), GSL_PREC_DOUBLE,  &result);
0152   if (status!=0) {
0153     std::cerr << "Warning, GSL function gsl_ellint_P_impl"
0154           << " return code" << status << std::endl;
0155     raise(SIGFPE);
0156   }
0157   return result.val;
0158 }
0159 } // end namespace EllipticIntegral
0160 } // end namespace Genfun