DECLARE SUB SPHK (N%, X!, NM%, SK!(), DK!()) '*************************************************************** '* Purpose: This program computes the spherical Bessel * '* functions yn(x) and yn'(x) using subroutine * '* SPHY * '* Input : x --- Argument of yn(x) ( x > 0 ) * '* n --- Order of yn(x) ( n = 0 to 250 ) * '* Output: SY(n) --- yn(x) * '* DY(n) --- yn'(x) * '* Example: x = 10.0 * '* n yn(x) yn'(x) * '* -------------------------------------------- * '* 0 .8390715291D-01 -.6279282638D-01 * '* 1 .6279282638D-01 .7134858763D-01 * '* 2 -.6506930499D-01 .8231361788D-01 * '* 3 -.9532747888D-01 -.2693831344D-01 * '* 4 -.1659930220D-02 -.9449751377D-01 * '* 5 .9383354168D-01 -.5796005523D-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 MSPHY OPTION BASE 0 DEFINT I-N DIM SY(250), DY(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 SPHY(N, X, NM, SY(), DY()) PRINT PRINT " n yn(x) yn'(x)" PRINT "--------------------------------------------" FOR K = 0 TO NM STEP NS PRINT K; TAB(12); SY(K); TAB(28); DY(K) NEXT K END Sub SPHY(N, X, NM, SY(), DY()) ' ====================================================== ' Purpose: Compute spherical Bessel functions yn(x) and ' their(derivatives) ' Input : x --- Argument of yn(x) ( x ע 0 ) ' n --- Order of yn(x) ( n = 0,1,תתת ) ' Output: SY(n) --- yn(x) ' DY(n) --- yn'(x) ' NM --- Highest order computed ' ====================================================== NM = N If X < 1.0E-60 Then For K = 0 To N SY(K) = -1.0E+300 DY(K) = 1.0E+300 Next K GoTo 20 End If SY(0) = -DCOS(X) / X SY(1) = (SY(0) - DSIN(X)) / X F0 = SY(0) F1 = SY(1) For K = 2 To N F = (2.0 * K - 1.0) * F1 / X - F0 SY(K) = F IF ABS(F) >= 1E+300) Then GoTo 10 F0 = F1 F1 = F NEXT K 10: NM = K - 1 DY(0) = (SIN(X) + COS(X) / X) / X For K = 1 To NM DY(K)=SY(K-1)-(K+1.0D0)*SY(K)/X Next K 20 End Sub