Post Reply 
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
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
Find all posts by this user
Quote this message in a reply
11-01-2017, 09:09 AM
Post: #2
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
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




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