📄 tbessk_cpp.txt
字号:
/**********************************************************
* This program tests the subroutine BESSK to calculate the*
* modified Bessel function of the third kind of order N *
* for any real positive argument X. *
* ------------------------------------------------------- *
* SAMPLE RUN: *
* *
* X = 1.209700 *
* *
* N = 0 *
* Y = 0.314324 *
* *
* N = 1 *
* Y = 0.428051 *
* *
* N = 2 *
* Y = 1.022022 *
* *
* N = 3 *
* Y = 3.807473 *
* *
* N = 4 *
* Y = 19.906737 *
* *
* ------------------------------------------------------- *
* Reference: From Numath Library By Tuan Dang Trong *
* in Fortran 77. *
* *
* C++ Version By J-P Moreau, Paris. *
**********************************************************/
#include <stdio.h>
#include <math.h>
#define BIG 1e30
double X, Y;
int i, N;
double BESSK0(double X);
double BESSK1(double X);
double BESSI0(double X);
double BESSI1(double X);
double BESSK(int N, double X) {
/* -----------------------------------------------------------------------
! CE SOUS-PROGRAMME CALCULE LA FONCTION BESSEL MODIFIFIEE 3E ESPECE
! D'ORDRE N ENTIER POUR TOUT X REEL POSITIF > 0. ON UTILISE ICI LA
! FORMULE DE RECURRENCE CLASSIQUE EN PARTANT DE BESSK0 ET BESSK1.
!
! THIS ROUTINE CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD
! KIND OF INTEGER ORDER, N FOR ANY POSITIVE REAL ARGUMENT, X. THE
! CLASSICAL RECURSION FORMULA IS USED, STARTING FROM BESSK0 AND BESSK1.
! ------------------------------------------------------------------------
! REFERENCE:
! C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
! MATHEMATICAL TABLES, VOL.5, 1962.
! ----------------------------------------------------------------------*/
double TOX,BK,BKM,BKP;
int J;
if (N == 0) return BESSK0(X);
if (N == 1) return BESSK1(X);
if (X == 0.0) return BIG; //arbitrary big value
TOX = 2.0/X;
BK = BESSK1(X);
BKM = BESSK0(X);
for (J=1; J<N; J++) {
BKP = BKM + J*TOX*BK;
BKM = BK;
BK = BKP;
}
return BK;
}
double BESSK0(double X) {
/* ---------------------------------------------------------------------
! CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0
! POUR TOUT X REEL NON NUL.
!
! CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF
! ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X.
! --------------------------------------------------------------------*/
double Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,tmp;
P1=-0.57721566; P2= 0.42278420; P3=0.23069756; P4=0.3488590e-1;
P5= 0.262698e-2; P6=0.10750e-3; P7=0.74e-5;
Q1= 1.25331414; Q2=-0.7832358e-1; Q3=0.2189568e-1; Q4=-0.1062446e-1;
Q5= 0.587872e-2; Q6=-0.251540e-2; Q7=0.53208e-3;
if (X == 0.0) return BIG; //arbitrary big value
if (X <= 2.0) {
Y=X*X/4.0;
AX=-log(X/2.0)*BESSI0(X);
tmp = AX+(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))));
return tmp;
}
else {
Y=2.0/X;
AX=exp(-X)/sqrt(X);
tmp = AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))));
return tmp;
}
}
double BESSK1(double X) {
/* ------------------------------------------------------------------------
! CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1
! POUR TOUT X REEL POSITF NON NUL.
!
! CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF
! ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X.
! -----------------------------------------------------------------------*/
double Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7;
P1=1.0; P2=0.15443144; P3=-0.67278579; P4=-0.18156897;
P5=-0.1919402e-1; P6=-0.110404e-2; P7=-0.4686e-4;
Q1=1.25331414; Q2=0.23498619; Q3=-0.3655620e-1; Q4=0.1504268e-1;
Q5=-0.780353e-2; Q6=0.325614e-2; Q7=-0.68245e-3;
if (X == 0.0) return BIG;
if (X <= 2.0) {
Y=X*X/4.0;
AX=log(X/2.0)*BESSI1(X);
return (AX+(1.0/X)*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))));
}
else {
Y=2.0/X;
AX=exp(-X)/sqrt(X);
return (AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))));
}
}
// Bessel Function of the 1st kind of order zero
double BESSI0(double X) {
double Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX;
P1=1.0; P2=3.5156229; P3=3.0899424; P4=1.2067429;
P5=0.2659732; P6=0.360768e-1; P7=0.45813e-2;
Q1=0.39894228; Q2=0.1328592e-1; Q3=0.225319e-2; Q4=-0.157565e-2;
Q5=0.916281e-2; Q6=-0.2057706e-1; Q7=0.2635537e-1;
Q8=-0.1647633e-1; Q9=0.392377e-2;
if (fabs(X) < 3.75) {
Y=pow((X/3.75),2);
return (P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))));
}
else {
AX=fabs(X);
Y=3.75/AX;
BX=exp(AX)/sqrt(AX);
AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))));
return AX*BX;
}
}
// Bessel Function of the 1st kind of order one
double BESSI1(double X) {
double Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX;
P1=0.5; P2=0.87890594; P3=0.51498869; P4=0.15084934;
P5=0.2658733e-1; P6=0.301532e-2; P7=0.32411e-3;
Q1=0.39894228; Q2=-0.3988024e-1; Q3=-0.362018e-2;
Q4=0.163801e-2;Q5=-0.1031555e-1; Q6= 0.2282967e-1;
Q7=-0.2895312e-1; Q8=0.1787654e-1; Q9=-0.420059e-2;
if (fabs(X) < 3.75) {
Y=pow((X/3.75),2);
return (X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))));
}
else {
AX=fabs(X);
Y=3.75/AX;
BX=exp(AX)/sqrt(AX);
AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))));
return AX*BX;
}
}
void main() {
N=5; // number of integer orders
X=1.2097; // value of real argument
printf("\n X = %f\n\n", X);
for (i=0; i<N; i++) {
Y=BESSK(i,X);
printf(" N = %d\n", i);
printf(" Y = %f\n\n", Y);
}
printf("\n");
}
// End of file Tbessk.cpp
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -