CLHEP/GenericFunctions/Bessel.icc

0001 // -*- C++ -*-
0002 // $Id:
0003 #include "gsl/gsl_sf_bessel.h"
0004 #include <cmath>
0005 #include <signal.h>
0006 #include <assert.h>
0007
0008
0009 #define GF_DBL_EPSILON        2.2204460492503131e-16
0010
0011
0012 namespace Genfun {
0013 namespace IntegralOrder {
0014
0015 FUNCTION_OBJECT_IMP(Bessel)
0016
0017 inline
0018 Bessel::Bessel(Type type, unsigned int order):
0019   _type(type),_order(order)
0020 {
0021 }
0022
0023 inline
0024 Bessel::~Bessel() {
0025 }
0026
0027 inline
0028 Bessel::Bessel(const Bessel & right):
0029   _type(right._type),
0030   _order(right._order)
0031 {
0032 }
0033
0034 inline
0035 double Bessel::operator() (double x) const {
0036   gsl_sf_result result;
0037   if (_type==J) {
0038     int status = gsl_sf_bessel_Jn_e(_order, x,  &result);
0039     if (status!=0) {
0040       std::cerr << "Warning, GSL function gsl_sf_bessel_Jn_impl"
0041         << " return code" << status << std::endl;
0042       raise(SIGFPE);
0043     }
0044     return result.val;
0045   }
0046   else if (_type==Y) {
0047     int status = gsl_sf_bessel_Yn_e(_order, x,  &result);
0048     if (status!=0) {
0049       std::cerr << "Warning, GSL function gsl_sf_bessel_Yn_impl"
0050         << " return code" << status << std::endl;
0051       raise(SIGFPE);
0052     }
0053     return result.val;
0054   }
0055   else {
0056     return 0;
0057   }
0058 }
0059
0060 } // end namespace IntegralOrder
0061
0062 namespace FractionalOrder {
0063
0064 FUNCTION_OBJECT_IMP(Bessel)
0065
0066 inline
0067 Bessel::Bessel(Type type):
0068   _type(type),
0069   _order("Order", 0.0,-10,10)
0070 {
0071 }
0072
0073 inline
0074 Bessel::~Bessel() {
0075 }
0076
0077 inline
0078 Bessel::Bessel(const Bessel & right):
0079   _type(right._type),
0080   _order(right._order)
0081 {
0082 }
0083
0084
0085 inline
0086 Parameter & Bessel::order() {
0087   return _order;
0088 }
0089
0090 inline
0091 const Parameter & Bessel::order() const {
0092   return _order;
0093 }
0094
0095
0096 inline
0097 double Bessel::operator() (double x) const {
0098   gsl_sf_result result;
0099   if (_type==J) {
0100     int status = gsl_sf_bessel_Jnu_e(_order.getValue(), x,  &result);
0101     if (status!=0) {
0102       std::cerr << "Warning, GSL function gsl_sf_bessel_Jnu_impl"
0103         << " return code" << status << std::endl;
0104       raise(SIGFPE);
0105     }
0106     return result.val;
0107   }
0108   else if (_type==Y) {
0109     int status = gsl_sf_bessel_Ynu_e(_order.getValue(), x,  &result);
0110     if (status!=0) {
0111       std::cerr << "Warning, GSL function gsl_sf_bessel_Ynu_impl"
0112         << " return code" << status << std::endl;
0113       raise(SIGFPE);
0114     }
0115     return result.val;
0116   }
0117   return result.val;
0118 }
0119
0120
0121 } // end namespace FractionalOrder
0122
0123 } // end namespace Genfun