drd – carlson-elliptic-integrals

public function drd(x, y, z, ier)

Compute an approximation for the incomplete or complete elliptic integral of the 2nd kind: $R_D(x,y,z) = \frac{3}{2} \int_{0}^{\infty} (t+x)^{-1/2} (t+y)^{-1/2} (t+z)^{-3/2} dt$ Where $x\ge0$ , $y\ge0$ , $x+y>0$ , and $z>0$ .

If $x=0$ or $y=0$ , the integral is complete.

The duplication theorem is iterated until the variables are nearly equal, and the function is then expanded in Taylor series to fifth order.

DRD Special Comments

$\begin{array}{rl} R_D(x,y,z) + R_D(y,z,x) + R_D(z,x,y) = \frac{3}{\sqrt{x y z}} & x>0, y>0, z>0 \end{array}$

Special functions via DRD and DRF

Legendre form of ELLIPTIC INTEGRAL of 2nd kind:

$E(\phi,k) = \sin \phi R_F(\cos^2 \phi,1-k^2 \sin^2 \phi,1) -\frac{k^2}{3} \sin^3 \phi R_D(\cos^2 \phi,1-k^2 \sin^2 \phi,1)$ When $\phi = \pi /2$ the integral is complete: $\begin{array}{rcl} E(k) &=& R_F(0,1-k^2 ,1) - \frac{k^2}{3} R_D(0,1-k^2 ,1) \\ &=& \int_{0}^{\pi/2} \sqrt{1-k^2 \sin^2 \phi} d \phi \end{array}$
Bulirsch form of ELLIPTIC INTEGRAL of 2nd kind:

$\mathrm{EL2}(x,k_c,a,b) = ax R_F(1,1+k_c^2 x^2 ,1+x^2 ) + \frac{1}{3}(b-a) x^3 R_D(1,1+k_c^2 x^2 ,1+x^2 )$
Legendre form of alternative ELLIPTIC INTEGRAL of 2nd kind:

$\begin{array}{rcl} D(q,k) &=& \int_{0}^{q} \sin^2 p (1-k^2 \sin^2 p)^{-1/2} dp \\ &=& \frac{1}{3} (\sin^3 q) R_D(\cos^2 q,1-k^2 \sin^2 q,1) \end{array}$
Lemniscate constant B:

$\begin{array}{rcl} B &=& \int_{0}^{1} s^2 (1-s^4)^{-1/2} ds \\ &=& \frac{1}{3} R_D (0,2,1) \end{array}$
Heuman's LAMBDA function:

$\begin{array}{rcl} \frac{\pi}{2} \Lambda_0(a,b) &=& \sin b \left(R_F(0,\cos^2 a,1)-\frac{1}{3} \sin^2 a R_D(0,\cos^2 a,1) \right) R_F(\cos^2 b,1-\cos^2 a \sin^2 b,1) \\ & &-\frac{1}{3} \cos^2 a \sin^3 b R_F(0,\cos^2 a,1) R_D(\cos^2 b,1-\cos^2 a \sin^2 b,1) \end{array}$
Jacobi ZETA function:

$\begin{array}{rcl} Z(b,k) &=& \frac{k^2}{3} \sin b R_F(\cos^2 b,1-k^2 \sin^2 b,1) R_D(0,1-k^2 ,1)/R_F(0,1-k^2 ,1) \\ & & -\frac{k^2}{3} \sin^3 b R_D(\cos^2 b,1-k^2 \sin^2 b,1) \end{array}$

Authors

Carlson, B. C. Ames Laboratory-DOE, Iowa State University, Ames, IA 50011
Notis, E. M., Ames Laboratory-DOE, Iowa State University, Ames, IA 50011
Pexton, R. L., Lawrence Livermore National Laboratory, Livermore, CA 94550

References

B. C. Carlson and E. M. Notis, Algorithms for incomplete elliptic integrals, ACM Transactions on Mathematical Software 7, 3 (September 1981), pp. 398-403.
B. C. Carlson, Computing elliptic integrals by duplication, Numerische Mathematik 33, (1979), pp. 1-16.
B. C. Carlson, Elliptic integrals of the first kind, SIAM Journal of Mathematical Analysis 8, (1977), pp. 231-242.

History

790801 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)
900510 Modify calls to XERMSG to put in standard form. (RWC)
920501 Reformatted the REFERENCES section. (WRB)
Jan 2016, Refactored SLATEC routine into modern Fortran. (Jacob Williams)

Warning

Changes in the program may improve speed at the expense of robustness.

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	x	nonnegative variable ( $x+y>0$ )
real(kind=wp),	intent(in)	::	y	nonnegative variable ( $x+y>0$ )
real(kind=wp),	intent(in)	::	z	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`: `min(x,y) < 0` `IER = 2`: `min(x + y, z ) < LOLIM` `IER = 3`: `max(x,y,z) > UPLIM`

Return Value real(kind=wp)

Source Code

drd

Source Code

    real(wp) function drd(x,y,z,ier)

    implicit none

    real(wp),intent(in) :: x    !! nonnegative variable (\(x+y>0\))
    real(wp),intent(in) :: y    !! nonnegative variable (\(x+y>0\))
    real(wp),intent(in) :: z    !! 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`: `min(x,y) < 0`
                                !! * `IER = 2`: `min(x + y, z ) < LOLIM`
                                !! * `IER = 3`: `max(x,y,z) > UPLIM`

    character(len=16) :: xern3 , xern4 , xern5 , xern6
    real(wp) :: epslon, ea , eb , ec , ed , ef , lamda
    real(wp) :: mu , power4 , sigma , s1 , s2 , xn , xndev
    real(wp) :: xnroot , yn , yndev , ynroot , zn , zndev , znroot

    real(wp),parameter :: errtol = (d1mach(3)/3.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
        !! `3 * ERRTOL ** 6 / (1-ERRTOL) ** 3/2`.
        !!
        !! The accuracy of the computed approximation to the integral
        !! 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 truncation
        !! error there will be round-off error, but in practice the
        !! total error from both sources is usually less than the
        !! amount given in the table.
        !!
        !! Sample choices:
        !! (ERRTOL, Relative truncation error less than):
        !! (1.0e-3, 4.0e-18),
        !! (3.0e-3, 3.0e-15),
        !! (1.0e-2, 4.0e-12),
        !! (3.0e-2, 3.0e-9),
        !! (1.0e-1, 4.0e-6)
        !!
        !! 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  = 2.0_wp/(d1mach(2))**(2.0_wp/3.0_wp) !! Lower limit of valid arguments
    real(wp),parameter :: tuplim = (0.10_wp*errtol)**(1.0_wp/3.0_wp)/&
                                    d1mach(1)**(1.0_wp/3.0_wp)
    real(wp),parameter :: uplim  = tuplim**2  !! Upper limit of valid arguments
    real(wp),parameter :: c1     = 3.0_wp/14.0_wp
    real(wp),parameter :: c2     = 1.0_wp/6.0_wp
    real(wp),parameter :: c3     = 9.0_wp/22.0_wp
    real(wp),parameter :: c4     = 3.0_wp/26.0_wp

    ! initialize:
    drd = 0.0_wp

    ! check for errors:
    if ( min(x,y)<0.0_wp ) then
        ier = 1
        write (xern3,'(1pe15.6)') x
        write (xern4,'(1pe15.6)') y
        write(error_unit,'(a)') 'drd: min(x,y)<0 where x = '//xern3// &
                                ' and y = '//xern4
        return
    endif

    if ( max(x,y,z)>uplim ) then
        ier = 3
        write (xern3,'(1pe15.6)') x
        write (xern4,'(1pe15.6)') y
        write (xern5,'(1pe15.6)') z
        write (xern6,'(1pe15.6)') uplim
        write(error_unit,'(a)') 'drd: max(x,y,z)>uplim where x = '// &
                                xern3//' y = '//xern4//' z = '//xern5// &
                                ' and uplim = '//xern6
        return
    endif

    if ( min(x+y,z)<lolim ) then
        ier = 2
        write (xern3,'(1pe15.6)') x
        write (xern4,'(1pe15.6)') y
        write (xern5,'(1pe15.6)') z
        write (xern6,'(1pe15.6)') lolim
        write(error_unit,'(a)') 'drd: min(x+y,z)<lolim where x = '// &
                                xern3//' y = '//xern4//' z = '//xern5// &
                                ' and lolim = '//xern6
        return
    endif

    ier    = 0
    xn     = x
    yn     = y
    zn     = z
    sigma  = 0.0_wp
    power4 = 1.0_wp

    do
        mu     = (xn+yn+3.0_wp*zn)*0.20_wp
        xndev  = (mu-xn)/mu
        yndev  = (mu-yn)/mu
        zndev  = (mu-zn)/mu
        epslon = max(abs(xndev),abs(yndev),abs(zndev))
        if ( epslon<errtol ) exit
        xnroot = sqrt(xn)
        ynroot = sqrt(yn)
        znroot = sqrt(zn)
        lamda  = xnroot*(ynroot+znroot) + ynroot*znroot
        sigma  = sigma + power4/(znroot*(zn+lamda))
        power4 = power4*0.250_wp
        xn     = (xn+lamda)*0.250_wp
        yn     = (yn+lamda)*0.250_wp
        zn     = (zn+lamda)*0.250_wp
    end do

    ea  = xndev*yndev
    eb  = zndev*zndev
    ec  = ea - eb
    ed  = ea - 6.0_wp*eb
    ef  = ed + ec + ec
    s1  = ed*(-c1+0.250_wp*c3*ed-1.50_wp*c4*zndev*ef)
    s2  = zndev*(c2*ef+zndev*(-c3*ec+zndev*c4*ea))
    drd = 3.0_wp*sigma + power4*(1.0_wp+s1+s2)/(mu*sqrt(mu))

    end function drd

drd Function

Contents

Source Code

public function drd(x, y, z, ier)

DRD Special Comments

Special functions via DRD and DRF

Authors

References

History

Warning

Arguments

Return Value real(kind=wp)

Contents

Source Code

Source Code