!*************************************************************** !* 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". * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !*************************************************************** PROGRAM MSPHY IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SY(0:250),DY(0:250) WRITE(*,*)'Please enter n and x ' READ(*,*)N,X WRITE(*,30)N,X IF (N.LE.10) THEN NS=1 ELSE WRITE(*,*)'Please enter order step Ns' READ(*,*)NS ENDIF CALL SPHY(N,X,NM,SY,DY) WRITE(*,*) WRITE(*,*)' n yn(x) yn''(x)' WRITE(*,*)'--------------------------------------------' DO 10 K=0,NM,NS 10 WRITE(*,20)K,SY(K),DY(K) 20 FORMAT(1X,I3,2D20.10) 30 FORMAT(3X,6HNmax =,I3,', ',2Hx=,F6.1) END SUBROUTINE 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 ! ====================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SY(0:N),DY(0:N) NM=N IF (X.LT.1.0D-60) THEN DO 10 K=0,N SY(K)=-1.0D+300 10 DY(K)=1.0D+300 RETURN ENDIF SY(0)=-DCOS(X)/X SY(1)=(SY(0)-DSIN(X))/X F0=SY(0) F1=SY(1) DO 15 K=2,N F=(2.0D0*K-1.0D0)*F1/X-F0 SY(K)=F IF (DABS(F).GE.1.0D+300) GO TO 20 F0=F1 15 F1=F 20 NM=K-1 DY(0)=(DSIN(X)+DCOS(X)/X)/X DO 25 K=1,NM 25 DY(K)=SY(K-1)-(K+1.0D0)*SY(K)/X RETURN END !end of file msphy.f90