MODULE double_complex_erf ! Double precision complex error function ! Adapted from the Naval Surface Warfare Center Mathematics Library ! by Alan.Miller @ vic.cmis.csiro.au ! http://www.ozemail.com.au/~milleraj ! Latest revision - 18 September 2000 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) CONTAINS FUNCTION dnorm(x, y) RESULT(fn_val) ! Replaces the statement function anorm in the F77 code. REAL (dp), INTENT(IN) :: x, y REAL (dp) :: fn_val fn_val = MAX(ABS(x), ABS(y)) RETURN END FUNCTION dnorm SUBROUTINE dcrec (x, y, u, v) !----------------------------------------------------------------------- ! COMPLEX RECIPROCAL U + I*V = 1/(X + I*Y) !----------------------------------------------------------------------- REAL (dp), INTENT(IN) :: x, y REAL (dp), INTENT(OUT) :: u, v REAL (dp) :: d, t IF (ABS(x) > ABS(y)) GO TO 10 t = x/y d = y + t*x u = t/d v = -1._dp/d RETURN 10 t = y/x d = x + t*y u = 1._dp/d v = -t/d RETURN END SUBROUTINE dcrec SUBROUTINE cdivid (ar, ai, br, bi, cr, ci) !----------------------------------------------------------------------- ! REAL (dp) COMPLEX DIVISION C = A/B AVOIDING OVERFLOW !----------------------------------------------------------------------- REAL (dp), INTENT(IN) :: ar, ai, br, bi REAL (dp), INTENT(OUT) :: cr, ci REAL (dp) :: d, t, u, v IF (ABS(br) <= ABS(bi)) GO TO 10 t = bi/br d = br + t*bi u = (ar + ai*t)/d v = (ai - ar*t)/d cr = u ci = v RETURN 10 IF (bi == 0._dp) GO TO 20 t = br/bi d = bi + t*br u = (ar*t + ai)/d v = (ai*t - ar)/d cr = u ci = v RETURN ! DIVISION BY ZERO. C = INFINITY 20 cr = HUGE(1.0_dp) ci = cr RETURN END SUBROUTINE cdivid SUBROUTINE dcerf (mo, z, w) !----------------------------------------------------------------------- ! COMPUTATION OF THE COMPLEX ERROR FUNCTION ! ----------------- ! W = ERF(Z) IF MO = 0 ! W = ERFC(Z) OTHERWISE !----------------------------------------------------------------------- INTEGER, INTENT(IN) :: mo REAL (dp), INTENT(IN) :: z(2) REAL (dp), INTENT(OUT) :: w(2) REAL (dp) :: m, n, n2, n4, np1 REAL (dp) :: c2, d, d2, e, eps, r, sn, tol, x, y REAL (dp) :: a0(2), an(2), b0(2), bn(2) REAL (dp) :: g0(2), gn(2), h0(2), hn(2) REAL (dp) :: qf(2), sm(2), sz(2), tm(2), ts(2), w0(2), wn(2) !------------------------ ! C = 1/SQRT(PI) !------------------------ REAL (dp), PARAMETER :: c = .56418958354775628694807945156077_dp !------------------------ ! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE ! SMALLEST NUMBER SUCH THAT 1._dp + EPS > 1._dp. eps = EPSILON(1.0_dp) !------------------------ x = z(1) y = z(2) sn = 1._dp IF (x >= 0._dp) GO TO 10 x = -x y = -y sn = -1._dp 10 r = x*x + y*y sz(1) = x*x - y*y sz(2) = 2._dp*x*y IF (r <= 1._dp) GO TO 20 IF (r >= 144._dp) GO TO 100 IF (ABS(y) > 31.8_dp*x) GO TO 50 IF (ABS(y) > 7.0_dp*x .AND. r < 64._dp) GO TO 50 IF (ABS(y) > 3.2_dp*x .AND. r < 49._dp) GO TO 50 IF (ABS(y) > 2.0_dp*x .AND. r < 36._dp) GO TO 50 IF (ABS(y) > 1.2_dp*x .AND. r < 25._dp) GO TO 50 IF (ABS(y) > 0.9_dp*x .AND. r < 16._dp) GO TO 50 IF (r >= 6.25_dp) GO TO 80 IF (ABS(y) > 0.6_dp*x) GO TO 50 IF (r >= 4.0_dp) GO TO 40 d = x - 2._dp IF (d*d + y*y < 1._dp) GO TO 40 GO TO 50 ! TAYLOR SERIES 20 c2 = c + c tm(1) = c2*x tm(2) = c2*y sm(1) = tm(1) sm(2) = tm(2) tol = 2._dp*eps m = 0._dp 21 m = m + 1._dp d = m + m + 1._dp ts(1) = tm(1)*sz(1) - tm(2)*sz(2) ts(2) = tm(1)*sz(2) + tm(2)*sz(1) tm(1) = -ts(1)/m tm(2) = -ts(2)/m ts(1) = tm(1)/d ts(2) = tm(2)/d sm(1) = sm(1) + ts(1) sm(2) = sm(2) + ts(2) IF (dnorm(ts(1),ts(2)) > tol*dnorm(sm(1),sm(2))) GO TO 21 IF (mo /= 0) GO TO 30 w(1) = sn*sm(1) w(2) = sn*sm(2) RETURN 30 IF (sn == 1._dp) GO TO 31 w(1) = 1._dp + sm(1) w(2) = sm(2) RETURN 31 w(1) = 0.5_dp + (0.5_dp - sm(1)) w(2) = -sm(2) RETURN ! TAYLOR SERIES AROUND Z0 = 2 40 tm(1) = x tm(2) = y CALL erfcm2 (0, tm, w) IF (mo /= 0) GO TO 41 w(1) = sn*(0.5_dp + (0.5_dp - w(1))) w(2) = - sn*w(2) RETURN 41 IF (sn > 0._dp) RETURN w(1) = 2._dp - w(1) w(2) = - w(2) RETURN ! PADE APPROXIMATION FOR THE TAYLOR SERIES ! FOR (EXP(Z*Z)/Z)*ERF(Z) 50 d = 4._dp IF (r > 16._dp) d = 16._dp IF (r > 64._dp) d = 64._dp d2 = d*d CALL dcrec (sz(1), sz(2), w(1), w(2)) a0(1) = 1._dp a0(2) = 0._dp an(1) = (w(1) + 4._dp/15._dp)*d an(2) = w(2)*d b0(1) = 1._dp b0(2) = 0._dp bn(1) = (w(1) - 0.4_dp)*d bn(2) = w(2)*d CALL cdivid (an(1), an(2), bn(1), bn(2), wn(1), wn(2)) tol = 10._dp*eps n4 = 0._dp 60 n4 = n4 + 4._dp e = (n4 + 1._dp)*(n4 + 5._dp) tm(1) = d*(w(1) - 2._dp/e) tm(2) = d*w(2) e = d2*(n4*(n4 + 2.0))/((n4 - 1.0)*(n4 + 3.0)*(n4 + 1.0)**2) qf(1) = (tm(1)*an(1) - tm(2)*an(2)) + e*a0(1) qf(2) = (tm(1)*an(2) + tm(2)*an(1)) + e*a0(2) a0(1) = an(1) a0(2) = an(2) an(1) = qf(1) an(2) = qf(2) qf(1) = (tm(1)*bn(1) - tm(2)*bn(2)) + e*b0(1) qf(2) = (tm(1)*bn(2) + tm(2)*bn(1)) + e*b0(2) b0(1) = bn(1) b0(2) = bn(2) bn(1) = qf(1) bn(2) = qf(2) w0(1) = wn(1) w0(2) = wn(2) CALL cdivid (an(1), an(2), bn(1), bn(2), wn(1), wn(2)) IF (dnorm(wn(1) - w0(1), wn(2) - w0(2)) > tol*dnorm(wn(1), wn(2))) GO TO 60 c2 = c + c sm(1) = c2*(x*wn(1) - y*wn(2)) sm(2) = c2*(x*wn(2) + y*wn(1)) e = EXP(-sz(1)) qf(1) = e*COS(-sz(2)) qf(2) = e*SIN(-sz(2)) tm(1) = qf(1)*sm(1) - qf(2)*sm(2) tm(2) = qf(1)*sm(2) + qf(2)*sm(1) w(1) = sn*tm(1) w(2) = sn*tm(2) IF (mo == 0) RETURN w(1) = 1._dp - w(1) w(2) = - w(2) RETURN ! PADE APPROXIMATION FOR THE ASYMPTOTIC EXPANSION ! FOR Z*EXP(Z*Z)*ERFC(Z) 80 d = 4._dp*r IF (r < 16._dp) d = 16._dp*r d2 = d*d tm(1) = sz(1) + sz(1) tm(2) = sz(2) + sz(2) g0(1) = 1._dp g0(2) = 0._dp gn(1) = (2._dp + tm(1))/d gn(2) = tm(2)/d h0(1) = 1._dp h0(2) = 0._dp tm(1) = 3._dp + tm(1) hn(1) = tm(1)/d hn(2) = tm(2)/d CALL cdivid (gn(1), gn(2), hn(1), hn(2), wn(1), wn(2)) np1 = 1._dp tol = 10._dp*eps 90 n = np1 np1 = n + 1._dp n2 = n + n e = (n2*(n2 + 1._dp))/d2 tm(1) = tm(1) + 4._dp qf(1) = (tm(1)*gn(1) - tm(2)*gn(2))/d - e*g0(1) qf(2) = (tm(1)*gn(2) + tm(2)*gn(1))/d - e*g0(2) g0(1) = gn(1) g0(2) = gn(2) gn(1) = qf(1) gn(2) = qf(2) qf(1) = (tm(1)*hn(1) - tm(2)*hn(2))/d - e*h0(1) qf(2) = (tm(1)*hn(2) + tm(2)*hn(1))/d - e*h0(2) h0(1) = hn(1) h0(2) = hn(2) hn(1) = qf(1) hn(2) = qf(2) w0(1) = wn(1) w0(2) = wn(2) CALL cdivid (gn(1), gn(2), hn(1), hn(2), wn(1), wn(2)) IF (dnorm(wn(1) - w0(1), wn(2) - w0(2)) > tol*dnorm(wn(1), wn(2))) GO TO 90 tm(1) = x*hn(1) - y*hn(2) tm(2) = x*hn(2) + y*hn(1) CALL cdivid (c*gn(1), c*gn(2), tm(1), tm(2), sm(1), sm(2)) GO TO 130 ! ASYMPTOTIC EXPANSION 100 CALL dcrec (x, y, tm(1), tm(2)) sm(1) = tm(1) sm(2) = tm(2) qf(1) = tm(1)*tm(1) - tm(2)*tm(2) qf(2) = 2._dp*tm(1)*tm(2) tol = 2._dp*eps d = -0.5_dp 110 d = d + 1._dp ts(1) = tm(1)*qf(1) - tm(2)*qf(2) ts(2) = tm(1)*qf(2) + tm(2)*qf(1) tm(1) = -d*ts(1) tm(2) = -d*ts(2) sm(1) = sm(1) + tm(1) sm(2) = sm(2) + tm(2) IF (dnorm(tm(1), tm(2)) > tol*dnorm(sm(1), sm(2))) GO TO 110 sm(1) = c*sm(1) sm(2) = c*sm(2) IF (x < 1.d-2) GO TO 200 ! TERMINATION 130 e = EXP(-sz(1)) qf(1) = e*COS(-sz(2)) qf(2) = e*SIN(-sz(2)) ts(1) = qf(1)*sm(1) - qf(2)*sm(2) ts(2) = qf(1)*sm(2) + qf(2)*sm(1) sm(1) = ts(1) sm(2) = ts(2) IF (mo /= 0) GO TO 140 w(1) = sn*(0.5_dp + (0.5_dp - sm(1))) w(2) = - sn*sm(2) RETURN 140 IF (sn == 1._dp) GO TO 141 w(1) = 2._dp - sm(1) w(2) = -sm(2) RETURN 141 w(1) = sm(1) w(2) = sm(2) RETURN ! MODIFIED ASYMPTOTIC EXPANSION 200 e = EXP(-sz(1)) qf(1) = e*COS(-sz(2)) qf(2) = e*SIN(-sz(2)) w(1) = qf(1)*sm(1) - qf(2)*sm(2) w(2) = qf(1)*sm(2) + qf(2)*sm(1) IF (mo == 0) GO TO 210 w(1) = 1._dp + sn*w(1) w(2) = sn*w(2) RETURN 210 IF (sn < 0.0) RETURN w(1) = - w(1) w(2) = - w(2) RETURN END SUBROUTINE dcerf SUBROUTINE erfcm2 (mo, z, w) !----------------------------------------------------------------------- ! CALCULATION OF ERFC(Z) USING THE TAYLOR SERIES ! AROUND Z0 = 2 !----------------------------------------------------------------------- INTEGER, INTENT(IN) :: mo REAL (dp), INTENT(IN) :: z(2) REAL (dp), INTENT(OUT) :: w(2) REAL (dp) :: eps, h(2), t(2), tol, x, y INTEGER :: j, n !------------------------------ ! C = (2/SQRT(PI))*EXP(-4) ! E = ERFC(2) !------------------------------ REAL (dp), PARAMETER :: c = .020666985354092053857068941306585476_dp, & e = .0046777349810472658379307436327470714_dp !------------------------------ REAL (dp) :: a(63) = (/ .20000000000000000000000000000000000D+01, & .23333333333333333333333333333333333D+01, .16666666666666666666666666666666667D+01, & .63333333333333333333333333333333333D+00, -.22222222222222222222222222222222222D-01, & -.16349206349206349206349206349206349D+00, -.76984126984126984126984126984126984D-01, & -.24250440917107583774250440917107584D-02, .12716049382716049382716049382716049D-01, & .50208433541766875100208433541766875D-02, -.25305969750414194858639303083747528D-03, & -.78593217482106370995259884148773038D-03, -.19118154038788959423880058800693721D-03, & .46324144207742091339974937858535742D-04, .33885549097189308829520469732109944D-04, & .28637897646612243562134629672756034D-05, -.29071891082127275370004560446169188D-05, & -.89674405786490646425523560263096103D-06, .96069103941908684338469767911200105D-07, & .99432863129093191401848891268744113D-07, .97610310501460621303387795457283579D-08, & -.65557500375673133822289344530892436D-08, -.18706782059105426900361744016236561D-08, & .20329898993447386223176373714372370D-09, .16941915827254374668448114614201210D-09, & .10619149520827430973786114446699534D-10, -.10136148256511788733365237088810952D-10, & -.21042890133669970575386166675721692D-11, .37186985840699828780916522245407087D-12, & .17921843632701679986488128324051002D-12, -.89823991804248069863542565948598397D-16, & -.10533182313660970970232171410372199D-13, -.12340742690978398320850088252659714D-14, & .44315624546581333350568244777175883D-15, .11584041639989442481950487524296214D-15, & -.10765703619385988116658460442219647D-16, -.70653158723054941879586082239984222D-17, & -.18708903154917138727191793341667090D-18, .32549879318817103966053527398133297D-18, & .40654116689599228385911733319215613D-19, -.11250074516817311101947327325293424D-19, & -.28923865378584966737386008432031980D-20, .23653053641701517160704870522922706D-21, & .14665384680061888088099002254334292D-21, .26971039707314316218154193225264469D-23, & -.58753834789274356433279671015522650D-23, -.59960357240498652932299485494869633D-24, & .18586826578121663981412155416486531D-24, .38364131854692721721867481914852428D-25, & -.41342210492630142578080062451711039D-26, -.17646283105274988992381528904600860D-26, & .19828685934364181151988692232131608D-28, .65592252170840353572672782446212733D-28, & .40626551379996340638338449938639730D-29, -.20097984104191034713653294173834095D-29, & -.28104226475997460044096389060743131D-30, .48705319298749358709127987806547949D-31, & .12664655832830787747161769929972617D-31, -.75168312488894341862391776330113688D-33, & -.45760473722605993842481669806804415D-33, -.56725491529575395930156379514718000D-35, & .13932664042920082608489441616061541D-34, .10452448992516358449586503951463322D-35 /) !------------------------------ ! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE ! SMALLEST NUMBER SUCH THAT 1._dp + EPS .GT. 1._dp eps = EPSILON(1.0_dp) !------------------------------ tol = eps*1.d+12 h(1) = 1._dp + (1._dp - z(1)) h(2) = - z(2) x = 1._dp y = 0._dp w(1) = a(30) w(2) = 0._dp DO n = 31,63 t(1) = x*h(1) - y*h(2) t(2) = x*h(2) + y*h(1) x = t(1) y = t(2) t(1) = a(n)*x t(2) = a(n)*y w(1) = w(1) + t(1) w(2) = w(2) + t(2) IF (dnorm(t(1), t(2)) <= tol*dnorm(w(1), w(2))) EXIT END DO DO j = 1,29 n = 30 - j x = h(1)*w(1) - h(2)*w(2) w(2) = h(1)*w(2) + h(2)*w(1) w(1) = a(n) + x END DO x = h(1)*w(1) - h(2)*w(2) w(2) = h(1)*w(2) + h(2)*w(1) w(1) = 1._dp + x x = c*(h(1)*w(1) - h(2)*w(2)) w(2) = c*(h(1)*w(2) + h(2)*w(1)) w(1) = e + x IF (mo == 0) RETURN ! COMPUTE EXP(Z*Z)*ERFC(Z) x = z(1)*z(1) - z(2)*z(2) y = 2._dp*z(1)*z(2) x = EXP(x) t(1) = x*COS(y) t(2) = x*SIN(y) x = t(1)*w(1) - t(2)*w(2) y = t(1)*w(2) + t(2)*w(1) w(1) = x w(2) = y RETURN END SUBROUTINE erfcm2 END MODULE double_complex_erf