estarlight/src/bessel.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 28                          $: revision of last commit
0024 // $Author:: bgrube                   $: author of last commit
0025 // $Date:: 2010-12-10 18:30:01 +0000 #$: date of last commit
0026 //
0027 // Description:
0028 //    Bessel functions taken from ROOT
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "bessel.h"
0039
0040
0041 using namespace std;
0042
0043
0044 //______________________________________________________________________________
0045 double bessel::besI0(double x)
0046 {
0047   //FROM ROOT...
0048   // Parameters of the polynomial approximation
0049   const double p1=1.0,          p2=3.5156229,    p3=3.0899424,
0050     p4=1.2067492,    p5=0.2659732,    p6=3.60768e-2,  p7=4.5813e-3;
0051
0052   const double q1= 0.39894228,  q2= 1.328592e-2, q3= 2.25319e-3,
0053     q4=-1.57565e-3,  q5= 9.16281e-3,  q6=-2.057706e-2,
0054     q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3;
0055
0056   const double k1 = 3.75;
0057   double ax = fabs(x);
0058
0059   double y=0., result=0.;
0060
0061   if (ax < k1) {
0062     double xx = x/k1;
0063     y = xx*xx;
0064     result = p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))));
0065   } else {
0066     y = k1/ax;
0067     result = (exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
0068   }
0069   return result;
0070 }
0071
0072
0073 //______________________________________________________________________________
0074 double bessel::dbesk0(double x)
0075 {
0076    // Compute the modified Bessel function K_0(x) for positive real x.  //should be k0?
0077    //
0078    //  M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
0079    //     Applied Mathematics Series vol. 55 (1964), Washington.
0080    //
0081    //--- NvE 12-mar-2000 UU-SAP Utrecht
0082
0083    // Parameters of the polynomial approximation
0084   const double p1= -0.57721566,   p2= 0.42278420,  p3=0.23069756,
0085                p4=3.488590e-2,  p5=2.62698e-3, p6=1.0750e-4,  p7=7.4e-6;
0086
0087    const double q1= 1.25331414,  q2= -7.832358e-2,  q3=2.189568e-2,
0088                 q4= -1.062446e-2, q5=5.87872e-3,  q6= -2.51540e-3,  q7= 5.3208e-4;
0089
0090    if (x <= 0) {
0091      cout << "BesselK0 *K0* Invalid argument x = " << x << endl;  //Should be k0?
0092      return 0;
0093    }
0094
0095    double y=0.,result=0.;
0096
0097    if (x <= 2) {
0098      y = x*x/4.;
0099      result = (-log(x/2.)*bessel::besI0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
0100    } else {
0101      y = 2./x;
0102      result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
0103    }
0104    return result;
0105 }
0106
0107
0108 //______________________________________________________________________________
0109 double bessel::besI1(double x)
0110 {
0111    // Compute the modified Bessel function I_1(x) for any real x.
0112    //
0113    //  M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
0114    //     Applied Mathematics Series vol. 55 (1964), Washington.
0115    //
0116    //--- NvE 12-mar-2000 UU-SAP Utrecht
0117
0118    // Parameters of the polynomial approximation
0119    const double p1=0.5,          p2=0.87890594,   p3=0.51498869,
0120                 p4=0.15084934,   p5=2.658733e-2,  p6=3.01532e-3,  p7=3.2411e-4;
0121
0122    const double q1= 0.39894228,  q2=-3.988024e-2, q3=-3.62018e-3,
0123                   q4= 1.63801e-3,  q5=-1.031555e-2, q6= 2.282967e-2,
0124                   q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3;
0125
0126    const double k1 = 3.75;
0127    double ax = fabs(x);
0128
0129    double y=0., result=0.;
0130
0131    if (ax < k1) {
0132      double xx = x/k1;
0133      y = xx*xx;
0134      result = x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
0135    } else {
0136      y = k1/ax;
0137      result = (exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
0138      if (x < 0) result = -result;
0139    }
0140    return result;
0141 }
0142
0143
0144 //______________________________________________________________________________
0145 double bessel::dbesk1(double x)
0146 {
0147    // Compute the modified Bessel function K_1(x) for positive real x.
0148    //
0149    //  M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
0150    //     Applied Mathematics Series vol. 55 (1964), Washington.
0151    //
0152    //--- NvE 12-mar-2000 UU-SAP Utrecht
0153
0154    // Parameters of the polynomial approximation
0155    const double p1= 1.,          p2= 0.15443144,  p3=-0.67278579,
0156                 p4=-0.18156897,  p5=-1.919402e-2, p6=-1.10404e-3,  p7=-4.686e-5;
0157
0158    const double q1= 1.25331414,  q2= 0.23498619,  q3=-3.655620e-2,
0159                 q4= 1.504268e-2, q5=-7.80353e-3,  q6= 3.25614e-3,  q7=-6.8245e-4;
0160
0161    if (x <= 0) {
0162      cout << "bessel:dbesk1 *K1* Invalid argument x = " << x << endl;
0163      return 0;
0164    }
0165
0166    double y=0.,result=0.;
0167
0168    if (x <= 2) {
0169      y = x*x/4.;
0170      result = (log(x/2.)*bessel::besI1(x))+(1./x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
0171    } else {
0172      y = 2./x;
0173      result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
0174    }
0175    return result;
0176 }