Line data    Source code

       1             :       MODULE m_rcerf
       2             :       use m_juDFT
       3             : c*********************************************************************
       4             : c     calculates  real( erf(x+iy) ) for z=x+iy in the first quadrant.
       5             : c             m. weinert   may 1987
       6             : c     new declaration part
       7             : c                       s. bl"ugel, IFF, Nov.97
       8             : 
       9             :       PRIVATE
      10             :       PUBLIC rcerf
      11             :       CONTAINS
      12             : c*********************************************************************
      13      860045 :       REAL FUNCTION rcerf(x,y)
      14             :       IMPLICIT NONE
      15             : C     ..
      16             : C     .. Scalar Arguments ..
      17             :       REAL    x,y
      18             : C     ..
      19             : C     .. Local Scalars ..
      20             :       COMPLEX z
      21             :       REAL wr,wi
      22             : c
      23      860045 :       z = cmplx(x,y)
      24             : c
      25             : c--->    calculate w(z')=exp(-z'**2)*erfc(-iz') for -iz'=z
      26             : c
      27      860045 :       CALL wofz(y,x,wr,wi)
      28             : c
      29             : c--->    erf(x+iy)=1-exp(-z**2)*conjg( w(y+ix) )
      30             : c
      31      860045 :       RCERF = 1.0 - real( exp( -z*z ) * cmplx( wr,-wi ) )
      32             : 
      33      860045 :       END FUNCTION rcerf
      34             : 
      35             : c*********************************************************************
      36      860045 :       SUBROUTINE wofz(x,y,wr,wi)
      37             : c     calculates  w(z) = exp(-z**2) erfc(-iz)   for z = x + iy in the
      38             : c     first quadrant of the complex plane. based on acm algorithm 363
      39             : c     by w. gautschi, comm. acm 12, 635 (1969). the accuracy is about
      40             : c     10**-10.
      41             : c                       m. weinert    may 1987
      42             : c     new declaration part
      43             : c                       s. bl"ugel, IFF, Nov.97
      44             : c*********************************************************************
      45             :       USE m_constants
      46             :       IMPLICIT NONE
      47             : C     ..
      48             : C     .. Scalar Arguments ..
      49             :       REAL    x,y,wr,wi
      50             : C     ..
      51             : C     .. Local Scalars ..
      52             :       INTEGER icap,nu,n
      53             :       REAL    c,h,h2,plam,r1,r2,s,s1,s2,tsqpi,t1,t2
      54             :       LOGICAL bol
      55             : 
      56             : c--->    2/sqrt(pi) required for the normalization
      57      860045 :       tsqpi = 2.0/sqrt(pi_const) 
      58             : c
      59             : c--->    stop if not in first quadrant 
      60      860045 :       IF ( x<0.0 .OR. y<0 )  CALL juDFT_error("wofz",calledby ="rcerf")
      61             : c
      62      860045 :       IF ( y.LT.4.29 .AND. x.LT.5.33 ) THEN
      63      525773 :         s=(1.-y/4.29)*sqrt(1.-x*x/28.41)
      64      525773 :         h=1.6*s
      65      525773 :         h2=2*h
      66      525773 :         icap=6+23*s
      67      525773 :         nu=9+21*s
      68      525773 :         plam=h2**icap
      69      525773 :         bol=.true.
      70             :       ELSE
      71             :         h=0.0
      72             :         icap=0
      73             :         nu=8
      74             :         bol=.false.
      75             :       END IF
      76             : c
      77      860045 :       r1=0.0
      78      860045 :       r2=0.0
      79      860045 :       s1=0.0
      80      860045 :       s2=0.0
      81    12545438 :       DO n = nu,0,-1
      82    11685393 :          t1=y+h+(n+1)*r1
      83    11685393 :          t2=x-(n+1)*r2
      84    11685393 :          c=0.5/(t1*t1+t2*t2)
      85    11685393 :          r1=c*t1
      86    11685393 :          r2=c*t2
      87    12545438 :          IF ( bol .AND. n.LE.icap ) then
      88     7469513 :            t1=plam+s1
      89     7469513 :            s1=r1*t1-r2*s2
      90     7469513 :            s2=r2*t1+r1*s2
      91     7469513 :            plam=plam/h2
      92             :          END IF
      93             :       ENDDO
      94      860045 :       IF ( bol ) THEN
      95      525773 :         wr=tsqpi*s1
      96      525773 :         wi=tsqpi*s2
      97             :       ELSE
      98      334272 :         wr=tsqpi*r1
      99      334272 :         wi=tsqpi*r2
     100             :       END IF
     101      860045 :       IF ( y.EQ.0.0 ) wr=exp(-x*x)
     102             : 
     103      860045 :       END SUBROUTINE wofz
     104             : 
     105             :       END MODULE m_rcerf