drf Function

  1. carlson_elliptic_module.f90
  2. carlson_elliptic_module
  3. drf

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

Compute an approximation for the incomplete or complete elliptic integral of the 1st kind: Where , , , and at most one of them is .

If , , or , 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.

DRF Special Comments

Special functions via DRF

  • Legendre form of ELLIPTIC INTEGRAL of 1st kind:

  • Bulirsch form of ELLIPTIC INTEGRAL of 1st kind:

  • Lemniscate constant A:

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)

Arguments

TypeIntentOptionalAttributesName
real(kind=wp), intent(in) :: x

nonnegative variable
real(kind=wp), intent(in) :: y

nonnegative variable
real(kind=wp), intent(in) :: z

nonnegative 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,z) < 0
  • IER = 2:min(x+y,x+z,y+z) < LOLIM
  • IER = 3: max(x,y,z) > UPLIM

Return Value real(kind=wp)


Contents

Source Code

drf

Source Code

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

    implicit none

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

    character(len=16) :: xern3 , xern4 , xern5 , xern6
    real(wp) :: epslon, e2 , e3 , lamda
    real(wp) :: mu , s , xn , xndev
    real(wp) :: xnroot , yn , yndev , ynroot , zn , zndev , znroot

    real(wp),parameter :: errtol = (4.0_wp*d1mach(3))**(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
        !! `ERRTOL ** 6 / (4 * (1-ERRTOL)`.
        !!
        !! 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, 3.0e-19),
        !! (3.0e-3, 2.0e-16),
        !! (1.0e-2, 3.0e-13),
        !! (3.0e-2, 2.0e-10),
        !! (1.0e-1, 3.0e-7)
        !!
        !! 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/24.0_wp
    real(wp),parameter :: c2     = 3.0_wp/44.0_wp
    real(wp),parameter :: c3     = 1.0_wp/14.0_wp

    ! initialize:
    drf = 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)') 'drf: min(x,y,z)<0 where x = '// &
                xern3//' y = '//xern4//' and z = '//xern5
        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)') 'drf: max(x,y,z)>uplim where x = '// &
                xern3//' y = '//xern4//' z = '//xern5// &
                ' and uplim = '//xern6
        return
    endif

    if ( min(x+y,x+z,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)') 'drf: min(x+y,x+z,y+z)<lolim where x = '//xern3// &
                ' y = '//xern4//' z = '//xern5//' and lolim = '//xern6
        return
    endif

    ier = 0
    xn  = x
    yn  = y
    zn  = z

    do
        mu     = (xn+yn+zn)/3.0_wp
        xndev  = 2.0_wp - (mu+xn)/mu
        yndev  = 2.0_wp - (mu+yn)/mu
        zndev  = 2.0_wp - (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
        xn     = (xn+lamda)*0.250_wp
        yn     = (yn+lamda)*0.250_wp
        zn     = (zn+lamda)*0.250_wp
    end do
    e2  = xndev*yndev - zndev*zndev
    e3  = xndev*yndev*zndev
    s   = 1.0_wp + (c1*e2-0.10_wp-c2*e3)*e2 + c3*e3
    drf = s/sqrt(mu)

    end function drf