!******************************************************************
!*       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".          *
!*                                                                *
!*                            F90 Release By J-P Moreau, Paris.   *
!*                                   (www.jpmoreau.fr)            *
!******************************************************************
        PROGRAM MSPHK
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION SK(0:250),DK(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 SPHK(N,X,NM,SK,DK)
        WRITE(*,*)
        WRITE(*,*)'  n          kn(x)               kn''(x)'
        WRITE(*,*)'--------------------------------------------'
        DO 10 K=0,NM,NS
10         WRITE(*,20)K,SK(K),DK(K)
20      FORMAT(1X,I3,2D20.10)
30      FORMAT(3X,'Nmax =',I3,',     ','x =',F6.1)
        END


        SUBROUTINE 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
!       =====================================================

        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION SK(0:N),DK(0:N)
        PI=3.141592653589793D0
        NM=N
        IF (X.LT.1.0D-60) THEN
           DO 10 K=0,N
              SK(K)=1.0D+300
10            DK(K)=-1.0D+300
           RETURN
        ENDIF
        SK(0)=0.5D0*PI/X*DEXP(-X)
        SK(1)=SK(0)*(1.0D0+1.0D0/X)
        F0=SK(0)
        F1=SK(1)
        DO 15 K=2,N
           F=(2.0D0*K-1.0D0)*F1/X+F0
           SK(K)=F
           IF (DABS(F).GT.1.0D+300) GO TO 20
           F0=F1
15         F1=F
20      NM=K-1
        DK(0)=-SK(1)
        DO 25 K=1,NM
25         DK(K)=-SK(K-1)-(K+1.0D0)/X*SK(K)
        RETURN
        END

!end of file msphk.f90