# include <float.h>
# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>

# include "ellipsoid.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    05 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;
}
/******************************************************************************/

double elliptic_inc_em ( double phi, double m )

/******************************************************************************/
/*
  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:

    05 January 2023

  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 )
  {
    printf ( "\n" );
    printf ( "elliptic_inc_em(): Fatal error!\n" ); 
    printf ( "  RF returned IERR = %d\n", ierr );
    exit ( 1 );
  }

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

  if ( ierr != 0 )
  {
    printf ( "\n" );
    printf ( "elliptic_inc_em(): Fatal error!\n" );
    printf ( "  RD returned IERR = %d\n", ierr );
    exit ( 1 );
  }

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

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

double elliptic_inc_fm ( double phi, double m )

/******************************************************************************/
/*
  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:

    05 January 2022

  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 )
  {
    printf ( "\n" );
    printf ( "elliptic_inc_fm(): Fatal error!\n" ); 
    printf ( "  RF returned IERR = %d\n", ierr );
    exit ( 1 );
  }

  value = sp * value;

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

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

/******************************************************************************/
/*
  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 incomplete 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 )
  {
    printf ( "\n" );
    printf ( "rd(): Fatal error!\n" );
    printf ( "  Invalid input arguments.\n" );
    printf ( "  X = %g\n", x );
    printf ( "  Y = %g\n", y );
    printf ( "  Z = %g\n", z );
    printf ( "\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;
  }

}
/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    05 January 2023

  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 incomplete 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 )
  {
    printf ( "\n" );
    printf ( "rf(): Fatal error!\n" );
    printf ( "  Invalid input arguments.\n" );
    printf ( "  X = %g\n", x );
    printf ( "  Y = %g\n", y );
    printf ( "  Z = %g\n", z );
    printf ( "\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;
  }

}