'*************************************************************** '* Purpose: This program computes the modified spherical * '* Bessel functions of the first kind in(x) and * '* in'(x) using subroutine SPHI * '* ----------------------------------------------------------- * '* Input : x --- Argument of in(x) (x > 0 ) * '* n --- Order of in(x) ( 0 to 250 ) * '* Output: SI(n) --- in(x) * '* DI(n) --- in'(x) * '* Example: x = 10.0 * '* n in(x) in'(x) * '* -------------------------------------------- * '* 0 .1101323287D+04 .9911909633D+03 * '* 1 .9911909633D+03 .9030850948D+03 * '* 2 .8039659985D+03 .7500011637D+03 * '* 3 .5892079640D+03 .5682828129D+03 * '* 4 .3915204237D+03 .3934477522D+03 * '* 5 .2368395827D+03 .2494166741D+03 * '* ----------------------------------------------------------- * '* 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) * '*************************************************************** OPTION BASE 0 DEFINT I-N DIM SI(250), DI(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 SPHI(N, X, NM, SI(), DI()) PRINT PRINT " n in(x) in'(x)" PRINT "--------------------------------------------" FOR K = 0 TO NM STEP NS PRINT K, SI(K), DI(K) NEXT K END SUB SPHI (N, X, NM, SI(), DI()) ' Purpose: Compute modified spherical Bessel functions ' of the first kind, in(x) and in'(x) ' Input : x --- Argument of in(x) ' n --- Order of in(x) ( n = 0,1,2,... ) ' Output: SI(n) --- in(x) ' DI(n) --- in'(x) ' NM --- Highest order computed ' Routines called: ' MSTA1 and MSTA2 for computing the starting ' point for backward recurrence ' ======================================================== NM = N IF ABS(X) < 1E-100 THEN FOR K = 0 TO N SI(K) = 0.0 DI(K) = 0.0 NEXT K SI(0) = 1.0 DI(1) = 0.333333333333333 GOTO 20 END IF SI(0) = SINH(X) / X SI(1) = -(SINH(X) / X - COSH(X)) / X SI0 = SI(0) IF N >= 2 THEN M = MSTA1(X, 200) IF M < N THEN NM = M ELSE M = MSTA2(X, N, 15) END IF F0 = 0.0 F1 = 1.0 - 100 FOR K = M TO 0 STEP -1 F = (2.0 * K + 3.0) * F1 / X + F0 IF K <= NM THEN SI(K) = F F0 = F1 F1 = F NEXT K CS = SI0 / F FOR K = 0 TO NM SI(K) = CS * SI(K) NEXT K END IF DI(0) = SI(1) FOR K = 1 TO NM DI(K) = SI(K - 1) - (K + 1.0) / X * SI(K) NEXT K 20 END SUB 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.0 - F0 / F1) F = ENVJ(NN, A0) - MP IF ABS(NN - N1) < 1 THEN GOTO 20 N0 = N1 F0 = F1 N1 = NN F1 = F NEXT IT 20 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 = 0.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.0 - F0 / F1) F = ENVJ(NN, A0) - OBJ IF ABS(NN - N1) < 1 THEN GOTO 20 N0 = N1 F0 = F1 N1 = NN F1 = F NEXT IT 20 MSTA2 = NN + 10 END FUNCTION FUNCTION ENVJ (N, X) ENVJ = 0.5 * LOG(6.28 * N) - N * LOG(1.36 * X / N) END FUNCTION FUNCTION COSH (X) COSH = (EXP(X) + EXP(-X)) / 2 END FUNCTION FUNCTION SINH (X) SINH = (EXP(X) - EXP(-X)) / 2 END FUNCTION 'end of file msphi.bas