Bessel Function Second Kind
11-01-2017, 12:20 AM (This post was last modified: 11-10-2017 01:18 AM by toml_12953.)
Post: #1
 toml_12953 Senior Member Posts: 1,221 Joined: Dec 2013
Bessel Function Second Kind
You can run Tbessy for an example where N = 2 and X = 0.75 or you can call BESSY(N,X) directly with your own values.

Code:
//************************************************************************ //*                                                                      * //*    Program to calculate the second kind Bessel function of integer   * //*    order N, for any REAL X, using the function BESSY(N,X).           * //*                                                                      * //* -------------------------------------------------------------------- * //*                                                                      * //*    SAMPLE RUN:                                                       * //*                                                                      * //*    (Calculate Bessel function for N:=2, X:=0.75).                    * //*                                                                      * //*    Second kind Bessel function of order  2 for X =  0.7500:          * //*                                                                      * //*         Y := -2.62974604                                             * //*                                                                      * //* -------------------------------------------------------------------- * //*   Reference: From Numath Library By Tuan Dang Trong in Fortran 77.   * //*                                                                      * //*                               C++ Release 1.0 By J-P Moreau, Paris.  * //*                                        (www.jpmoreau.fr)             * //*                                  Translated to HPPL by Tom Lake      * //***********************************************************************/ BESSJ0 (X) BEGIN //*********************************************************************** //      This subroutine calculates the First Kind Bessel Function of //      order 0, for any real number X. The polynomial approximation by //      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1. //      REFERENCES: //      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965. //      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES, //      VOL.5, 1962. //************************************************************************/     LOCAL P1:=1.0, P2:=-0.1098628627E-2, P3:=0.2734510407E-4;     LOCAL P4:=-0.2073370639E-5, P5:= 0.2093887211E-6;     LOCAL Q1:=-0.1562499995E-1, Q2:= 0.1430488765E-3, Q3:=-0.6911147651E-5;     LOCAL Q4:= 0.7621095161E-6, Q5:=-0.9349451520E-7;     LOCAL R1:= 57568490574.0, R2:=-13362590354.0, R3:=651619640.7;     LOCAL R4:=-11214424.18, R5:= 77392.33017, R6:=-184.9052456;     LOCAL S1:= 57568490411.0, S2:=1029532985.0, S3:=9494680.718;     LOCAL S4:= 59272.64853, S5:=267.8532712, S6:=1.0;     LOCAL AX,FR,FS,Z,FP1,FQ,XX,Y;     LOCAL TMP;     IF (X==0.0) THEN RETURN 1.0; END;     AX := ABS(X);     IF (AX < 8.0) THEN         Y := X*X;         FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));         FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));         TMP := FR/FS;          ELSE          Z := 8./AX;         Y := Z*Z;         XX := AX-0.785398164;         FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));         FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));         TMP := SQRT(0.636619772/AX)*(FP1*COS(XX)-Z*FQ*SIN(XX));     END;     RETURN TMP; END; Signu (X,Y) BEGIN     IF (Y<0.0) THEN          RETURN (-ABS(X));     ELSE          RETURN (ABS(X));     END; END; BESSJ1 (X) BEGIN   //********************************************************************** //      This subroutine calculates the First Kind Bessel Function of //     order 1, for any real number X. The polynomial approximation by //      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1. //      REFERENCES: //      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965. //      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES, //      VOL.5, 1962. //***********************************************************************/     LOCAL P1:=1.0, P2:=0.183105E-2, P3:=-0.3516396496E-4, P4:=0.2457520174E-5;     LOCAL P5:=-0.240337019E-6,  P6:=0.636619772;     LOCAL Q1:= 0.04687499995, Q2:=-0.2002690873E-3, Q3:=0.8449199096E-5;     LOCAL Q4:=-0.88228987E-6, Q5:= 0.105787412E-6;     LOCAL R1:= 72362614232.0, R2:=-7895059235.0, R3:=242396853.1;     LOCAL R4:=-2972611.439,   R5:=15704.48260,  R6:=-30.16036606;     LOCAL S1:=144725228442.0, S2:=2300535178.0, S3:=18583304.74;     LOCAL S4:=99447.43394,    S5:=376.9991397,  S6:=1.0;     LOCAL AX,FR,FS,Y,Z,FP1,FQ,XX;     LOCAL TMP;     AX := ABS(X);     IF (AX < 8.0) THEN         Y := X*X;         FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));         FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));         TMP := X*(FR/FS);           ELSE         Z := 8.0/AX;         Y := Z*Z;         XX := AX-2.35619491;         FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));         FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));         TMP := SQRT(P6/AX)*(COS(XX)*FP1-Z*SIN(XX)*FQ)*Signu(S6,X);     END;     RETURN TMP; END; BESSY0(X) BEGIN //* -------------------------------------------------------------------- //      This subroutine calculates the Second Kind Bessel Function of //      order 0, for any real number X. The polynomial approximation by //      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1. //      REFERENCES: //      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965. //      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES, //      VOL.5, 1962. //  --------------------------------------------------------------------- */     LOCAL P1:= 1.0, P2:=-0.1098628627E-2, P3:=0.2734510407E-4;     LOCAL P4:=-0.2073370639E-5, P5:= 0.2093887211E-6;     LOCAL Q1:=-0.1562499995E-1, Q2:= 0.1430488765E-3, Q3:=-0.6911147651E-5;     LOCAL Q4:= 0.7621095161E-6, Q5:=-0.9349451520E-7;     LOCAL R1:=-2957821389.0, R2:=7062834065.0, R3:=-512359803.6;     LOCAL R4:= 10879881.29,  R5:=-86327.92757, R6:=228.4622733;     LOCAL S1:= 40076544269.0, S2:=745249964.8, S3:=7189466.438;     LOCAL S4:= 47447.26470,   S5:=226.1030244, S6:=1.0;     LOCAL FS,FR,Z,FP1,FQ,XX,Y;     IF (X == 0.0) THEN RETURN -1e30; END;     IF (X < 8.0) THEN         Y := X*X;         FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));         FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));         RETURN (FR/FS+0.636619772*BESSJ0(X)*LN(X));     ELSE          Z := 8.0/X;         Y := Z*Z;         XX := X-0.785398164;         FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));         FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));         RETURN (SQRT(0.636619772/X)*(FP1*SIN(XX)+Z*FQ*COS(XX)));     END; END; BESSY1(X) BEGIN //* --------------------------------------------------------------------- //      This subroutine calculates the Second Kind Bessel Function of //      order 1, for any real number X. The polynomial approximation by //      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1. //      REFERENCES: //      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965. //      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES, //      VOL.5, 1962. //  ---------------------------------------------------------------------- */     LOCAL P1:= 1.0, P2:=0.183105E-2, P3:=-0.3516396496E-4;     LOCAL P4:= 0.2457520174E-5, P5:=-0.240337019E-6;     LOCAL Q1:= 0.04687499995, Q2:=-0.2002690873E-3, Q3:=0.8449199096E-5;     LOCAL Q4:=-0.88228987E-6, Q5:= 0.105787412E-6;     LOCAL R1:=-0.4900604943E13, R2:= 0.1275274390E13, R3:=-0.5153438139E11;     LOCAL R4:= 0.7349264551E9,  R5:=-0.4237922726E7,  R6:= 0.8511937935E4;     LOCAL S1:= 0.2499580570E14, S2:= 0.4244419664E12, S3:= 0.3733650367E10;     LOCAL S4:= 0.2245904002E8,  S5:= 0.1020426050E6,  S6:= 0.3549632885E3, S7:=1.0;     LOCAL FR,FS,Z,FP1,FQ,XX, Y;     IF (X == 0.0) THEN RETURN -1e30; END;     IF (X < 8.0) THEN         Y := X*X;         FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));         FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*(S6+Y*S7)))));         RETURN (X*(FR/FS)+0.636619772*(BESSJ1(X)*LN(X)-1.0/X));           ELSE          Z := 8./X;         Y := Z*Z;         XX := X-2.356194491;         FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));         FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));         RETURN (SQRT(0.636619772/X)*(SIN(XX)*FP1+Z*COS(XX)*FQ));     END; END; EXPORT BESSY(N,X)   BEGIN //* ----------------------------------------------------------------- //      This subroutine calculates the second kind Bessel Function of //      integer order N, for any real X. We use here the classical //      recursive formula.  //  ------------------------------------------------------------------ */     LOCAL TOX,BY,BYM,BYP,J;     IF (N == 0) THEN RETURN BESSY0(X); END;     IF (N == 1) THEN RETURN BESSY1(X); END;     IF (X == 0.0) THEN RETURN -1e30;  END;     TOX := 2.0/X;     BY  := BESSY1(X);     BYM := BESSY0(X);     FOR J:=1 TO N-1 DO         BYP := J*TOX*BY-BYM;         BYM := BY;         BY  := BYP;     END;     RETURN BY;   END; // -------------------------------------------------------------------------- BESSYP(N,X) BEGIN     IF (N == 0) THEN         RETURN (-BESSY(1,X));     ELSE          IF (X == 0.0) THEN             RETURN 1e-30;         ELSE             RETURN (BESSY(N-1,X)-(1.0*N/X)*BESSY(N,X));         END;     END; END; EXPORT Tbessy() BEGIN   LOCAL X,Y,N;   N:=2;   X:=0.75;   Y := BESSY(N,X);   PRINT("\n Second Kind Bessel Function of order " + N + ", for X= " + X + "\n\n");   PRINT("      Y = "+CAS(format(Y,"s9"))+"\n\n"); END; // End of file Tbessy.hpprgm

Tom L
I think therefore I am-Descartes
I think therefore you are-Gorgias
You're not here to think-Army Sergeant
11-01-2017, 09:09 AM
Post: #2
 salvomic Senior Member Posts: 1,366 Joined: Jan 2015
RE: Bessel Function Second Kind
Thank you Tom!

I tested it with some values in the literature and the function is very precise.

Salvo Micciché

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)