125 lines
3.3 KiB
Fortran
125 lines
3.3 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)
|
|
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)
|
|
implicit none
|
|
integer, intent(in) :: k
|
|
real(wp_) :: fact
|
|
integer :: i
|
|
! Factorial function
|
|
fact=zero
|
|
if(k<0) return
|
|
fact=one
|
|
if(k==0) return
|
|
do i=1,k
|
|
fact=fact*i
|
|
end do
|
|
end function fact
|
|
|
|
function gamm(xx)
|
|
implicit none
|
|
real(wp_) :: gamm
|
|
real(wp_), intent(in) :: xx
|
|
! Returns the value Gamma(xx) for xx > 0.
|
|
INTEGER :: j
|
|
real(wp_) :: ser,tmp,x,y
|
|
real(wp_), parameter :: stp=2.5066282746310005_wp_
|
|
real(wp_), dimension(6), parameter :: cof=(/76.18009172947146_wp_, &
|
|
-86.50532032941677_wp_,24.01409824083091_wp_,-1.231739572450155_wp_, &
|
|
.1208650973866179e-2_wp_,-.5395239384953e-5_wp_/)
|
|
x=xx
|
|
y=x
|
|
tmp=x+5.5_wp_
|
|
tmp=(x+0.5_wp_)*log(tmp)-tmp
|
|
ser=1.000000000190015_wp_
|
|
do j=1,6
|
|
y=y+1._wp_
|
|
ser=ser+cof(j)/y
|
|
end do
|
|
gamm=exp(tmp)*(stp*ser/x)
|
|
end function gamm
|
|
|
|
end module math |