I wrote a simply Fortran code that uses a Bessel function subroutine to solve a fluid-structure interaction problem. the bessel function found are 8-dimensional vector of complex*16 type. The first four elements of the found vectors give NAN and for the last four elements give fair values. the same code is written under intel fortran and it gives all the numerical values of the components of the bessel function vector.
INTEGER N, NP1, IER
DOUBLE PRECISION, PARAMETER :: EPS = 10.0D-12
COMPLEX*16 :: AUX, AUX2, AUX3
COMPLEX*16 :: VECT_AUX(8), VECT_AUX2(8), VECT_AUX3(8), RQ(8), LAMDAS(8)
CALL BESSJ(AUX,N,AUX2,EPS,IER)
VECT_AUX2(JK)=AUX2
CALL BESSJ(AUX,NP1,AUX3,EPS,IER)
VECT_AUX3(JK)=AUX3
RQ(JK)=1.0/(PN-AUX*(AUX3/AUX2)) ! SOLUTION OF BEESEL FUNCTION
The results given using simply fortran:
!====== VECT_AUX ======
( -52.078796118420847 , 62.859646088924379 )
( -52.078796118420847 , -62.859646088924379 )
( 52.078796118420847 , 62.859646088924379 )
( 52.078796118420847 , -62.859646088924379 )
( -4.8816792606421169 , 5.8903897259013362 )
( -4.8816792606421169 , -5.8903897259013362 )
( 4.8816792606421169 , 5.8903897259013362 )
( 4.8816792606421169 , -5.8903897259013362 )
!====== VECT_AUX2 ======
( NaN, NaN)
( NaN, NaN)
( NaN, NaN)
( NaN, NaN)
( 2.2275304232763871E-011, 1.4362161577399207E-011)
( 2.2275304232763871E-011, -1.4362161577399207E-011)
( -2.2275304232763871E-011, 1.4362161577399207E-011)
( -2.2275304232763871E-011, -1.4362161577399207E-011)
!====== VECT_AUX3 ======
( NaN, NaN)
( NaN, NaN)
( NaN, NaN)
( NaN, NaN)
( -3.6772018789810943E-012, 1.2444292792182487E-012)
( -3.6772018789810943E-012, -1.2444292792182487E-012)
( -3.6772018789810943E-012, -1.2444292792182487E-012)
( -3.6772018789810943E-012, 1.2444292792182487E-012)
! ============ RQ ============
( NaN, NaN)
( NaN, NaN)
( NaN, NaN)
( NaN, NaN)
( 0.18338870108203947 ,-0.87337225864044421 )
( 0.18338870108203947 , 0.87337225864044421 )
( 0.18338870108203947 , 0.87337225864044421 )
( 0.18338870108203947 ,-0.87337225864044421 )
whereas with intel fortran i get
! ====== VECT_AUX ======
-.520788E+02 0.628596E+02
-.520788E+02 -.628596E+02
0.520788E+02 0.628596E+02
0.520788E+02 -.628596E+02
-.488168E+01 0.589039E+01
-.488168E+01 -.589039E+01
0.488168E+01 0.589039E+01
0.488168E+01 -.589039E+01
! ====== VECT_AUX2 ======
0.313633E+25 -.327621E+25
0.313633E+25 0.327621E+25
-.313633E+25 -.327621E+25
-.313633E+25 0.327621E+25
0.222753E-10 0.143622E-10
0.222753E-10 -.143622E-10
-.222753E-10 0.143622E-10
-.222753E-10 -.143622E-10
! ====== VECT_AUX3 ======
0.203881E+25 0.291338E+25
0.203881E+25 -.291338E+25
0.203881E+25 -.291338E+25
0.203881E+25 0.291338E+25
-.367720E-11 0.124443E-11
-.367720E-11 -.124443E-11
-.367720E-11 -.124443E-11
-.367720E-11 0.124443E-11
! ============ RQ ============
0.969849E-02 -.737084E-02
0.969849E-02 0.737084E-02
0.969849E-02 0.737084E-02
0.969849E-02 -.737084E-02
0.395602E-01 -.171949E-02
0.395602E-01 0.171949E-02
0.395602E-01 0.171949E-02
0.395602E-01 -.171949E-02
I wonder why with simply fortran I could not get the same results for the first four components of the bessel vector.
PROGRAM BESSEL
INTEGER :: N, NP1, IER, JK
DOUBLE PRECISION, PARAMETER :: EPS = 10.0D-12
REAL*8 :: DENSITYI
COMPLEX*16 :: COMPI, AUX, AUX2, AUX3
COMPLEX*16 :: VECT_AUX(8), VECT_AUX2(8), VECT_AUX3(8), RQ(8), LAMDAS(8)
N = 25
NP1 = N + 1
COMPI = (0.0D0,1.0D0)
DENSITYI = 0.100000E+04 ! INPUT
LAMDAS(1) = ( 62.859646088924379 , 52.078796118420847 )
LAMDAS(2) = ( -62.859646088924379 , 52.078796118420847 )
LAMDAS(3) = ( 62.859646088924379 , -52.078796118420847 )
LAMDAS(4) = ( -62.859646088924379 , -52.078796118420847 )
LAMDAS(5) = ( 5.8903897259013362 , 4.8816792606421169 )
LAMDAS(6) = ( -5.8903897259013362 , 4.8816792606421169 )
LAMDAS(7) = ( 5.8903897259013362 , -4.8816792606421169 )
LAMDAS(8) = ( -5.8903897259013362 , -4.8816792606421169 )
DO 2000 JK=1,8
AUX= LAMDAS(JK)*COMPI
VECT_AUX(JK)=AUX
IF(DENSITYI.NE.0.0) GO TO 2010
RQ(JK)=(0.0,0.0)
GO TO 2000
2010 CALL BESSJ(AUX,N,AUX2,EPS,IER)
VECT_AUX2(JK)=AUX2
CALL BESSJ(AUX,NP1,AUX3,EPS,IER)
VECT_AUX3(JK)=AUX3
RQ(JK)=1.0/(N-AUX*(AUX3/AUX2)) ! SOLUTION OF BEESEL FUNCTION
2000 CONTINUE
PRINT*,'======= VECT_AUX ======'
DO JK=1,8
PRINT*, VECT_AUX(JK)
END DO
PRINT*,'======= VECT_AUX2 ======'
DO JK=1,8
PRINT*, VECT_AUX2(JK)
END DO
PRINT*,'======= VECT_AUX3 ======'
DO JK=1,8
PRINT*, VECT_AUX3(JK)
END DO
PRINT*,'======= RQ ======'
DO JK=1,8
PRINT*, RQ(JK)
END DO
END PROGRAM BESSEL
!************************************************************************
SUBROUTINE BESSJ (Z,N,BJ,EPS,IER)
! BUT = CALCULER LA FONCTION BESSEL J(Z)
! USAGE = CALL BESSJ (Z,N,BJ,EPS,IER)
! PARAMETRES
! Z = ARGUMENT (COMPLEXE)
! N = ORDRE DEMANDE. FAUT N GE 0
! BJ= LE RESULTAT DE LA FONCTION BESSEL J
! EPS=PRECISION DEMANDEE, LE CALCUL SE TERMINE LORSQUE L"ERREUR
! RELATIVE EPS A ETE ATTEINTE.
! IER=CODE D"ERREUR
! =0 CORRECT
! =1 N LT 0
! =3 PRECISION NON ATTEINTE
!
!
!***********************************************************************
IMPLICIT REAL*8(A-H,O-Z)
COMPLEX*16 Z, BJ, BJ1
!
!
BJ=(0.0,0.0)
IER=0
IF(N) 10,100,100
10 IER=1
RETURN
100 IF(DREAL(Z))200,50,200
! Z=(0.0,0.0)
50 IF(N.EQ.0) BJ=(1.0,0.0)
RETURN
! CALCUL PAR LA METHODE DES SERIES ASCENDANTES
!
200 DO 55 KK=1,99
K=KK-1
XX=DFLOAT(N+K+1)
BJ= BJ+(((-0.25*Z**2)**K) /(AGAMMA(DFLOAT(KK))*AGAMMA(XX)))
IF(KK-10)55,9,90
90 AUX=(CDABS(BJ)-CDABS(BJ1))/CDABS(BJ1)
IF(DABS(AUX).LE.EPS) GO TO 25
9 BJ1=BJ
55 CONTINUE
! PRECISION NON ATTEINTE CODE ERREUR 3
IER=3
RETURN
25 BJ= ((0.5*Z)**N)*BJ
RETURN
END
In fact, this part is extracted from a long program, which consists in determining the BESSEL functions to study the dynamic behaviour of structures when they are subjected to a potential flow.
FUNCTION AGAMMA(XX)
! ---------------------
! COMPUTES THE GAMMA FUNCTION FOR A GIVEN ARGUMENT
! XX=THE ARGUMENT FOR THE GAMMA FUNCTION.
!
!
IMPLICIT REAL*8(A-H,O-Z)
REAL (KIND=8):: AGAMMA !ADDED BY FARHAD FOR QUAD PERCISION
IF(XX-100.0)6,6,4
4 IER=2
AGAMMA=9.33262E155_8 !IT IS MODIFIED BY FARHAD
RETURN
6 X=XX
ERR=1.0D-6
IER=0
AGAMMA=1.0
IF(X-2.0)50,50,15
10 IF(X-2.0)110,110,15
15 X=X-1.0
AGAMMA=AGAMMA*X
GO TO 10
50 IF(X-1.0)60,120,110
!
! SEE IF X IS NEAR NEGATIVE INTEGER OR ZERO
!
60 IF(X-ERR)62,62,80
62 Y=DFLOAT(INT(X))-X
IF(DABS(Y)-ERR)130,130,64
64 IF(1.0-Y-ERR)130,130,70
!
! X NOT NEAR A NEGATIVE INTEGER OR ZERO
!
70 IF(X-1.0)80,80,110
80 AGAMMA=AGAMMA/X
X=X+1.0
GO TO 70
110 Y=X-1.0
GY= 1.0+Y*(-0.5771017+Y*(0.9858540+Y*(-0.8764218+Y*(0.8328212+ &
Y*(-0.5684729+Y*(0.2548205+Y*(-0.05149930)))))))
AGAMMA=AGAMMA*GY
120 RETURN
130 IER=1
RETURN
END