Complementary Error Function

c **********************************************************************
c     derfc.F used with permission from Naval Surface Warfare Group
!******************************************************************************
      function derfc (x)
!
!******************************************************************************
!
!! DERFC: the complementary error function
!
      double precision derfc
      double precision an, ax, c, eps, rpinv, t, x, w
      double precision a(21), b(44), e(44)
      double precision dpmpar, dcsevl
!
!     rpinv = 1/sqrt(pi)
!
  data rpinv /.56418958354775628694807945156077259d0/
!
  data a(1)  / .1283791670955125738961589031215d+00/, &
       a(2)  /-.3761263890318375246320529677070d+00/, &
       a(3)  / .1128379167095512573896158902931d+00/, &
       a(4)  /-.2686617064513125175943235372542d-01/, &
       a(5)  / .5223977625442187842111812447877d-02/, &
       a(6)  /-.8548327023450852832540164081187d-03/, &
       a(7)  / .1205533298178966425020717182498d-03/, &
       a(8)  /-.1492565035840625090430728526820d-04/, &
       a(9)  / .1646211436588924261080723578109d-05/, &
       a(10) /-.1636584469123468757408968429674d-06/
  data a(11) / .1480719281587021715400818627811d-07/, &
       a(12) /-.1229055530145120140800510155331d-08/, &
       a(13) / .9422759058437197017313055084212d-10/, &
       a(14) /-.6711366740969385085896257227159d-11/, &
       a(15) / .4463222608295664017461758843550d-12/, &
       a(16) /-.2783497395542995487275065856998d-13/, &
       a(17) / .1634095572365337143933023780777d-14/, &
       a(18) /-.9052845786901123985710019387938d-16/, &
       a(19) / .4708274559689744439341671426731d-17/, &
       a(20) /-.2187159356685015949749948252160d-18/, &
       a(21) / .7043407712019701609635599701333d-20/
!
  data b(1)  / .610143081923200417926465815756d+00/, &
       b(2)  /-.434841272712577471828182820888d+00/, &
       b(3)  / .176351193643605501125840298123d+00/, &
       b(4)  /-.607107956092494148600512158250d-01/, &
       b(5)  / .177120689956941144861471411910d-01/, &
       b(6)  /-.432111938556729381859986496800d-02/, &
       b(7)  / .854216676887098678819832055000d-03/, &
       b(8)  /-.127155090609162742628893940000d-03/, &
       b(9)  / .112481672436711894688470720000d-04/, &
       b(10) / .313063885421820972630152000000d-06/
  data b(11) /-.270988068537762022009086000000d-06/, &
       b(12) / .307376227014076884409590000000d-07/, &
       b(13) / .251562038481762293731400000000d-08/, &
       b(14) /-.102892992132031912759000000000d-08/, &
       b(15) / .299440521199499393630000000000d-10/, &
       b(16) / .260517896872669362900000000000d-10/, &
       b(17) /-.263483992417196938600000000000d-11/, &
       b(18) /-.643404509890636443000000000000d-12/, &
       b(19) / .112457401801663447000000000000d-12/, &
       b(20) / .172815333899860980000000000000d-13/
  data b(21) /-.426410169494237500000000000000d-14/, &
       b(22) /-.545371977880191000000000000000d-15/, &
       b(23) / .158697607761671000000000000000d-15/, &
       b(24) / .208998378443340000000000000000d-16/, &
       b(25) /-.590052686940900000000000000000d-17/, &
       b(26) /-.941893387554000000000000000000d-18/, &
       b(27) / .214977356470000000000000000000d-18/, &
       b(28) / .466609850080000000000000000000d-19/, &
       b(29) /-.724301186200000000000000000000d-20/, &
       b(30) /-.238796682400000000000000000000d-20/
  data b(31) / .191177535000000000000000000000d-21/, &
       b(32) / .120482568000000000000000000000d-21/, &
       b(33) /-.672377000000000000000000000000d-24/, &
       b(34) /-.574799700000000000000000000000d-23/, &
       b(35) /-.428493000000000000000000000000d-24/, &
       b(36) / .244856000000000000000000000000d-24/, &
       b(37) / .437930000000000000000000000000d-25/, &
       b(38) /-.815100000000000000000000000000d-26/, &
       b(39) /-.308900000000000000000000000000d-26/, &
       b(40) / .930000000000000000000000000000d-28/
  data b(41) / .174000000000000000000000000000d-27/, &
       b(42) / .160000000000000000000000000000d-28/, &
       b(43) /-.800000000000000000000000000000d-29/, &
       b(44) /-.200000000000000000000000000000d-29/
!
  data e(1)  / .107797785207238315116833591035d+01/, &
       e(2)  /-.265598904091486733721465009040d-01/, &
       e(3)  /-.148707314669809950960504633300d-02/, &
       e(4)  /-.138040145414143859607708920000d-03/, &
       e(5)  /-.112803033322874914985073660000d-04/, &
       e(6)  /-.117286984274372522405373900000d-05/, &
       e(7)  /-.103476150393304615537382000000d-06/, &
       e(8)  /-.118991140858924382544470000000d-07/, &
       e(9)  /-.101622254498949864047600000000d-08/, &
       e(10) /-.137895716146965692169000000000d-09/
  data e(11) /-.936961303373730333500000000000d-11/, &
       e(12) /-.191880958395952534900000000000d-11/, &
       e(13) /-.375730172019937070000000000000d-13/, &
       e(14) /-.370537260269833570000000000000d-13/, &
       e(15) / .262756542349037100000000000000d-14/, &
       e(16) /-.112132287643793300000000000000d-14/, &
       e(17) / .184136028922538000000000000000d-15/, &
       e(18) /-.491302565748860000000000000000d-16/, &
       e(19) / .107044551673730000000000000000d-16/, &
       e(20) /-.267189366240500000000000000000d-17/
  data e(21) / .649326867976000000000000000000d-18/, &
       e(22) /-.165399353183000000000000000000d-18/, &
       e(23) / .426056266040000000000000000000d-19/, &
       e(24) /-.112558407650000000000000000000d-19/, &
       e(25) / .302561744800000000000000000000d-20/, &
       e(26) /-.829042146000000000000000000000d-21/, &
       e(27) / .231049558000000000000000000000d-21/, &
       e(28) /-.654695110000000000000000000000d-22/, &
       e(29) / .188423140000000000000000000000d-22/, &
       e(30) /-.550434100000000000000000000000d-23/
  data e(31) / .163095000000000000000000000000d-23/, &
       e(32) /-.489860000000000000000000000000d-24/, &
       e(33) / .149054000000000000000000000000d-24/, &
       e(34) /-.459220000000000000000000000000d-25/, &
       e(35) / .143180000000000000000000000000d-25/, &
       e(36) /-.451600000000000000000000000000d-26/, &
       e(37) / .144000000000000000000000000000d-26/, &
       e(38) /-.464000000000000000000000000000d-27/, &
       e(39) / .151000000000000000000000000000d-27/, &
       e(40) /-.500000000000000000000000000000d-28/
  data e(41) / .170000000000000000000000000000d-28/, &
       e(42) /-.600000000000000000000000000000d-29/, &
       e(43) / .200000000000000000000000000000d-29/, &
       e(44) /-.100000000000000000000000000000d-29/
!
  eps = epsilon ( eps )
!
!                     dabs(x) <= 1
!
  ax = dabs(x)
  if (ax > 1.d0) go to 20
  t = x*x
  w = a(21)
  do 10 i = 1,20
     k = 21 - i
     w = t*w + a(k)
   10 continue
  derfc = 0.5d0 + (0.5d0 - x*(1.d0 + w))
  return
!
!                  1 < dabs(x)  <  2
!
   20 if (ax >= 2.d0) go to 30
  n = 44
  if (eps >= 1.d-20) n = 30
  t = (ax - 3.75d0)/(ax + 3.75d0)
  derfc = dcsevl(t, b, n)
   21 derfc = dexp(-x*x) * derfc
  if (x < 0.d0) derfc = 2.d0 - derfc
  return
!
!                  2 < dabs(x)  <  12
!
   30 if (x < -9.d0) go to 60
  if (x >= 12.d0) go to 40
  n = 44
  if (eps >= 1.d-20) n = 25
  t = (1.d0/x)**2
  w = (10.5d0*t - 1.d0)/(2.5d0*t + 1.d0)
  derfc = dcsevl(w, e, n) / ax
  go to 21
!
!                       x >= 12
!
   40 if (x > 50.d0) go to 70
  t = (1.d0/x)**2
  an = -0.5d0
  c  =  0.5d0
  w  =  0.0d0
   50    c = c + 1.d0
     an = - c*an*t
     w = w + an
     if (dabs(an) > eps) go to 50
  w = (-0.5d0 + w)*t + 1.d0
  derfc = dexp(-x*x) * ((rpinv*w)/ax)
  return
!
!             limit value for large negative x
!
   60 derfc = 2.d0
  return
!
!             limit value for large positive x
!
   70 derfc = 0.d0
  return
end

function dcsevl (x, a, n)
!
!*******************************************************************************
!
!! DCSEVL: evaluate the n term chebyshev series a at x.
!          only half of the first coefficient is used.
!
  double precision dcsevl
  double precision a(n),x,x2,s0,s1,s2
!
  if (n > 1) go to 10
     dcsevl = 0.5d0 * a(1)
     return
!
   10 x2 = x + x
  s0 = a(n)
  s1 = 0.d0
  do 20 i = 2,n
     s2 = s1
     s1 = s0
     k = n - i + 1
     s0 = x2*s1 - s2 + a(k)
   20 continue
  dcsevl = 0.5d0 * (s0 - s2)
  return
end