Implizite keine macht Programm ungenau

Question

Ich schreibe ein Programm, das eine gegebene Subroutine verwendet, um sphärische Bessel-Funktionen zu berechnen. Ich habe das Unterprogramm modifiziert, das eine Tabelle in eine Funktion gibt, die nur einen Wert gibt. Ich erkannte jedoch, dass wenn ich meine Funktion anrufe, ich IMPLICIT DOUBLE PRECISION (A-H,O-Z) in meinem Programm haben muss, oder es wird mir einen falschen Wert oder einen Fehler geben. Unten habe ich ein Beispielprogramm eingefügt, das korrekt funktioniert.

! n = 3, x = 2 should return ~ 6.07E-2 
program hello 
IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
    doubleprecision :: bessel, ans 
    WRITE(*,*)'Please enter n and x ' 
    READ(*,*)N,X 

    ans = bessel(N,X) 

    print *, ans 
end program 

    SUBROUTINE SPHJ(N,X,NM,SJ,DJ) 

!  ======================================================= 
!  Purpose: Compute spherical Bessel functions jn(x) and 
!    their derivatives 
!  Input : x --- Argument of jn(x) 
!    n --- Order of jn(x) (n = 0,1,???) 
!  Output: SJ(n) --- jn(x) 
!    DJ(n) --- jn'(x) 
!    NM --- Highest order computed 
!  Routines called: 
!    MSTA1 and MSTA2 for computing the starting 
!    point for backward recurrence 
!  ======================================================= 

     IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
     DIMENSION SJ(0:N),DJ(0:N) 
     NM=N 
     IF (DABS(X).EQ.1.0D-100) THEN 
      DO 10 K=0,N 
       SJ(K)=0.0D0 
10   DJ(K)=0.0D0 
      SJ(0)=1.0D0 
      DJ(1)=.3333333333333333D0 
      RETURN 
     ENDIF 
     SJ(0)=DSIN(X)/X 
     SJ(1)=(SJ(0)-DCOS(X))/X 
     IF (N.GE.2) THEN 
      SA=SJ(0) 
      SB=SJ(1) 
      M=MSTA1(X,200) 
      IF (M.LT.N) THEN 
       NM=M 
      ELSE 
       M=MSTA2(X,N,15) 
      ENDIF 
      F0=0.0D0 
      F1=1.0D0-100 
      DO 15 K=M,0,-1 
       F=(2.0D0*K+3.0D0)*F1/X-F0 
       IF (K.LE.NM) SJ(K)=F 
       F0=F1 
15   F1=F 
      IF (DABS(SA).GT.DABS(SB)) CS=SA/F 
      IF (DABS(SA).LE.DABS(SB)) CS=SB/F0 
      DO 20 K=0,NM 
20   SJ(K)=CS*SJ(K) 
     ENDIF 
     DJ(0)=(DCOS(X)-DSIN(X)/X)/X 
     DO 25 K=1,NM 
25   DJ(K)=SJ(K-1)-(K+1.0D0)*SJ(K)/X 
     RETURN 
     END 


     INTEGER FUNCTION MSTA1(X,MP) 

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

     IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
     A0=DABS(X) 
     N0=INT(1.1*A0)+1 
     F0=ENVJ(N0,A0)-MP 
     N1=N0+5 
     F1=ENVJ(N1,A0)-MP 
     DO 10 IT=1,20 
      NN=N1-(N1-N0)/(1.0D0-F0/F1) 
      F=ENVJ(NN,A0)-MP 
      IF(ABS(NN-N1).LT.1) GO TO 20 
      N0=N1 
      F0=F1 
      N1=NN 
10  F1=F 
20  MSTA1=NN 
     RETURN 
     END 


     INTEGER FUNCTION MSTA2(X,N,MP) 

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

     IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
     A0=DABS(X) 
     HMP=0.5D0*MP 
     EJN=ENVJ(N,A0) 
     IF (EJN.LE.HMP) THEN 
      OBJ=MP 
      N0=INT(1.1*A0) 
     ELSE 
      OBJ=HMP+EJN 
      N0=N 
     ENDIF 
     F0=ENVJ(N0,A0)-OBJ 
     N1=N0+5 
     F1=ENVJ(N1,A0)-OBJ 
     DO 10 IT=1,20 
      NN=N1-(N1-N0)/(1.0D0-F0/F1) 
      F=ENVJ(NN,A0)-OBJ 
      IF (ABS(NN-N1).LT.1) GO TO 20 
      N0=N1 
      F0=F1 
      N1=NN 
10   F1=F 
20  MSTA2=NN+10 
     RETURN 
     END 

     REAL*8 FUNCTION ENVJ(N,X) 
     DOUBLE PRECISION X 
     ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N) 
     RETURN 
     END 

!end of file msphj.f90 

doubleprecision function bessel(N,X) 
implicit doubleprecision(A-Z) 
    DIMENSION SJ(0:250),DJ(0:250) 
    integer :: N 

    CALL SPHJ(N,X,N,SJ,DJ) 

    bessel = SJ(N) 

end function 

Und hier ist ein Beispielprogramm, das nicht funktioniert, mit der gleichen Funktion.

program hello 
    IMPLICIT none 
    doubleprecision :: bessel, ans 
    integer :: N, X 
    WRITE(*,*)'Please enter n and x ' 
    READ(*,*)N,X 

    ans = bessel(N,X) 

    print *, ans 

end program 

Ich bin relativ neu in Fortran und verstehe nicht, warum mein Programm nicht funktioniert. Ich schätze jede Hilfe, die jemand zur Verfügung stellen kann.

Answer 1

Ich denke, das nicht funktionierende Beispielprogramm verwendet die gleiche Implementierung von bessel als Arbeitsprobe.

Wenn ja, gibt es einen Konflikt vom Typ zwischen der Art der N und X (integer im Nicht-Haupt-Programm arbeiten) und die entsprechenden Argumente in bessel die alle mit doppelter Genauigkeit gemäß der Spezifikation sind

implicit doubleprecision(A-Z) 

Alles in Bessel ist standardmäßig doubleprecision. Ihr Hauptprogramm muss N und X als doubleprecision definieren.

Die beste Lösung, wie ich im obigen Kommentar gesagt habe, ist explizite Typisierung überall zu verwenden.