{************************************************************** * 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) * * 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". * * * * Pascal Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * **************************************************************} Program msphi; Uses WinCrt; Type VEC = Array[0..250] of double; Const MM = 2.30258509; Var SI, DI: VEC; k, n, nm, ns: Integer; X: double; Function ENVJ(N:integer; X:double): double; Begin ENVJ := (0.5*Ln(6.28*N)-N*Ln(1.36*X/N))/MM End; Function MSTA1(X: double; MP: integer): Integer; { =================================================== ! 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 ! =================================================== } Label 20; Var A0,F,F0,F1:Double; IT,NN,N0,N1:Integer; Begin A0:=abs(X); N0:=Round(1.1*A0)+1; F0:=ENVJ(N0,A0)-MP; N1:=N0+5; F1:=ENVJ(N1,A0)-MP; for IT:=1 to 20 do begin NN:=Round(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 end; 20: MSTA1 := NN End; Function MSTA2(X:double; N, MP: integer): Integer; { =================================================== ! 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 ! =================================================== } Label 20; Var A0,EJN,F,F0,F1,HMP,OBJ:Double; IT,N0,N1,NN: Integer; Begin A0:=abs(X); HMP:=0.5*MP; EJN:=ENVJ(N,A0); if EJN <= HMP then begin OBJ:=MP; N0:=Round(1.1*A0) end else begin OBJ:=HMP+EJN; N0:=N end; F0:=ENVJ(N0,A0)-OBJ; N1:=N0+5; F1:=ENVJ(N1,A0)-OBJ; for IT:=1 to 20 do begin NN:=Round(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 end; 20: MSTA2 := NN+10 End; Function cosh(x:double): double; begin cosh:=(exp(x)+exp(-x))/2.0 end; Function sinh(x:double): double; begin sinh:=(exp(x)-exp(-x))/2.0 end; Procedure SPHI(N:integer; X:double; Var NM:integer; Var SI, DI: VEC); { ======================================================== ! 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 ! ======================================================== } Label 10; Var K, M: Integer; CS, F, F0, F1, SI0: Double; Begin NM:=N; if abs(X) < 1E-100 then begin for K:=0 to N do begin SI[K]:=0.0; DI[K]:=0.0 end; SI[0]:=1.0; DI[1]:=0.333333333333333; goto 10 end; SI[0]:=sinh(X)/X; SI[1]:=-(sinh(X)/X-cosh(X))/X; SI0:=SI[0]; if N >= 2 then begin M:=MSTA1(X,200); if M < N then NM:=M else M:=MSTA2(X,N,15); F0:=0.0; F1:=1.0-100; for K:=M Downto 0 do begin F:=(2.0*K+3.0)*F1/X+F0; if K <= NM then SI[K]:=F; F0:=F1; F1:=F end; CS:=SI0/F; for K:=0 to NM do SI[K] := SI[K]*CS end; DI[0]:=SI[1]; for K:=1 to NM do DI[K]:=SI[K-1]-(K+1.0)/X*SI[K]; 10: End; {main program} BEGIN InitWinCrt; writeln; write(' Please enter n and x: '); readln(n, X); if n <= 10 then ns:=1 else begin write(' Please enter order step Ns: '); readln(ns) end; SPHI(n,X,nm,SI,DI); writeln; writeln(' n in(x) in''(x)'); writeln('--------------------------------------------'); k:=0; Repeat writeln(k:3, SI[k]:18:6, DI[k]:18:6); k:=k+ns until k = nm+1; ReadKey; DoneWinCrt END. {end of file msphi.pas}