'************************************************************************
'*                                                                      *
'*    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.629746036321392                                       *
'*                                                                      *
'* -------------------------------------------------------------------- *
'*   Reference: From Numath Library By Tuan Dang Trong in Fortran 77.   *
'*                                                                      *
'*                             Basic Release 1.0 By J-P Moreau, Paris.  *
'*                                       (www.jpmoreau.fr)              *
'************************************************************************
'PROGRAM TBESSY

DEFINT I-N
DEFDBL A-H, O-Z

'   X, Y: Double
'   N: Integer

  CLS

  PRINT
  PRINT " BESSEL FUNCTIONS OF SECOND KIND"

  PRINT
  INPUT " ORDER = ", N
  PRINT
  INPUT "     X = ", X

  GOSUB 500  'CALL BESSY(N, X)

  PRINT
  PRINT " BESSY OF ORDER "; N; " = "; BESSY
  PRINT
  
END 'of main program

500  'FUNCTION BESSY (N, X)
' ----------------------------------------------------------------
'     This function calculates the second kind Bessel Function of
'     integer order N, for any real X. We use here the classical
'     recursive formula.
' ----------------------------------------------------------------
      IF N = 0 THEN
	GOSUB 800
	BESSY = BESSY0
	RETURN
      END IF
      IF N = 1 THEN
	GOSUB 900
	BESSY = BESSY1
	RETURN
      END IF
      IF X = 0# THEN
	BESSY = -1D+30
	RETURN
      END IF
      TOX = 2# / X
      GOSUB 900: BY = BESSY1
      GOSUB 800: BYM = BESSY0
      FOR J = 1 TO N - 1
	BYP = J * TOX * BY - BYM
	BYM = BY
	BY = BYP
      NEXT J
      BESSY = BY
RETURN

600   'FUNCTION BESSJ0 (X)
      
'     THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION
'     OF ORDER 0 FOR ANY REAL X. WE USE HERE THE POLYNOMIAL APPROXIMATION
'     BY SERIES OF CHEBYSHEV POLYNOMIALS 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.

      P1 = 1#: P2 = -.001098628627#: P3 = .00002734510407#
      P4 = -.000002073370639#: P5 = .0000002093887211#
      Q1 = -.01562499995#: Q2 = .0001430488765#: Q3 = -.000006911147651#
      Q4 = 7.621095161000001D-07: Q5 = -.0000000934945152#
      R1 = 57568490574#: R2 = -13362590354#: R3 = 651619640.7#
      R4 = -11214424.18#: R5 = 77392.33017#: R6 = -184.9052456#
      S1 = 57568490411#: S2 = 1029532985#: S3 = 9494680.718#
      S4 = 59272.64853#: S5 = 267.8532712#: S6 = 1#
      
      IF X = 0# THEN GOTO 1
      AX = ABS(X)
      IF AX < 8 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))))
	BESSJ0 = FR / FS
      ELSE
	Z = 8# / AX
	Y = Z * Z
	XX = AX - .785398164#
	FP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)))
	FQ = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)))
	BESSJ0 = SQR(.636619772# / AX) * (FP * COS(XX) - Z * FQ * SIN(XX))
      END IF
      RETURN
1     BESSJ0 = 1#
RETURN

' ----------------------------------------------------------------------
700   'FUNCTION BESSJ1 (X)
      
'     THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION
'     OF ORDER 0 FOR ANY REAL X. WE USE HERE THE POLYNOMIAL APPROXIMATION
'     BY SERIES OF CHEBYSHEV POLYNOMIALS 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.

      P1 = 1#: P2 = 1.83105E-03: P3 = -.00003516396496#
      P4 = .000002457520174#: P5 = -.000000240337019#: P6 = .636619772#
      Q1 = .04687499995#: Q2 = -.0002002690873#: Q3 = .000008449199096#
      Q4 = -.00000088228987#: Q5 = .000000105787412#
      R1 = 72362614232#: R2 = -7895059235#: R3 = 242396853.1#
      R4 = -2972611.439#: R5 = 15704.4826#: R6 = -30.16036606#
      S1 = 144725228442#: S2 = 2300535178#: S3 = 18583304.74#
      S4 = 99447.43394#: S5 = 376.9991397#: S6 = 1#

      AX = ABS(X)
      IF AX < 8 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))))
	BESSJ1 = X * (FR / FS)
      ELSE
	Z = 8# / AX
	Y = Z * Z
	XX = AX - 2.35619491#
	FP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)))
	FQ = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)))
	IF X >= 0# THEN
	  SIGN = ABS(S6)
	ELSE
	  SIGN = -ABS(S6)
	END IF
	BESSJ1 = SQR(P6 / AX) * (COS(XX) * FP - Z * SIN(XX) * FQ) * SIGN
      END IF
RETURN


800  'FUNCTION BESSY0 (X)

' ----------------------------------------------------------------------
'     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.
' ----------------------------------------------------------------------
      P1 = 1#: P2 = -.001098628627#: P3 = .00002734510407#
      P4 = -.000002073370639#: P5 = .0000002093887211#
      Q1 = -.01562499995#: Q2 = .0001430488765#: Q3 = -.000006911147651#
      Q4 = 7.621095161000001D-07: Q5 = -.0000000934945152#
      R1 = -2957821389#: R2 = 7062834065#: R3 = -512359803.6#
      R4 = 10879881.29#: R5 = -86327.92757#: R6 = 228.4622733#
      S1 = 40076544269#: S2 = 745249964.8#: S3 = 7189466.438#
      S4 = 47447.2647#:  S5 = 226.1030244#: S6 = 1#
      IF X = 0# THEN
	BESSY0 = -1D+30
	RETURN
      END IF
      IF X < 8# THEN
	Y = X * X
	FR1 = R1 + Y * (R2 + Y * (R3 + Y * (R4 + Y * (R5 + Y * R6))))
	FS1 = S1 + Y * (S2 + Y * (S3 + Y * (S4 + Y * (S5 + Y * S6))))
	GOSUB 600
	BESSY0 = FR1 / FS1 + .636619772# * BESSJ0 * LOG(X)
      ELSE
	Z = 8# / X
	Y = Z * Z
	XX = X - .785398164#
	FP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)))
	FQ = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)))
	BESSY0 = SQR(.636619772# / X) * (FP * SIN(XX) + Z * FQ * COS(XX))
      END IF
RETURN

900  'FUNCTION BESSY1 (X)
' -----------------------------------------------------------------------
'     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.
' -----------------------------------------------------------------------
      P1 = 1#: P2 = .00183105#: P3 = -.00003516396496#
      P4 = .000002457520174#: P5 = -.000000240337019#
      Q1 = .04687499995#: Q2 = -.0002002690873#: Q3 = .000008449199096#
      Q4 = -.00000088228987#: Q5 = .000000105787412#
      R1 = -4900604943000#: R2 = 1275274390000#: R3 = -51534381390#
      R4 = 734926455.1#:   R5 = -4237922.726#:  R6 = 8511.937935#
      S1 = 24995805700000#: S2 = 424441966400#: S3 = 3733650367#
      S4 = 22459040.02#:   S5 = 102042.605#:    S6 = 354.9632885#: S7 = 1#
      IF X = 0# THEN
	BESSY1 = -1D+30
	RETURN
      END IF
      IF X < 8# THEN
	Y = X * X
	FR1 = R1 + Y * (R2 + Y * (R3 + Y * (R4 + Y * (R5 + Y * R6))))
	TMP = S6 + Y * S7
	FS1 = S1 + Y * (S2 + Y * (S3 + Y * (S4 + Y * (S5 + Y * TMP))))
	GOSUB 700
	BESSY1 = X * (FR1 / FS1) + .636619772# * (BESSJ1 * LOG(X) - 1# / X)
      ELSE
	Z = 8# / X
	Y = Z * Z
	XX = X - 2.356194491#
	FP = P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)))
	FQ = Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)))
	BESSY1 = SQR(.636619772# / X) * (SIN(XX) * FP + Z * COS(XX) * FQ)
      END IF
RETURN


'End of file Tbessy.bas