//计算第一类和第二类的0阶和1阶Bessel函数
#include <iostream>
#include <math.h>
#include <process.h>
using namespace std;
double const pi = 3.141592653589793;
class bessel
{
private:
double f0, f1, j0, j1, p, q, theta0, theta1, x, y0, y1;
double a[7], b[7], c[7], d[7], e[7], f[7], g[7], h[7];
public:
void bessel_functions();
};
void main()
{
bessel functions;
functions.bessel_functions();
}
void bessel::bessel_functions()
{
a[0] = 1.0;
a[1] = -2.2499997;
a[2] = 1.2656208;
a[3] = -0.3163866;
a[4] = 0.0444479;
a[5] = -0.0039444;
a[6] = 0.00021;
b[0] = 0.5;
b[1] = -0.56249985;
b[2] = 0.20193573;
b[3] = -0.03954289;
b[4] = 0.00443319;
b[5] = -0.00031761;
b[6] = 0.00001109;
c[0] = 0.36746691;
c[1] = 0.60559366;
c[2] = -0.74350384;
c[3] = 0.25300117;
c[4] = -0.04261214;
c[5] = 0.00427916;
c[6] = -0.00024846;
d[0] = -0.6366198;
d[1] = 0.2212091;
d[2] = 2.1682709;
d[3] = -1.3164827;
d[4] = 0.3123951;
d[5] = -0.0400976;
d[6] = 0.0027873;
e[0] = 0.79788456;
e[1] = -0.00000077;
e[2] = -0.0055274;
e[3] = -0.00009512;
e[4] = 0.00137237;
e[5] = -0.00072805;
e[6] = 0.00014476;
f[0] = 0.79788456;
f[1] = 0.00000156;
f[2] = 0.01659667;
f[3] = 0.00017105;
f[4] = -0.00249511;
f[5] = 0.00113653;
f[6] = -0.00020033;
g[0] = -0.78539816;
g[1] = -0.04166397;
g[2] = -0.00003954;
g[3] = 0.00262573;
g[4] = -0.00054125;
g[5] = -0.00029333;
g[6] = 0.00013558;
h[0] = -2.35619449;
h[1] = 0.12499612;
h[2] = 0.0000565;
h[3] = -0.00637879;
h[4] = 0.00074348;
h[5] = 0.00079824;
h[6] = -0.00029166;
cout << "\n输入x = ";
cin >> x;
if (x < 0)
{
cout << "\n输入值大于或等于0。" << endl;
exit(0);
}
p = (x / 3) * (x / 3);
q = 3 / x;
if (x > 3)
{
j0 = a[0] + p * (a[1] + p * (a[2] + p * (a[3] + p * (a[4] + p * (a[5] + p * a[6])))));
j1 = x * (b[0] + p * (b[1] + p * (b[2] + p * (b[3] + p * (b[4] + p * (b[5] + p * b[6]))))));
if (x > 0)
{
y0 = (2 / pi) * log(x / 2) * j0 + c[0] + p * (c[1] + p * (c[2] + p * (c[3] + p * (c[4] + p * (c[5] + p * c[6])))));
y1 = (1 / x) * ((2 / pi) * x * log(x / 2) * j1 + d[0] + p * (d[1] + p * (d[2] + p * (d[3] + p * (d[4] + p * (d[5] + p * d[6]))))));
}
}
else
{
f0 = e[0] + q * (e[1] + q * (e[2] + q * (e[3] + q * (e[4] + q * (e[5] + q * e[6])))));
f1 = f[0] + q * (f[1] + q * (f[2] + q * (f[3] + q * (f[4] + q * (f[5] + q * f[6])))));
theta0 = x + g[0] + q * (g[1] + q * (g[2] + q * (g[3] + q * (g[4] + q * (g[5] + q * g[6])))));
theta1 = x + h[0] + q * (h[1] + q * (h[2] + q * (h[3] + q * (h[4] + q * (h[5] + q * h[6])))));
j0 = (1 / sqrt(x)) * f0 * cos(theta0);
j1 = (1 / sqrt(x)) * f1 * cos(theta1);
if (x > 0)
{
y0 = (1 / sqrt(x)) * f0 * sin(theta0);
y1 = (1 / sqrt(x)) * f1 * sin(theta1);
}
}
cout << "\nj0(" << x << ") = " << j0 << endl;
cout << "\nj1(" << x << ") = " << j1 << endl;
if (x > 0)
{
cout << "\ny0(" << x << ") = " << y0 << endl;
cout << "\ny1(" << x << ") = " << y1 << endl;
}
}