function ellipsoid_area ( a, b, c )

!*****************************************************************************80
!
!! ellipsoid_area() returns the surface area of an ellipsoid.
!
!  Discussion:
!
!    The ellipsoid may be represented by the equation
!
!      (x/a)^2 + (y/b)^2 + (z/c)^2 = 1
!
!    with a => b => c
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    06 January 2023
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    John D Cook,
!    Ellipsoid surface area,
!    6 July 2014,
!    https:!www.johndcook.com/blog/2014/07/06/ellipsoid-surface-area/
!
!  Input:
!
!    real ( kind = rk ) a, b, c: the semi-axes of the ellipsoid, with a => b => c
!
!  Output:
!
!    real ( kind = rk ) ellipsoid_area: the surface area of the ellipsoid.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ) ellipsoid_area
  real ( kind = rk ) elliptic_inc_em
  real ( kind = rk ) elliptic_inc_fm
  real ( kind = rk ) m
  real ( kind = rk ) phi
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) t
  real ( kind = rk ) temp
  real ( kind = rk ) temp2

  a = abs ( a )
  b = abs ( b )
  c = abs ( c )
!
!  Force a => b => c;
!
  if ( a < b ) then
    t = a
    a = b
    b = t
  end if

  if ( a < c ) then
    t = a
    a = c
    c = t
  end if

  if ( b < c ) then
    t = b
    b = c
    c = t
  end if

  phi = acos ( c / a )

  if ( a * a - c * c == 0.0D+00 ) then
    m = 1.0D+00
  else
    m = ( a*a * ( b*b - c*c ) ) / ( b*b * ( a*a - c*c ) )
  end if

  temp = elliptic_inc_em ( phi, m ) * ( sin ( phi ) ) ** 2 &
       + elliptic_inc_fm ( phi, m ) * ( cos ( phi ) ) ** 2

  if ( sin ( phi ) == 0.0D+00 ) then
    temp2 = 1.0D+00
  else
    temp2 = temp / sin ( phi )
  end if

  ellipsoid_area = 2.0D+00 * r8_pi * ( c*c + a * b * temp2 )
  
  return
end
function ellipsoid_volume ( a, b, c )

!*****************************************************************************80
!
!! ellipsoid_volume() returns the volume of an ellipsoid.
!
!  Discussion:
!
!    The ellipsoid may be represented by the equation
!
!      (x/a)^2 + (y/b)^2 + (z/c)^2 = 1
!
!    with a => b => c
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    06 January 2023
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk ) a, b, c: the semi-axes of the ellipsoid, with a => b => c
!
!  Output:
!
!    real ( kind = rk ) ellipsoid_volume: the volume of the ellipsoid.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ) ellipsoid_volume
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00

  a = abs ( a )
  b = abs ( b )
  c = abs ( c )

  ellipsoid_volume = ( 4.0D+00 / 3.0D+00 ) * r8_pi * ( a * b * c )
  
  return
end
function elliptic_inc_em ( phi, m )

!*****************************************************************************80
!
!! elliptic_inc_em() evaluates the incomplete elliptic integral E(PHI,M).
!
!  Discussion:
!
!    The value is computed using Carlson elliptic integrals:
!
!      E(phi,m) = 
!                sin ( phi )   RF ( cos^2 ( phi ), 1-m sin^2 ( phi ), 1 ) 
!        - 1/3 m sin^3 ( phi ) RD ( cos^2 ( phi ), 1-m sin^2 ( phi ), 1 ).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 May 2018
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk ) PHI, K, the argument.
!    0 <= PHI <= PI/2.
!    0 <= M * sin^2(PHI) <= 1.
!
!  Output:
!
!    real ( kind = rk ) ELLIPTIC_INC_EM, the function value.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) cp
  real ( kind = rk ) elliptic_inc_em
  real ( kind = rk ) errtol
  integer ierr
  real ( kind = rk ) m
  real ( kind = rk ) phi
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) rd
  real ( kind = rk ) rf
  real ( kind = rk ) sp
  real ( kind = rk ) value
  real ( kind = rk ) value1
  real ( kind = rk ) value2
  real ( kind = rk ) x
  real ( kind = rk ) y
  real ( kind = rk ) z

  cp = cos ( phi )
  sp = sin ( phi )
  x = cp * cp
  y = 1.0D+00 - m * sp * sp
  z = 1.0D+00
  errtol = 1.0D-03

  value1 = rf ( x, y, z, errtol, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'elliptic_inc_em(): Fatal error!'
    write ( *, '(a,i2)' ) '  RF returned IERR = ', ierr
    stop 1
  end if

  value2 = rd ( x, y, z, errtol, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'elliptic_inc_em(): Fatal error!'
    write ( *, '(a,i2)' ) '  RD returned IERR = ', ierr
    stop 1
  end if

  value = sp * value1 - m * sp ** 3 * value2 / 3.0D+00

  elliptic_inc_em = value

  return
end
function elliptic_inc_fm ( phi, m )

!*****************************************************************************80
!
!! elliptic_inc_fm() evaluates the incomplete elliptic integral F(PHI,M).
!
!  Discussion:
!
!    The value is computed using Carlson elliptic integrals:
!
!      F(phi,m) = sin(phi) * RF ( cos^2 ( phi ), 1-m sin^2 ( phi ), 1 )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 June 2018
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk ) PHI, M, the argument.
!    0 <= PHI <= PI/2.
!    0 <= M * sin^2(PHI) <= 1.
!
!  Output:
!
!    real ( kind = rk ) ELLIPTIC_INC_FM, the function value.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) cp
  real ( kind = rk ) elliptic_inc_fm
  real ( kind = rk ) errtol
  integer ierr
  real ( kind = rk ) m
  real ( kind = rk ) phi
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) rf
  real ( kind = rk ) sp
  real ( kind = rk ) value
  real ( kind = rk ) x
  real ( kind = rk ) y
  real ( kind = rk ) z

  cp = cos ( phi )
  sp = sin ( phi )
  x = cp * cp
  y = 1.0D+00 - m * sp ** 2
  z = 1.0D+00
  errtol = 1.0D-03

  value = rf ( x, y, z, errtol, ierr )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'elliptic_inc_fm(): Fatal error!'
    write ( *, '(a,i2)' ) '  RF returned IERR = ', ierr
    stop 1
  end if

  elliptic_inc_fm = sp * value

  return
end
function rd ( x, y, z, errtol, ierr )

!*****************************************************************************80
!
!! rd() computes an incomplete elliptic integral of the second kind, RD(X,Y,Z).
!
!  Discussion:
!
!    This function computes an incomplete elliptic integral of the second kind.
!
!    RD(X,Y,Z) = Integral ( 0 <= T < oo )
!
!                    -1/2     -1/2     -3/2
!          (3/2)(T+X)    (T+Y)    (T+Z)    DT,
!
!    where X and Y are nonnegative, X + Y is positive, and Z is positive.
!
!    If X or Y is zero, 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.
!
!    Check:
!
!      RD(X,Y,Z) + RD(Y,Z,X) + RD(Z,X,Y) = 3 / sqrt ( X * Y * Z ),
!      where X, Y, and Z are positive.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 May 2018
!
!  Author:
!
!    Original FORTRAN77 version by Bille Carlson, Elaine Notis.
!    This version by John Burkardt.
!
!  Reference:
!
!    Bille Carlson,
!    Computing Elliptic Integrals by Duplication,
!    Numerische Mathematik,
!    Volume 33, 1979, pages 1-16.
!
!    Bille Carlson, Elaine Notis,
!    Algorithm 577, Algorithms for Incomplete Elliptic Integrals,
!    ACM Transactions on Mathematical Software,
!    Volume 7, Number 3, pages 398-403, September 1981.
!
!  Input:
!
!    real ( kind = rk ) X, Y, Z, the arguments in the integral.
!
!    real ( kind = rk ) ERRTOL, the error tolerance.
!    The relative error due to truncation is less than
!      3 * ERRTOL ^ 6 / (1-ERRTOL) ^ 3/2.
!    Sample choices:
!      ERRTOL   Relative truncation error less than
!      1.D-3    4.D-18
!      3.D-3    3.D-15
!      1.D-2    4.D-12
!      3.D-2    3.D-9
!      1.D-1    4.D-6
!
!  Output:
!
!    integer IERR, the error flag.
!    0, no error occurred.
!    1, abnormal termination.
!
!    real ( kind = rk ) rd: the value of the elliptic integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) c1
  real ( kind = rk ) c2
  real ( kind = rk ) c3
  real ( kind = rk ) c4
  real ( kind = rk ) ea
  real ( kind = rk ) eb
  real ( kind = rk ) ec
  real ( kind = rk ) ed
  real ( kind = rk ) ef
  real ( kind = rk ) epslon
  real ( kind = rk ) errtol
  integer ierr
  real ( kind = rk ) lamda
  real ( kind = rk ) lolim
  real ( kind = rk ) mu
  real ( kind = rk ) power4
  real ( kind = rk ) rd
  real ( kind = rk ) sigma
  real ( kind = rk ) s1
  real ( kind = rk ) s2
  real ( kind = rk ) uplim
  real ( kind = rk ) x
  real ( kind = rk ) xn
  real ( kind = rk ) xndev
  real ( kind = rk ) xnroot
  real ( kind = rk ) y
  real ( kind = rk ) yn
  real ( kind = rk ) yndev
  real ( kind = rk ) ynroot
  real ( kind = rk ) z
  real ( kind = rk ) zn
  real ( kind = rk ) zndev
  real ( kind = rk ) znroot
!
!  lolim and uplim determine the range of valid arguments.
!  lolim is not less than 2 / (machine maximum) ^ (2/3).
!  uplim is not greater than (0.1 * errtol / machine
!  minimum) ^ (2/3), where errtol is described below.
!  in the following table it is assumed that errtol will
!  never be chosen smaller than 1.d-5.
!
  save lolim
  save uplim

  data lolim /6.D-51/
  data uplim /1.D+48/

  if ( &
    x < 0.0D+00 .or. &
    y < 0.0D+00 .or. &
    x + y < lolim .or. &
    z < lolim .or. &
    uplim < x .or. &
    uplim < y .or. &
    uplim < z ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'rd(): Fatal error!'
    write ( *, '(a)' ) '  Invalid input arguments.'
    write ( *, '(a,d23.16)' ) '  X = ', x
    write ( *, '(a,d23.16)' ) '  Y = ', y
    write ( *, '(a,d23.16)' ) '  Z = ', z
    write ( *, '(a)' ) ''
    ierr = 1
    rd = 0.0D+00
    return
  end if

  ierr = 0
  xn = x
  yn = y
  zn = z
  sigma = 0.0d0
  power4 = 1.0d0

  do

    mu = ( xn + yn + 3.0d0 * zn ) * 0.2d0
    xndev = ( mu - xn ) / mu
    yndev = ( mu - yn ) / mu
    zndev = ( mu - zn ) / mu
    epslon = max ( abs ( xndev ), abs ( yndev ), abs ( zndev ) )

    if ( epslon < errtol ) then
      c1 = 3.0d0 / 14.0d0
      c2 = 1.0d0 / 6.0d0
      c3 = 9.0d0 / 22.0d0
      c4 = 3.0d0 / 26.0d0
      ea = xndev * yndev
      eb = zndev * zndev
      ec = ea - eb
      ed = ea - 6.0d0 * eb
      ef = ed + ec + ec
      s1 = ed * ( - c1 + 0.25d0 * c3 * ed - 1.5d0 * c4 * zndev * ef )
      s2 = zndev * ( c2 * ef + zndev * ( - c3 * ec + zndev * c4 * ea ) )
      rd = 3.0d0 * sigma + power4 * ( 1.0d0 + s1 + s2 ) / ( mu * sqrt ( mu ) )

      return
    end if

    xnroot = sqrt ( xn )
    ynroot = sqrt ( yn )
    znroot = sqrt ( zn )
    lamda = xnroot * ( ynroot + znroot ) + ynroot * znroot
    sigma = sigma + power4 / ( znroot * ( zn + lamda ) )
    power4 = power4 * 0.25d0
    xn = ( xn + lamda ) * 0.25d0
    yn = ( yn + lamda ) * 0.25d0
    zn = ( zn + lamda ) * 0.25d0

  end do

end
function rf ( x, y, z, errtol, ierr )

!*****************************************************************************80
!
!! rf() computes an incomplete elliptic integral of the first kind, RF(X,Y,Z).
!
!  Discussion:
!
!    This function computes the incomplete elliptic integral of the first kind.
!
!    RF(X,Y,Z) = Integral ( 0 <= T < oo )
!
!                    -1/2     -1/2     -1/2
!          (1/2)(T+X)    (T+Y)    (T+Z)    DT,
!
!    where X, Y, and Z are nonnegative and at most one of them is zero.
!
!    If X or Y or Z is zero, 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.
!
!    Check by addition theorem:
!
!      RF(X,X+Z,X+W) + RF(Y,Y+Z,Y+W) = RF(0,Z,W),
!      where X, Y, Z, W are positive and X * Y = Z * W.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 May 2018
!
!  Author:
!
!    Original FORTRAN77 version by Bille Carlson, Elaine Notis.
!    This version by John Burkardt.
!
!  Reference:
!
!    Bille Carlson,
!    Computing Elliptic Integrals by Duplication,
!    Numerische Mathematik,
!    Volume 33, 1979, pages 1-16.
!
!    Bille Carlson, Elaine Notis,
!    Algorithm 577, Algorithms for Incomplete Elliptic Integrals,
!    ACM Transactions on Mathematical Software,
!    Volume 7, Number 3, pages 398-403, September 1981.
!
!  Input:
!
!    real ( kind = rk ) X, Y, Z, the arguments in the integral.
!
!    real ( kind = rk ) ERRTOL, the error tolerance.
!    Relative error due to truncation is less than
!      ERRTOL ^ 6 / (4 * (1 - ERRTOL)).
!    Sample choices:
!      ERRTOL   Relative truncation error less than
!      1.D-3    3.D-19
!      3.D-3    2.D-16
!      1.D-2    3.D-13
!      3.D-2    2.D-10
!      1.D-1    3.D-7
!
!  Output:
!
!    integer IERR, the error flag.
!    0, no error occurred.
!    1, abnormal termination.
!
!    real ( kind = rk ) rf: the value of the elliptic integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) c1
  real ( kind = rk ) c2
  real ( kind = rk ) c3
  real ( kind = rk ) e2
  real ( kind = rk ) e3
  real ( kind = rk ) epslon
  real ( kind = rk ) errtol
  integer ierr
  real ( kind = rk ) lamda
  real ( kind = rk ) lolim
  real ( kind = rk ) mu
  real ( kind = rk ) rf
  real ( kind = rk ) s
  real ( kind = rk ) uplim
  real ( kind = rk ) x
  real ( kind = rk ) xn
  real ( kind = rk ) xndev
  real ( kind = rk ) xnroot
  real ( kind = rk ) y
  real ( kind = rk ) yn
  real ( kind = rk ) yndev
  real ( kind = rk ) ynroot
  real ( kind = rk ) z
  real ( kind = rk ) zn
  real ( kind = rk ) zndev
  real ( kind = rk ) znroot
!
!  lolim and uplim determine the range of valid arguments.
!  lolim is not less than the machine minimum multiplied by 5.
!  uplim is not greater than the machine maximum divided by 5.
!
  save lolim
  save uplim

  data lolim /3.D-78/
  data uplim /1.D+75/

  if ( &
    x < 0.0D+00 .or. &
    y < 0.0D+00 .or. &
    z < 0.0D+00 .or. &
    x + y < lolim .or. &
    x + z < lolim .or. &
    y + z < lolim .or. &
    uplim <= x .or. &
    uplim <= y .or. &
    uplim <= z ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'rf(): Fatal error!'
    write ( *, '(a)' ) '  Invalid input arguments.'
    write ( *, '(a,d23.16)' ) '  X = ', x
    write ( *, '(a,d23.16)' ) '  Y = ', y
    write ( *, '(a,d23.16)' ) '  Z = ', z
    write ( *, '(a)' ) ''
    ierr = 1
    rf = 0.0D+00
    return
  end if

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

  do

    mu = ( xn + yn + zn ) / 3.0d0
    xndev = 2.0d0 - ( mu + xn ) / mu
    yndev = 2.0d0 - ( mu + yn ) / mu
    zndev = 2.0d0 - ( mu + zn ) / mu
    epslon = max ( abs ( xndev ), abs ( yndev ), abs ( zndev ) )

    if ( epslon < errtol ) then
      c1 = 1.0d0 / 24.0d0
      c2 = 3.0d0 / 44.0d0
      c3 = 1.0d0 / 14.0d0
      e2 = xndev * yndev - zndev * zndev
      e3 = xndev * yndev * zndev
      s = 1.0d0 + ( c1 * e2 - 0.1d0 - c2 * e3 ) * e2 + c3 * e3
      rf = s / sqrt ( mu )
      return
    end if

    xnroot = sqrt ( xn )
    ynroot = sqrt ( yn )
    znroot = sqrt ( zn )
    lamda = xnroot * ( ynroot + znroot ) + ynroot * znroot
    xn = ( xn + lamda ) * 0.25d0
    yn = ( yn + lamda ) * 0.25d0
    zn = ( zn + lamda ) * 0.25d0

  end do

end