FIGINI LORENZO
3eab989d7b
Function catand returned an incorrect value for its real part. To be possibly replaced with the intrinsic function atan.
101 lines
2.6 KiB
Fortran
101 lines
2.6 KiB
Fortran
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
|