Bessel函数

//计算第一类和第二类的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;
 }
}