2016-07-13 3 views
1

私は、球面ベッセル関数を計算するために与えられたサブルーチンを使用するプログラムを書いています。私は、テーブルに1つの値しか与えない関数を与えるサブルーチンを修正しました。しかし、私は自分の関数を呼び出すときに私のプログラムにIMPLICIT DOUBLE PRECISION (A-H,O-Z)を持たなければならない、あるいはそれが私に間違った価値やエラーを与えることに気付きました。以下に、正しく動作するサンプルプログラムを含めました。Implicitはプログラムを不正確にします

! 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 

ここでは、同じ機能を使用して機能しないサンプルプログラムを示します。

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 

私は比較的新しいFortranを使用しているため、自分のプログラムが動作しない理由を理解できません。誰もが提供できる助けに感謝します。

+0

動作しないサンプルプログラムは、6.07E-02とは対照的に1.00を返します。 –

+0

プログラムのどの部分でも暗黙のうちに必要なのでしょうか?暗黙の型を使うのは良い習慣ではありません。すべての変数に明示的な型を与えて、どこにでも暗黙的にnoneを入れます。 – innoSPG

+0

マーク:サンプルプログラムは、上のプログラムで使用されたベッセルの実装を使用します。 –

答えて

2

私は、動作していないサンプルプログラムが、同じサンプル実装であるbesselを使用していると思います。その場合

は、仕様

implicit doubleprecision(A-Z) 
当たり全て倍精度である besselNXinteger非作動メインプログラムにおける)と対応する引数の型との間の種類の競合があります

ベッセル内のすべてがデフォルトでdoubleprecisionです。メインプログラムではNXdoubleprecisionと定義する必要があります。

上記のコメントで私が言ったように、最良の解決策はどこにでも明示的な入力を使用することです。

+1

ありがとうございます。私はXをdoubleprecisionとして明示的に定義しましたが、これは問題を修正したようです。これで正しい値が返されるようになりました。どうもありがとうございました。 –

+0

問題が解決する場合は、回答を受け入れて、今後同様の問題が発生した場合に役立ちます。 – innoSPG

関連する問題