I am trying to compute the Bessel function of the second kind (Bessel_y) using the SLATEC's Amos library available on Netlib. Here is the SLATEC code I use. Below I have pasted my test program that calls SLATEC routine CBESY.
PROGRAM BESSELTEST
IMPLICIT NONE
REAL:: FNU
INTEGER, PARAMETER :: N = 2, KODE = 1
COMPLEX,ALLOCATABLE :: CWRK (:), CY (:)
COMPLEX:: Z, ci
INTEGER :: NZ, IERR
ALLOCATE(CWRK(N), CY(N))
ci = cmplx (0.0, 1.0)
FNU = 0.0e0
Z = CMPLX(0.3e0, 0.4e0)
CALL CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)
WRITE(*,*) 'CY: ', CY
WRITE(*,*) 'IERR: ', IERR
STOP
END PROGRAM
And here is the output of the above program:
CY: ( 5.78591091E-39, 5.80327020E-39) ( 0.0000000 , 0.0000000 )
IERR: 4
Ierr = 4 meaning there is some problem with the input itself. To be precise, the IERR = 4 means the following as per the header info in CBESY.f file:
! IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
! TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
! CANCE BY ARGUMENT REDUCTION
Clearly, CABS(Z) (which is 0.50) or FNU + N - 1 (which is 1.0) are not too large but still the routine CBESY throws the error message number 4 as above. The CY array should have following values for the argument given in above code:
CY(1) = -0.4983 + 0.6700i
CY(2) = -1.0149 + 0.9485i
These values are computed using Matlab.
I can't figure out what's the problem when I call CBESY from SLATEC library. Any clues? Much thanks for the suggestions/help.
PS: if it is of any help, I used gfortran to compile, link and then create the SLATEC library file ( the .a file ) which I keep in the same directory as my test program above.
shell command to execute above code:
gfortran -c BesselTest.f95
gfortran -o a *.o libslatec.a
a
GD.
