cerror Subroutine

subroutine cerror ( z, cer )

  !*****************************************************************************80
  !
  !! CERROR computes the error function for a complex argument.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    15 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  ! 
  !  Parameters:
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, complex ( kind = 8 ) CER, the function value.
  !
  implicit none

  real ( kind = 8 ) a0
  complex ( kind = 8 ) c0
  complex ( kind = 8 ) cer
  complex ( kind = 8 ) cl
  complex ( kind = 8 ) cr
  complex ( kind = 8 ) cs
  integer ( kind = 4 ) k
  real ( kind = 8 ) pi
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z1

  a0 = abs ( z )
  c0 = exp ( - z * z )
  pi = 3.141592653589793D+00
  z1 = z

  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
     z1 = - z
  end if

  if ( a0 <= 5.8D+00 ) then    

     cs = z1
     cr = z1
     do k = 1, 120
        cr = cr * z1 * z1 / ( k + 0.5D+00 )
        cs = cs + cr
        if ( abs ( cr / cs ) < 1.0D-15 ) then
           exit
        end if
     end do

     cer = 2.0D+00 * c0 * cs / sqrt ( pi )

  else

     cl = 1.0D+00 / z1              
     cr = cl
     do k = 1, 13
        cr = -cr * ( k - 0.5D+00 ) / ( z1 * z1 )
        cl = cl + cr
        if ( abs ( cr / cl ) < 1.0D-15 ) then
           exit
        end if
     end do

     cer = 1.0D+00 - c0 * cl / sqrt ( pi )

  end if

  if ( real ( z, kind = 8 ) < 0.0D+00 ) then
     cer = -cer
  end if

  return
end subroutine cerror