module math

  use const_and_precisions, only : wp_, zero, one
  implicit none

contains

  function catand(z)
!***begin prologue  catan
!***purpose  compute the complex arc tangent.
!***library   slatec (fnlib)
!***category  c4a
!***type      complex (catan-c)
!***keywords  arc tangent, elementary functions, fnlib, trigonometric
!***author  fullerton, w., (lanl)
!***description
!
! catan(z) calculates the complex trigonometric arc tangent of z.
! the result is in units of radians, and the real part is in the first
! or fourth quadrant.
!
!***references  (none)
!***routines called  (none)
!***revision history  (yymmdd)
!   770801  date written
!   890531  changed all specific intrinsics to generic.  (wrb)
!   890531  revision date from version 3.2
!   891214  prologue converted to version 4.0 format.  (bab)
!   900315  calls to xerror changed to calls to xermsg.  (thj)
!   900326  removed duplicate information from description section.
!       (wrb)
!***end prologue  catan
    use const_and_precisions, only : comp_eps, pi2=>pihalf, czero, cunit
    implicit none
    complex(wp_) :: catand
    complex(wp_), intent(in) :: z
    complex(wp_) :: z2
    real(wp_) :: r,x,y,r2,xans,yans,twoi
    integer :: i
    logical, save :: first=.true.
    integer, save :: nterms
    real(wp_), save :: rmin, rmax, sqeps
!***first executable statement  catan
    if (first) then
! nterms = log(eps)/log(rbnd) where rbnd = 0.1
      nterms = int(-0.4343_wp_*log(0.5_wp_*comp_eps) + 1.0_wp_)
      sqeps = sqrt(comp_eps)
      rmin = sqrt (1.5_wp_*comp_eps)
      rmax = 2.0_wp_/comp_eps
    endif
    first = .false.
!
    r = abs(z)
    if (r<=0.1_wp_) then
!
      catand = z
      if (r<rmin) return
!
      catand = czero
      z2 = z*z
      do i=1,nterms
        twoi = 2*(nterms-i) + 1
        catand = 1.0_wp_/twoi - z2*catand
      end do
      catand = z*catand
!
    else if (r<=rmax) then
      x = real(z)
      y = aimag(z)
      r2 = r*r
      if (r2==one.and.x==zero) print*,'catand, z is +i or -i'
      if (abs(r2-one)<=sqeps) then
        if (abs(cunit+z*z) < sqeps) &
         print*,'catand, answer lt half precision, z**2 close to -1'
!
      end if
      xans = 0.5_wp_*atan2(2.0_wp_*x, one-r2)
      yans = 0.25_wp_*log((r2+2.0_wp_*y+one)/(r2-2.0_wp_*y+one))
      catand = cmplx(xans, yans, wp_)
!
    else
      catand = cmplx(pi2, zero, wp_)
      if (real(z)<zero) catand = cmplx(-pi2, zero, wp_)
    end if
  end function catand


  function fact(k)
    ! Factorial function k!
    ! Note: the result is real

    implicit none
    integer, intent(in) :: k
    real(wp_)           :: fact

    fact = gamma(real(1 + k, kind=wp_))

   end function fact

end module math