# include <cfloat>
# include <cmath>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "ellipsoid.hpp"

//****************************************************************************80

double ellipsoid_area ( double a, double b, double c )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    double a, b, c: the semi-axes of the ellipsoid, with a => b => c
//
//  Output:
//
//    double ellipsoid_area: the surface area of the ellipsoid.
//
{
  double m;
  double phi;
  const double r8_pi = 3.141592653589793;
  double t;
  double temp;
  double temp2;
  double value;

  a = fabs ( a );
  b = fabs ( b );
  c = fabs ( c );
/*
  Force a => b => c;
*/
  if ( a < b )
  {
    t = a;
    a = b;
    b = t;
  }

  if ( a < c )
  {
    t = a;
    a = c;
    c = t;
  }

  if ( b < c )
  {
    t = b;
    b = c;
    c = t;
  }

  phi = acos ( c / a );

  if ( a * a - c * c == 0 )
  {
    m = 1.0;
  }
  else
  {
    m = ( a*a * ( b*b - c*c ) ) / ( b*b * ( a*a - c*c ) );
  }

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

  if ( sin ( phi ) == 0.0 )
  {
    temp2 = 1.0;
  }
  else
  {
    temp2 = temp / sin ( phi );
  }

  value = 2.0 * r8_pi * ( c*c + a * b * temp2 );
  
  return value;
}
//****************************************************************************80

double ellipsoid_volume ( double a, double b, double c )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    double a, b, c: the semi-axes of the ellipsoid, with a => b => c
//
//  Output:
//
//    double ellipsoid_volume: the volume of the ellipsoid.
//
{
  const double r8_pi = 3.141592653589793;
  double value;

  a = fabs ( a );
  b = fabs ( b );
  c = fabs ( c );

  value = ( 4.0 / 3.0 ) * r8_pi * ( a * b * c );
  
  return value;
}
//****************************************************************************80

double elliptic_inc_em ( double phi, double m )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    25 June 2018
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double PHI, M, the arguments.
//    0 <= PHI <= PI/2.
//    0 <= M * sin^2(PHI) <= 1.
//
//  Output:
//
//    double ELLIPTIC_INC_EM, the function value.
//
{
  double cp;
  double errtol;
  int ierr;
  double sp;
  double value;
  double value1;
  double value2;
  double x;
  double y;
  double z;

  cp = cos ( phi );
  sp = sin ( phi );
  x = cp * cp;
  y = 1.0 - m * sp * sp;
  z = 1.0;
  errtol = 1.0E-03;

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

  if ( ierr != 0 )
  {
    cout << "\n";
    cout << "elliptic_inc_em(): Fatal error!\n"; 
    cout << "  rf() returned IERR = " << ierr << "\n";
    exit ( 1 );
  }

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

  if ( ierr != 0 )
  {
    cout << "\n";
    cout << "elliptic_inc_em (): Fatal error!\n";
    cout << "  rd() returned IERR = " << ierr << "\n";
    exit ( 1 );
  }

  value = sp * value1 - m * sp * sp * sp * value2 / 3.0;

  return value;
}
//****************************************************************************80

double elliptic_inc_fm ( double phi, double m )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    25 June 2018
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double PHI, M, the arguments.
//    0 <= PHI <= PI/2.
//    0 <= M * sin^2(PHI) <= 1.
//
//  Output:
//
//    double ELLIPTIC_INC_FM, the function value.
//
{
  double cp;
  double errtol;
  int ierr;
  double sp;
  double value;
  double x;
  double y;
  double z;

  cp = cos ( phi );
  sp = sin ( phi );
  x = cp * cp;
  y = 1.0 - m * sp * sp;
  z = 1.0;
  errtol = 1.0E-03;

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

  if ( ierr != 0 )
  {
    cout << "\n";
    cout << "elliptic_inc_fm(): Fatal error!\n"; 
    cout << "  rf() returned IERR = " << ierr << "\n";
    exit ( 1 );
  }

  value = sp * value;

  return value;
}
//****************************************************************************80

double rd ( double x, double y, double z, double errtol, int &ierr )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    02 June 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:
//
//    double X, Y, Z, the arguments in the integral.
//
//    double 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:
//
//    int &IERR, the error flag.
//    0, no error occurred.
//    1, abnormal termination.
//
//    double rd: the value of the elliptic integral.
//
{
  double c1;
  double c2;
  double c3;
  double c4;
  double ea;
  double eb;
  double ec;
  double ed;
  double ef;
  double epslon;
  double lamda;
  const double lolim = 6.0E-51;
  double mu;
  double power4;
  double sigma;
  double s1;
  double s2;
  const double uplim = 1.0E+48;
  double value;
  double xn;
  double xndev;
  double xnroot;
  double yn;
  double yndev;
  double ynroot;
  double zn;
  double zndev;
  double 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.
//
  if ( 
    x < 0.0 || 
    y < 0.0 || 
    x + y < lolim || 
    z < lolim || 
    uplim < x || 
    uplim < y || 
    uplim < z )
  {
    cout << "\n";
    cout << "rd(): Error!\n";
    cout << "  Invalid input arguments.\n";
    cout << "  X = " << x << "\n";
    cout << "  Y = " << y << "\n";
    cout << "  Z = " << z << "\n";
    cout << "\n";
    ierr = 1;
    value = 0.0;
    return value;
  }

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

  while ( true )
  {
    mu = ( xn + yn + 3.0 * zn ) * 0.2;
    xndev = ( mu - xn ) / mu;
    yndev = ( mu - yn ) / mu;
    zndev = ( mu - zn ) / mu;
    epslon = fmax ( fabs ( xndev ), 
      fmax ( fabs ( yndev ), fabs ( zndev ) ) );

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

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

}
//****************************************************************************80

double rf ( double x, double y, double z, double errtol, int &ierr )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    02 June 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:
//
//    double X, Y, Z, the arguments in the integral.
//
//    double 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:
//
//    int &IERR, the error flag.
//    0, no error occurred.
//    1, abnormal termination.
//
//    double rf: the value of the elliptic integral.
//
{
  double c1;
  double c2;
  double c3;
  double e2;
  double e3;
  double epslon;
  double lamda;
  const double lolim = 3.0E-78;
  double mu;
  double s;
  const double uplim = 1.0E+75;
  double value;
  double xn;
  double xndev;
  double xnroot;
  double yn;
  double yndev;
  double ynroot;
  double zn;
  double zndev;
  double 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.
//
  if ( 
    x < 0.0 || 
    y < 0.0 || 
    z < 0.0 || 
    x + y < lolim || 
    x + z < lolim || 
    y + z < lolim || 
    uplim <= x || 
    uplim <= y || 
    uplim <= z )
  {
    cout << "\n";
    cout << "rf(): Error!\n";
    cout << "  Invalid input arguments.\n";
    cout << "  X = " << x << "\n";
    cout << "  Y = " << y << "\n";
    cout << "  Z = " << z << "\n";
    cout << "\n";
    ierr = 1;
    value = 0.0;
    return value;
  }

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

  while ( true )
  {
    mu = ( xn + yn + zn ) / 3.0;
    xndev = 2.0 - ( mu + xn ) / mu;
    yndev = 2.0 - ( mu + yn ) / mu;
    zndev = 2.0 - ( mu + zn ) / mu;
    epslon = fmax ( fabs ( xndev ), 
      fmax ( fabs ( yndev ), fabs ( zndev ) ) );

    if ( epslon < errtol )
    {
      c1 = 1.0 / 24.0;
      c2 = 3.0 / 44.0;
      c3 = 1.0 / 14.0;
      e2 = xndev * yndev - zndev * zndev;
      e3 = xndev * yndev * zndev;
      s = 1.0 + ( c1 * e2 - 0.1 - c2 * e3 ) * e2 + c3 * e3;
      value = s / sqrt ( mu );
      return value;
    }

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

}