'************************************************************** '* CALCULATE THE Kth ZERO OF THE DERIVATIVE OF BESSEL * '* FUNCTION OF ORDER N, J(N,X) * '* ---------------------------------------------------------- * '* SAMPLE RUN: * '* (Calculate the 10th zero of the derivative of J(1,X) ). * '* * '* X= 30.60192297351452 * '* * '* BESSJP(X)= 4.236689395576559D-17 * '* * '* ---------------------------------------------------------- * '* Reference: From Numath Library By Tuan Dang Trong in * '* Fortran 77 [BIBLI 18]. * '* * '* Basic Release By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '************************************************************** DEFDBL A-H, O-Z DEFINT I-N ITRUE = 1 N = 1 K = 10 GOSUB 1000 'call ZEROJP(N,K) CLS PRINT PRINT " X="; ZEROJP PRINT X = ZEROJP: GOSUB 1500 PRINT " BESSJP(X)="; BESSJP PRINT END 500 'Sign(S6, X) IF X < 0 THEN Sign = -ABS(S6) ELSE Sign = ABS(S6) END IF RETURN 1000 'ZEROJP(N, K) '-------------------------------------------------------------------- ' CALCULATE THE Kth ZERO OF THE DERIVATIVE OF BESSEL FUNCTION ' OF ORDER N, J(N,X) '-------------------------------------------------------------------- ' CALLING MODE: ' RES = ZEROJP(N,K) ' ' INPUTS: ' N ORDER OF BESSEL FUNCTION J (INTEGER >= 0) I*4 ' K RANK OF ZERO (INTEGER > 0) I*4 ' OUTPUT: ' ZEROJP R*8 ' REFERENCE: ' ABRAMOWITZ M. & STEGUN IRENE A. ' HANDBOOK OF MATHEMATICAL FUNCTIONS '-------------------------------------------------------------------- ' Labels: 10, 20, 25 ' boolean IMPROV PI = 4# * ATN(1#) TOL = .0000001: NITMX = 15 C1 = .8086165#: C2 = .07249: C3 = .05097: C4 = .0094 IMPROV = ITRUE FFN = 1# * N 'FN is a reserved word in Basic FK = 1# * K IF K > 1 THEN GOTO 10 ' if N = 0 and K = 1 IF N = 0 THEN ZEROJP = 0# RETURN ' TCHEBYCHEV'S SERIES FOR K <= N ELSE F1 = FFN ^ (1# / 3#) F2 = F1 * F1 * FFN ZEROJP = FFN + C1 * F1 + (C2 / F1) - (C3 / FFN) + (C4 / F2) GOTO 20 END IF ' MAC MAHON'S SERIES FOR K >> N 10 B0 = (FK + .5 * FFN - .75) * PI B1 = 8# * B0 B2 = B1 * B1 B3 = 3# * B1 * B2 B5 = 5# * B3 * B2 B7 = 7# * B5 * B2 T0 = 4# * FFN * FFN T1 = T0 + 3# T3 = 4# * ((7# * T0 + 82#) * T0 - 9#) T5 = 32# * (((83# * T0 + 2075#) * T0 - 3039#) * T0 + 3537#) T7 = 64# * ((((6949# * T0 + 296492#) * T0 - 1248002#) * T0 + 7414380#) * T0 - 5853627#) ZEROJP = B0 - (T1 / B1) - (T3 / B3) - (T5 / B5) - (T7 / B7) 20 IF IMPROV <> 0 THEN ' IMPROVE SOLUTION BY SECANT METHOD WHEN K > N ' AND IMPROV = TRUE P0 = .9 * ZEROJP P1 = ZEROJP IER = 0 NEV = 2 X = P0: GOSUB 1500 Q0 = BESSJP X = P1: GOSUB 1500 Q1 = BESSJP FOR IT = 1 TO NITMX P = P1 - Q1 * (P1 - P0) / (Q1 - Q0) DP = P - P1 IF IT = 1 THEN GOTO 25 IF ABS(DP) < TOL THEN ZEROJP = P RETURN END IF 25 NEV = NEV + 1 P0 = P1 Q0 = Q1 P1 = P X = P1: GOSUB 1500 Q1 = BESSJP NEXT IT IER = 1 PRINT " ** ZEROJP ** NITMX EXCEEDED." END IF RETURN 1500 'BESSJP(N,X) ' ---------------------------------------------------------------------- ' NAME : BESSJP ' DATE : 06/01/1982 ' IV : 1 ' IE : 1 ' AUTHOR: DANG TRONG TUAN ' ...................................................................... ' ' FIRST DERIVATIVE OF FIRST KIND BESSEL FUNCTION OF ORDER N, FOR REAL X ' ' MODULE BESSJP ' ...................................................................... ' ' THIS SUBROUTINE CALCULATES THE FIRST DERIVATIVE OF FIRST KIND BESSEL ' FUNCTION OF ORDER N, FOR REAL X. ' ' ...................................................................... ' ' I VARIABLE DIMENSION/TYPE DESCRIPTION (INPUTS) ' N I*4 ORDER OF FUNCTION ' X R*8 ABSCISSA OF FUNCTION BESSJP(N,X) ' ' O VARIABLE,DIMENSION/TYPE DESCRIPTION (OUTPUT) ' ' BESSJP R*8 FUNCTION EVALUATION AT X '....................................................................... ' CALLED SUBROUTINE ' ' BESSJ FIRST KIND BESSEL FUNCTION ' ' ---------------------------------------------------------------------- IF N = 0 THEN NN = 1: GOSUB 1600 BESSJP = -BESSJ ELSEIF X = 0# THEN BESSJP = 0# ELSE NN = N - 1: GOSUB 1600 BESSJP = BESSJ NN = N: GOSUB 1600 BESSJP = BESSJP - (1# * N / X) * BESSJ END IF RETURN 1600 'BESSJ(NN, X) '------------------------------------------------------------------------ ' THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ' NN, INTEGER FOR ANY REAL X. THE CLASSICAL RECURRENCE FORMULA IS USED ' HERE, WHEN X IS > N, FOR X < N, THE MILLER'S ALGORITHM ALLOWS ' AVOIDING OVERFLOWS. ' REFERENCE: ' C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, ' MATHEMATICAL TABLES, VOL.5, 1962. '------------------------------------------------------------------------ IACC = 40: BIGNO = 1E+10: BIGNI = 1E-10 IF NN = 0 THEN GOSUB 1610 BESSJ = BESSJ0 RETURN END IF IF NN = 1 THEN GOSUB 1620 BESSJ = BESSJ1 RETURN END IF IF X = 0# THEN BESSJ = 0# RETURN END IF TOX = 2# / X IF X > 1# * NN THEN GOSUB 1610: BJM = BESSJ0 GOSUB 1620: BJ = BESSJ1 FOR J = 1 TO NN - 1 BJP = J * TOX * BJ - BJM BJM = BJ BJ = BJP NEXT J BESSJ = BJ RETURN ELSE M = INT(2 * ((NN + INT(SQR(IACC * NN))) / 2)) BSJ = 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 BSJ = BSJ * BIGNI SUM = SUM * BIGNI END IF IF JSUM <> 0 THEN SUM = SUM + BJ JSUM = 1 - JSUM IF J = NN THEN BSJ = BJP NEXT J SUM = 2# * SUM - BJ BESSJ = BSJ / SUM END IF RETURN 1610 'BESSJ0(X) '---------------------------------------------------------------------- ' THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ' ZERO FOR ANY REAL X. THE CHEBYSHEV'S POLYNOMIAL SERIES IS USED ' FOR 0