File indexing completed on 2024-09-27 07:03:38
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
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
0048
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
0077
0078
0079
0080
0081
0082
0083
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;
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
0112
0113
0114
0115
0116
0117
0118
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
0148
0149
0150
0151
0152
0153
0154
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 }