DECLARE SUB SPHK (N%, X!, NM%, SK!(), DK!()) '****************************************************************** '* Purpose: This program computes the modified spherical * '* Bessel functions kn(x) and kn'(x) using * '* subroutine SPHK * '* Input : x --- Argument of kn(x) ( x > 0 ) * '* n --- Order of kn(x) ( n = 0 to 250 ) * '* Output: SK(n) --- kn(x) * '* DK(n) --- kn'(x) * '* Example: x= 10.0 * '* n kn(x) kn'(x) * '* -------------------------------------------- * '* 0 .7131404291D-05 -.7844544720D-05 * '* 1 .7844544720D-05 -.8700313235D-05 * '* 2 .9484767707D-05 -.1068997503D-04 * '* 3 .1258692857D-04 -.1451953914D-04 * '* 4 .1829561771D-04 -.2173473743D-04 * '* 5 .2905298451D-04 -.3572740841D-04 * '* -------------------------------------------------------------- * '* 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 SK(250), DK(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 SPHK(N, X, NM, SK(), DK()) PRINT PRINT " n kn(x) kn'(x)" PRINT "--------------------------------------------" FOR K = 0 TO NM STEP NS PRINT K; TAB(12); SK(K); TAB(28); DK(K) NEXT K END SUB SPHK (N, X, NM, SK(), DK()) ' ===================================================== ' Purpose: Compute modified spherical Bessel functions ' of the second kind, kn(x) and kn'(x) ' Input : x --- Argument of kn(x) ( x ò 0 ) ' n --- Order of kn(x) ( n = 0,1,2,... ) ' Output: SK(n) --- kn(x) ' DK(n) --- kn'(x) ' NM --- Highest order computed ' ===================================================== PI = 3.141592653589793# NM = N IF X < 0! THEN FOR K = 0 TO N SK(K) = 1.0E+300 DK(K) = -1.0E+300 NEXT K GOTO 20 END IF SK(0) = .5 * PI / X * EXP(-X) SK(1) = SK(0) * (1! + 1! / X) F0 = SK(0) F1 = SK(1) FOR K = 2 TO N F = (2! * K - 1!) * F1 / X + F0 SK(K) = F If ABS(F) > 1.0E+300 Then GoTo 10 F0 = F1 F1 = F NEXT K 10 NM = K - 1 DK(0) = -SK(1) FOR K = 1 TO NM DK(K) = -SK(K - 1) - (K + 1!) / X * SK(K) NEXT K 20 End Sub 'end of file msphk.bas