DECLARE SUB SPHJ (N%, X!, NM%, SJ!(), DJ!()) DECLARE FUNCTION ENVJ! (N%, X!) DECLARE FUNCTION MSTA1% (X!, MP%) DECLARE FUNCTION MSTA2% (X!, N%, MP%) '**************************************************************** '* Purpose: This program computes the spherical Bessel * '* functions jn(x) and jn'(x) using subroutine * '* SPHJ * '* Input : x --- Argument of jn(x) ( x > 0) * '* n --- Order of jn(x) ( n = 0 to 250 ) * '* Output: SJ(n) --- jn(x) * '* DJ(n) --- jn'(x) * '* Example: x =10.0 * '* n jn(x) jn'(x) * '* -------------------------------------------- * '* 0 -.5440211109D-01 -.7846694180D-01 * '* 1 .7846694180D-01 -.7009549945D-01 * '* 2 .7794219363D-01 .5508428371D-01 * '* 3 -.3949584498D-01 .9374053162D-01 * '* 4 -.1055892851D+00 .1329879757D-01 * '* 5 -.5553451162D-01 -.7226857814D-01 * '* ------------------------------------------------------------ * '* REFERENCE: "Fortran Routines for Computation of Special * '* Functions, * '* jin.ece.uiuc.edu/routines/routines.html". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '**************************************************************** 'PROGRAM MSPHJ OPTION BASE 0 DEFINT I-N DIM SJ(250), DJ(250) PRINT INPUT " Please enter n and x: ", N, X IF N <= 10 THEN NS = 1 ELSE INPUT " Please enter order step Ns: ", NS END IF CALL SPHJ(N, X, NM, SJ(), DJ()) PRINT PRINT " n jn(x) jn'(x)" PRINT "--------------------------------------------" FOR K = 0 TO NM STEP NS PRINT K; TAB(12); SJ(K); TAB(28); DJ(K) NEXT K END 'end of file msphj.bas FUNCTION ENVJ (N, X) ENVJ = .5 * LOG(6.28 * N) - N * LOG(1.36 * X / N) END FUNCTION FUNCTION MSTA1 (X, MP) ' =================================================== ' Purpose: Determine the starting point for backward ' recurrence such that the magnitude of ' Jn(x) at that point is about 10^(-MP) ' Input : x --- Argument of Jn(x) ' MP --- Value of magnitude ' Output: MSTA1 --- Starting point ' =================================================== A0 = ABS(X) N0 = INT(1.1 * A0) + 1 F0 = ENVJ(N0, A0) - MP N1 = N0 + 5 F1 = ENVJ(N1, A0) - MP FOR IT = 1 TO 20 NN = N1 - (N1 - N0) / (1! - F0 / F1) F = ENVJ(NN, A0) - MP IF ABS(NN - N1) < 1 THEN GOTO 30 N0 = N1 F0 = F1 N1 = NN F1 = F NEXT IT 30 MSTA1 = NN END FUNCTION FUNCTION MSTA2 (X, N, MP) ' =================================================== ' Purpose: Determine the starting point for backward ' recurrence such that all Jn(x) has MP ' significant digits ' Input : x --- Argument of Jn(x) ' n --- Order of Jn(x) ' MP --- Significant digit ' Output: MSTA2 --- Starting point ' =================================================== A0 = ABS(X) HMP = .5 * MP EJN = ENVJ(N, A0) IF (EJN.LE.HMP) THEN OBJ = MP N0 = INT(1.1 * A0) ELSE OBJ = HMP + EJN N0 = N END IF F0 = ENVJ(N0, A0) - OBJ N1 = N0 + 5 F1 = ENVJ(N1, A0) - OBJ FOR IT = 1 TO 20 NN = N1 - (N1 - N0) / (1! - F0 / F1) F = ENVJ(NN, A0) - OBJ IF ABS(NN - N1) < 1 THEN GOTO 10 N0 = N1 F0 = F1 N1 = NN F1 = F NEXT IT 10 MSTA2 = NN + 10 END FUNCTION SUB SPHJ (N, X, NM, SJ(), DJ()) ' ======================================================= ' Purpose: Compute spherical Bessel functions jn(x) and ' their derivatives ' Input : x --- Argument of jn(x) ' n --- Order of jn(x) ( n = 0,1,תתת ) ' Output: SJ(n) --- jn(x) ' DJ(n) --- jn'(x) ' NM --- Highest order computed ' Routines called: ' MSTA1 and MSTA2 for computing the starting ' point for backward recurrence ' ======================================================= NM = N IF ABS(X) < 0! THEN FOR K = 0 TO N SJ(K) = 0! DJ(K) = 0! NEXT K SJ(0) = 1! DJ(1) = .3333333333333333# GOTO 20 END IF SJ(0) = SIN(X) / X SJ(1) = (SJ(0) - COS(X)) / X IF N >= 2 THEN SA = SJ(0) SB = SJ(1) M = MSTA1(X, 200) IF M < N THEN NM = M ELSE M = MSTA2(X, N, 15) END IF F0 = 0! F1 = 1! - 100 FOR K = M TO 0 STEP -1 F = (2! * K + 3!) * F1 / X - F0 IF K <= NM THEN SJ(K) = F F0 = F1 F1 = F NEXT K IF ABS(SA) > ABS(SB) THEN CS = SA / F IF ABS(SA) <= ABS(SB) THEN CS = SB / F0 FOR K = 0 TO NM SJ(K) = CS * SJ(K) NEXT K END IF DJ(0) = (COS(X) - SIN(X) / X) / X FOR K = 1 TO NM DJ(K) = SJ(K - 1) - (K + 1!) * SJ(K) / X NEXT K 20 END SUB