' --------------------------------------------------------------------------
' Program to calculate the value of the first kind Bessel function of integer
' order N using the subroutine BESSJ.
' --------------------------------------------------------------------------
' SAMPLE RUN:
'
' (Calculate the 1st kind Bessel function of order 2 for X=8.41724462111).
'
'  ORDER = 2
'
'  X = 8.41724462111
'
'  BESSJ OF ORDER 2 = 1.25443232...D-07
'
' --------------------------------------------------------------------------
' Reference: From Numath Library By Tuan Dang Trong in Fortran 77.
'
'                               Basic Release 1.0 By J-P Moreau, Paris
'                                         (www.jpmoreau.fr)
' --------------------------------------------------------------------------
'PROGRAM TBESSJ

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

  CLS

  PRINT
  PRINT " BESSEL FUNCTIONS OF FIRST KIND"

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

  GOSUB 500  'CALL BESSJ

  PRINT
  PRINT " BESSJ OF ORDER "; N; " = "; BESSJ
  PRINT
  
END

' ----------------------------------------------------------------------
500   'FUNCTION BESSJ (N,X)

'     THIS FUNCTION RETURNS THE VALUE OF THE FIRST KIND BESSEL FUNCTION
'     OF ORDER N, INTEGER FOR ANY REAL X. WE USE HERE THE CLASSICAL
'     RECURRENT FORMULA, WHEN  X > N. FOR X < N, THE MILLER'S ALGORITHM
'     IS USED TO AVOID OVERFLOWS.
'     REFERENCE :
'     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
'     MATHEMATICAL TABLES, VOL.5, 1962.

      IACC = 40
      BIGNO = 1E+10
      BIGNI = 1E-10
      
      IF N = 0 THEN
	'BESSJ = BESSJ0(X)
	GOSUB 600: BESSJ = BESSJ0
	RETURN
      END IF
      IF N = 1 THEN
	'BESSJ = BESSJ1(X)
	GOSUB 700: BESSJ = BESSJ1
	RETURN
      END IF
      IF X = 0# THEN
	BESSJ = 0#
	RETURN
      END IF
      
      TOX = 2# / X
      IF X > 1# * N THEN
	'BJM = BESSJ0(X)
	GOSUB 600: BJM = BESSJ0
	'BJ  = BESSJ1(X)
	GOSUB 700: BJ = BESSJ1
	FOR J = 1 TO N - 1
	  BJP = J * TOX * BJ - BJM
	  BJM = BJ
	  BJ = BJP
	NEXT J
	BESSJ = BJ
      ELSE
	M = 2 * ((N + INT(SQR(1# * IACC * N))) / 2)
	BESSJ = 0#
	JSUM = 0
	SUM = 0#
	BJP = 0#
	BJ = 1#
	FOR J = M TO 1 STEP -1
	  BJM = J * TOX * BJ - BJP
	  BJP = BJ
	  BJ = BJM
	  IF ABS(BJ) > BIGNO THEN
	    BJ = BJ * BIGNI
	    BJP = BJP * BIGNI
	    BESSJ = BESSJ * BIGNI
	    SUM = SUM * BIGNI
	  END IF
	  IF JSUM <> 0 THEN SUM = SUM + BJ
	  JSUM = 1 - JSUM
	  IF J = N THEN BESSJ = BJP
	NEXT J
	SUM = 2# * SUM - BJ
	BESSJ = BESSJ / SUM
      END IF

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

' end of file tbessj.bas