Compute an approximation of the Carlson elliptic integral: where and .
The duplication theorem is iterated until the variables are nearly equal, and the function is then expanded in Taylor series to fifth order. Logarithmic, inverse circular, and inverse hyperbolic functions can be expressed in terms of DRC.
Changes in the program may improve speed at the expense of robustness.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | x | nonnegative variable |
||
| real(kind=wp), | intent(in) | :: | y | positive variable |
||
| integer, | intent(out) | :: | ier | indicates normal or abnormal termination:
|
real(wp) function drc(x,y,ier)
implicit none
real(wp),intent(in) :: x !! nonnegative variable
real(wp),intent(in) :: y !! positive variable
integer,intent(out) :: ier !! indicates normal or abnormal termination:
!!
!! * `IER = 0`: Normal and reliable termination of the
!! routine. It is assumed that the requested
!! accuracy has been achieved.
!! * `IER > 0`: Abnormal termination of the routine:
!! * `IER = 1`: `x<0 or y<=0`
!! * `IER = 2`: `x+y<LOLIM`
!! * `IER = 3`: `max(x,y) > UPLIM`
character(len=16) :: xern3 , xern4 , xern5
real(wp) :: lamda, mu , s , sn , xn , yn
real(wp),parameter :: errtol = (d1mach(3)/16.0_wp)**(1.0_wp/6.0_wp)
!! Determines the accuracy of the answer.
!!
!! The value assigned by the routine will result
!! in solution precision within 1-2 decimals of
!! machine precision.
!!
!! Relative error due to truncation is less than
!! `16 * ERRTOL ** 6 / (1 - 2 * ERRTOL)`.
!!
!! Sample choices:
!! (ERRTOL, Relative truncation error less than):
!! (1.0e-3, 2.0e-17),
!! (3.0e-3, 2.0e-14),
!! (1.0e-2, 2.0e-11),
!! (3.0e-2, 2.0e-8),
!! (1.0e-1, 2.0e-5)
!!
!! The accuracy of the computed approximation to the inte-
!! gral can be controlled by choosing the value of ERRTOL.
!! Truncation of a Taylor series after terms of fifth order
!! introduces an error less than the amount shown in the
!! second column of the following table for each value of
!! ERRTOL in the first column. In addition to the trunca-
!! tion error there will be round-off error, but in prac-
!! tice the total error from both sources is usually less
!! than the amount given in the table.
!!
!! Decreasing ERRTOL by a factor of 10 yields six more
!! decimal digits of accuracy at the expense of one or
!! two more iterations of the duplication theorem.
real(wp),parameter :: lolim = 5.0_wp*d1mach(1) !! Lower limit of valid arguments
real(wp),parameter :: uplim = d1mach(2)/5.0_wp !! Upper limit of valid arguments
real(wp),parameter :: c1 = 1.0_wp/7.0_wp
real(wp),parameter :: c2 = 9.0_wp/22.0_wp
!initialize:
drc = 0.0_wp
! check for errors:
if ( x<0.0_wp .or. y<=0.0_wp ) then
ier = 1
write (xern3,'(1pe15.6)') x
write (xern4,'(1pe15.6)') y
write(error_unit,'(a)') &
'drc: x<0 .or. y<=0 where x = '//xern3//' and y = '//xern4
return
endif
if ( max(x,y)>uplim ) then
ier = 3
write (xern3,'(1pe15.6)') x
write (xern4,'(1pe15.6)') y
write (xern5,'(1pe15.6)') uplim
write(error_unit,'(a)') &
'drc: max(x,y)>uplim where x = '//&
xern3//' y = '//xern4//' and uplim = '//xern5
return
endif
if ( x+y<lolim ) then
ier = 2
write (xern3,'(1pe15.6)') x
write (xern4,'(1pe15.6)') y
write (xern5,'(1pe15.6)') lolim
write(error_unit,'(a)') &
'drc: x+y<lolim where x = '//xern3//&
' y = '//xern4//' and lolim = '//xern5
return
endif
ier = 0
xn = x
yn = y
do
mu = (xn+yn+yn)/3.0_wp
sn = (yn+mu)/mu - 2.0_wp
if ( abs(sn)<errtol ) exit
lamda = 2.0_wp*sqrt(xn)*sqrt(yn) + yn
xn = (xn+lamda)*0.250_wp
yn = (yn+lamda)*0.250_wp
end do
s = sn*sn*(0.30_wp+sn*(c1+sn*(0.3750_wp+sn*c2)))
drc = (1.0_wp+s)/sqrt(mu)
end function drc