C ---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      parameter(amin=0.0,amax=2.0,astep=0.1)

      a = amin
      do 1200 ic=1,1000
        if(a.gt.amax) goto 1250
        write(6,6100) a, erfc_(a)
        a = a + astep
 1200 continue
 1250 continue
      read(5,6200) dummy

 6100 format(' erfc(',f7.4,') = ',f9.6)
 6200 format(a)

      stop
      end

C =====================================================================
C ---------------------------------------------------------------------
      real*8 function erfc_(a)
C ---------------------------------------------------------------------
C  calculate erfc(a)
C  using formula 1 when a is small, and formula 2 when a is large
C ---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      parameter(alo=1.5,ahi=2.0)

      external fml1,fml2

      if(a.lt.0) then
        print *,'ERROR: illegal a=', a
        erfc_ = 0
      else if(a.eq.0) then
        erfc_ = 1.d0
      else if(a.le.alo) then
        ndiv = 1 + 5*a
        erfc_ = extpl(a,fml1,ndiv)
      else if(a.ge.ahi) then
        ndiv = 1 + 5.0/a
        erfc_ = extpl(a,fml2,ndiv)
      else
        ndiv = 1 + 5*a
        v1 = extpl(a,fml1,ndiv)
        ndiv = 1 + 5.0/a
        v2 = extpl(a,fml2,ndiv)
        erfc_ = v1 + (v2-v1)*(a-alo)/(ahi-alo)
      endif

      return
      end

C ---------------------------------------------------------------------
      real*8 function extpl(a,f,nd0)
C ---------------------------------------------------------------------
C  calculate f(a) by extrapolation,
C  i.e. as limiting value of f(a,nd) with nd (=#div) going to infinity
C  where f(a,nd) is intergated evaluation of f(a) with mesh dx
C ---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      parameter(MAXTRY=20,ERRMAX=1.e-10)
      real*8 x(0:MAXTRY), p(0:MAXTRY), q(0:MAXTRY), r(0:MAXTRY)

      nd = nd0
      d  = 1.d0/nd

      do 2100 i=0, MAXTRY-1

        fa = f(a,nd)
        x(i) = 1.0/nd
	p(i) = fa
	q(i) = fa
	r(i) = fa
        qq = q(i)

        do 2150 j=i-1, 0, -1
      	  h    = x(i) - x(j)
      	  dtmp = (r(j)-qq) / h
      	  qq   = dtmp * x(j)
          r(j) = dtmp * x(i)
          p(j) = p(j+1) + r(j)
 2150   continue

C       now p(0) approximates f(a,0) with err=r(0)

        extpl = p(0)
	ep    = abs(r(0))
	if(ep<ERRMAX) then
          return
        else
          nd = (3*nd+1)/2
        endif

 2100 continue

      print *,'warning: give up with err = ', ep
      return
      end

C ---------------------------------------------------------------------
      real*8 function fml1(a,nd)
C ---------------------------------------------------------------------
C  evaluate erfc(a)
C  by integrating [ 1/(1+y*y) - exp(-y*y) ] ( 0 < y < a )
C  this is good when a is not large
C ---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)

      PI = 4.d0*datan(1.d0)
      CN = 2.d0/sqrt(PI)

      sum = 0.d0

      do 3100 j=1, nd-1
        y = a*dble(j)/nd
        s = y*y
        sum = sum + 1.d0/(1+s) - exp(-s)
 3100 continue

      y = a
      s = y*y
      sum = sum + 0.5d0*( 1.d0/(1+s) - exp(-s) )

      fml1 =  1.d0 + CN*( (sum*a/nd) - atan(a) )
      return
      end

C ---------------------------------------------------------------------
      real*8 function fml2(a,nd)
C ---------------------------------------------------------------------
C  evaluate erfc(a)
C  by integrating [ exp(-1/y*y)/(y*y) ] ( 0 < y < 1/a )
C  this is good when a is large
C ---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)

      PI = 4.d0*datan(1.d0)
      CN = 2.d0/sqrt(PI)

      sum = 0.d0

      do 3200 j=1, nd-1
        y = (1.d0/a)*(dble(j)/nd)
        s = 1.d0/(y*y)
        sum = sum + s*exp(-s)
 3200 continue

      y = 1.d0/a
      s = 1.d0/(y*y)
      sum = sum + 0.5d0*s*exp(-s)

      fml2 =  CN*sum/(nd*a)
      return
      end

C ---------------------------------------------------------------------

