drj – carlson-elliptic-integrals

public function drj(x, y, z, p, ier)

Compute an approximation for the incomplete or complete elliptic integral of the 3rd kind: $R_J(x,y,z,p) = \frac{3}{2} \int_{0}^{\infty} (t+x)^{-1/2} (t+y)^{-1/2} (t+z)^{-1/2} (t+p)^{-1} dt$ where $x\ge0$ , $y\ge0$ , and $z\ge0$ , and at most one of them $=0$ , and $p>0$ .

If $x=0$ or $y=0$ or $z=0$ , then 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.

DRJ Special Comments

$R_J(x,x+z,x+w,x+p) + R_J(y,y+z,y+w,y+p) + (a-b) R_J(a,b,b,a) + \frac{3}{\sqrt{a}} = R_J(0,z,w,p)$ where: $x>0, y>0, z>0, w>0, p>0$ and: $\begin{array}{rl} x y &= z w \\ a &= p^2 (x+y+z+w) \\ b &= p (p+x) (p+y) \\ b - a &= p (p-z) (p-w) \end{array}$ The sum of the third and fourth terms on the left side is $3 R_C(a,b)$ .

Special functions via DRJ and DRF

Legendre form of ELLIPTIC INTEGRAL of 3rd kind:

$\begin{array}{rcl} P(\phi,k,n) &=& \int_{0}^{\phi} (1+n \sin^2 \theta )^{-1} (1-k^2 \sin^2 \theta )^{-1/2} d \theta \\ &=& \sin \phi R_F(\cos^2 \phi, 1-k^2 sin^2 \phi,1) -\frac{n}{3} \sin^3 \phi R_J(\cos^2 \phi,1-k^2 \sin^2 \phi, 1,1+n \sin^2 \phi) \end{array}$

Bulirsch form of ELLIPTIC INTEGRAL of 3rd kind:

$\begin{array}{rcl} \mathrm{EL3}(x,k_c,p) &=& x R_F(1,1+k_c^2 x^2 ,1+x^2 ) + \frac{1}{3}(1-p) x^3 R_J(1,1+k_c^2 x^2 ,1+x^2 ,1+p x^2 ) \\ \mathrm{CEL}(k_c,p,a,b) &=& a R_F(0,k_c^2 ,1) + \frac{1}{3}(b-pa) R_J(0,k_c^2 ,1,p) \end{array}$

Heuman's LAMBDA function:

$\begin{array}{rcl} L(a,b,p) &=& ( \cos^2 a \sin b \cos b /(1-\cos^2 a \sin^2 b )^{1/2} ) (\sin p R_F(\cos^2 p ,1-\sin^2 a \sin^2 p ,1) \\ & ~ & + (\sin^2 a \sin^3 p /(3(1-\cos^2 a \sin^2 b ))) R_J(\cos^2 p ,1-\sin^2 a \sin^2 p ,1,1- \sin^2 a \sin^2 p /(1-\cos^2 a \sin^2 b ))) \\ \frac{\pi}{2} \Lambda_0(a,b) &=& L(a,b,\pi/2) \\ &=& \cos^2 a \sin b \cos b (1-\cos^2 a \sin^2 b )^{-1/2} R_F(0,\cos^2 a ,1) \\ &~& + (1/3) \sin^2 a \cos^2 a \sin b \cos b (1-\cos^2 a \sin^2 b )^{-3/2} R_J(0,\cos^2 a ,1,\cos^2 a \cos^2 b /(1-\cos^2 a \sin^2 b )) \end{array}$

Jacobi ZETA function:

$Z(b,k) = \frac{k^2}{3} \sin b \cos b (1-k^2 \sin^2 b)^{1/2} \frac{R_J(0,1-k^2 ,1,1-k^2 \sin^2 b)}{R_F (0,1-k^2 ,1)}$

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)
891009 Removed unreferenced statement labels. (WRB)
891009 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 Changed calls to XERMSG to standard form, and some editorial changes. (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
real(kind=wp),	intent(in)	::	y	nonnegative variable
real(kind=wp),	intent(in)	::	z	nonnegative variable
real(kind=wp),	intent(in)	::	p	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 = 1`: `min(x,y,z) < 0.0_wp` `IER = 2`: `min(x+y,x+z,y+z,p) < LOLIM` `IER = 3`: `max(x,y,z,p) > UPLIM`

Return Value real(kind=wp)

Calls

Help

Source Code

drj

Source Code

    real(wp) function drj(x,y,z,p,ier)

    implicit none

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

    character(len=16) xern3 , xern4 , xern5 , xern6 , xern7
    real(wp) :: alfa , beta , ea , eb , ec , e2 , e3, epslon
    real(wp) :: lamda , mu , pn , pndev
    real(wp) :: power4 , sigma , s1 , s2 , s3 , 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 of the series for DRJ
        !! 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.

    ! LOLIM and UPLIM determine the valid range of X, Y, Z, and P
    real(wp),parameter :: lolim  = (5.0_wp*d1mach(1))**(1.0_wp/3.0_wp)
        !! not less than the cube root of the value
        !! of LOLIM used in the routine for [[DRC]].
    real(wp),parameter :: uplim  = 0.30_wp*(d1mach(2)/5.0_wp)**(1.0_wp/3.0_wp)
        !! not greater than 0.3 times the cube root of
        !! the value of UPLIM used in the routine for [[DRC]].

    real(wp),parameter :: c1     = 3.0_wp/14.0_wp
    real(wp),parameter :: c2     = 1.0_wp/3.0_wp
    real(wp),parameter :: c3     = 3.0_wp/22.0_wp
    real(wp),parameter :: c4     = 3.0_wp/26.0_wp

    !initialize:
    drj = 0.0_wp

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

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

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

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

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

    ea  = xndev*(yndev+zndev) + yndev*zndev
    eb  = xndev*yndev*zndev
    ec  = pndev*pndev
    e2  = ea - 3.0_wp*ec
    e3  = eb + 2.0_wp*pndev*(ea-ec)
    s1  = 1.0_wp + e2*(-c1+0.750_wp*c3*e2-1.50_wp*c4*e3)
    s2  = eb*(0.50_wp*c2+pndev*(-c3-c3+pndev*c4))
    s3  = pndev*ea*(c2-pndev*c3) - c2*pndev*ec
    drj = 3.0_wp*sigma + power4*(s1+s2+s3)/(mu*sqrt(mu))

    end function drj

drj Function

Contents

Source Code

public function drj(x, y, z, p, ier)

DRJ Special Comments

Special functions via DRJ and DRF

Authors

References

History

Warning

Arguments

Return Value real(kind=wp)

Calls

Contents

Source Code

Source Code