legendreIntegrals.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           legendreIntegrals.cpp
// author(s):      Jan Novak
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "legendreIntegrals.h"
#include "gelib.h"
#include "macros.h"


namespace mumech {

//********************************************************************************************************
//                                            PUBLIC FUNCTIONS
//********************************************************************************************************

//********************************************************************************************************
// description: Function calculating the Legendre's integral of the first kind
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Legendre elliptic integral of the 1st kind F(@Theta, @k), evaluated using Carlson's function RF.
// The argument ranges are 0 <= Theta <= PI/2, 0 <= k*sin(Theta) <= 1.
//********************************************************************************************************
double legendreIntegrals :: ellf (double phi, double ak)
{
  double s = sin (phi);
  
  return ( s * legendreIntegrals::rf( SQR( cos( phi ) ), ( 1.0 - s * ak ) * ( 1.0 + s * ak ), 1.0 ) );
}

//********************************************************************************************************
// description: Function calculating the Legendre's integral of the seccond kind
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Legendre elliptic integral of the 2nd kind E(@Theta, @k), evaluated using Carlson's functions RD and RF.
// The argument ranges are 0 <= Theta <= PI/2, 0 <= k*sin(Theta) <= 1.
//********************************************************************************************************
double legendreIntegrals :: elle (double phi, double ak)
{
  double s, cc, q;
  s = sin( phi );
  cc = SQR( cos( phi ) );
  q = ( 1.0 - s * ak ) * ( 1.0 + s * ak );
  
  return ( s * ( legendreIntegrals::rf( cc, q, 1.0 ) - ( SQR( s * ak ) ) * legendreIntegrals::rd( cc, q, 1.0 )/3.0 ) );
}

//********************************************************************************************************
// description: Function calculating the Legendre's integral of the third kind
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Legendre elliptic integral of the 3rd kind PI(@Theta, @n, @k), evaluated using Carlson's functions RJ and
// RF. (Note that the sign convention on n is opposite that of Abramowitz and Stegun.)
// The argument ranges are 0 <= Theta <= PI/2, 0 <= k*sin(Theta) <= 1.
//********************************************************************************************************
double legendreIntegrals :: ellpi (double phi, double en, double ak)
{
  double s, enss, cc, q;
  s = sin( phi );
  enss = en * s * s;
  cc = SQR( cos( phi ) );
  q = ( 1.0 - s * ak ) * ( 1.0 + s * ak );
  
  return ( s * ( legendreIntegrals::rf( cc, q, 1.0 ) - enss * legendreIntegrals::rj( cc, q, 1.0, 1.0 + enss )/3.0 ) );
}


//********************************************************************************************************
//                                            PROTECTED FUNCTIONS
//********************************************************************************************************

//********************************************************************************************************
// description: Solution of Carlson's integral of the first kind - numerical recepies
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Computes Carlson's elliptic integral of the first kind, RF(@x, @y, @z). @x, @y, and @z must be
// nonnegative, and at most one can be zero. TINY must be at least 5 times the machine underflow limit,
// BIG at most one fifth the machine overflow limit.
//********************************************************************************************************

#define RF_ERRTOL 0.0025
//#define RF_TINY 1.5e-38  //as in Numer. recepies
//#define RF_BIG 3.0e37    //as in Numer. recepies
#define RF_TINY ( DBL_MIN * 5. )
#define RF_BIG  ( DBL_MAX / 5. )
#define RF_THIRD ( 1.0/3.0 )
#define RF_C1 ( 1.0/24.0 )
#define RF_C2 0.1
#define RF_C3 ( 3.0/44.0 )
#define RF_C4 ( 1.0/14.0 )

double legendreIntegrals::rf( double x, double y, double z )
{
  double alamb, ave, delx, dely, delz, e2, e3, sqrtx, sqrty, sqrtz, xt, yt, zt;

  if ( MIN( MIN( x, y ), z ) < 0.0 ||
       MIN( MIN( x + y, x + z ), y + z ) < RF_TINY ||
       MAX( MAX( x, y ), z ) > RF_BIG ){
    _errorr1( "legendreIntegrals::rf: invalid argument(s)\n" );
  }
  
  xt = x;
  yt = y;
  zt = z;
  
  do {
    sqrtx = sqrt( xt );
    sqrty = sqrt( yt );
    sqrtz = sqrt( zt );
    alamb = sqrtx * ( sqrty + sqrtz ) + sqrty * sqrtz;
    xt = 0.25 * ( xt + alamb );
    yt = 0.25 * ( yt + alamb );
    zt = 0.25 * ( zt + alamb );
    ave = RF_THIRD * ( xt + yt + zt );
    delx = ( ave - xt )/ave;
    dely = ( ave - yt )/ave;
    delz = ( ave - zt )/ave;
  }while( MAX( MAX( ABS( delx ), ABS( dely ) ), ABS( delz ) ) > RF_ERRTOL );
  
  e2 = delx * dely - delz * delz;
  e3 = delx * dely * delz;
  
  return( 1.0 + ( RF_C1 * e2 - RF_C2 - RF_C3 * e3 ) * e2 + RF_C4 * e3 )/sqrt( ave );
}//end of function: rf()

//********************************************************************************************************
// description: Solution of Carlson's integral of the seccond kind - numerical recepies
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Computes Carlson's elliptic integral of the second kind, RD(@x, @y, @z). @x and @y must be nonnegative,
// and at most one can be zero. @z must be positive. TINY must be at least twice the
// negative 2/3 power of the machine overflow limit. BIG must be at most 0.1*ERRTOL times
// the negative 2/3 power of the machine underflow limit.
//********************************************************************************************************

#define RD_ERRTOL 0.0015
//#define RD_TINY 1.0e-25  //as in Numer. recepies
//#define RD_BIG 4.5e21    //as in Numer. recepies
#define RD_TINY ( 2. * pow( DBL_MAX, (-2./3.) ) )
#define RD_BIG  ( 0.1 * RD_ERRTOL * pow( DBL_MIN, (-2./3.) ) )
#define RD_C1 ( 3.0/14.0 )
#define RD_C2 ( 1.0/6.0 )
#define RD_C3 ( 9.0/22.0 )
#define RD_C4 ( 3.0/26.0 )
#define RD_C5 ( 0.25 * RD_C3 )
#define RD_C6 ( 1.5 * RD_C4 )

double legendreIntegrals::rd( double x, double y, double z )
{
  double alamb, ave, delx, dely, delz, ea, eb, ec, ed, ee, fac, sqrtx, sqrty, sqrtz, sum, xt, yt, zt;
  
  if( MIN( x, y ) < 0.0 || MIN( x + y, z ) < RD_TINY || MAX( MAX( x, y ), z ) > RD_BIG ){
    _errorr( "legendreIntegrals::rd: invalid argument(s)\n" );
  }
  
  xt = x;
  yt = y;
  zt = z;
  sum = 0.0;
  fac = 1.0;

  do{
    sqrtx = sqrt( xt );
    sqrty = sqrt( yt );
    sqrtz = sqrt( zt );
    alamb = sqrtx * ( sqrty + sqrtz ) + sqrty * sqrtz;
    sum += fac/( sqrtz *( zt + alamb ) );
    fac = 0.25 * fac;
    xt = 0.25 * ( xt + alamb );
    yt = 0.25 * ( yt + alamb );
    zt = 0.25 * ( zt + alamb );
    ave = 0.2 * ( xt + yt + 3.0 * zt );
    delx = ( ave - xt )/ave;
    dely = ( ave - yt )/ave;
    delz = ( ave - zt )/ave;
  }while( MAX( MAX( ABS( delx ), ABS( dely ) ), ABS( delz ) ) > RD_ERRTOL );

  ea = delx * dely;
  eb = delz * delz;
  ec = ea - eb;
  ed = ea - 6.0 * eb;
  ee = ed + ec + ec;
  
  return( 3.0 * sum + fac * ( 1.0 + ed * ( -RD_C1 + RD_C5 * ed - RD_C6 * delz * ee )
                  + delz * ( RD_C2 * ee + delz * ( -RD_C3 * ec + delz * RD_C4 * ea ) ) )/( ave * sqrt( ave ) ) );
}//end of function: rd()


//********************************************************************************************************
// description: Solution of Carlson's integral of the third kind - numerical recepies
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Computes Carlson's elliptic integral of the third kind, RJ(@x, @y, @z, @p). @x, @y, and @z must be
// nonnegative, and at most one can be zero. @p must be nonzero. If p < 0, the Cauchy principal
// value is returned. TINY must be at least twice the cube root of the machine underflow limit,
// BIG at most one fifth the cube root of the machine overflow limit.
//********************************************************************************************************

#define RJ_ERRTOL 0.0015
//#define RJ_TINY 2.5e-13  //as in Numer. recepies
//#define RJ_BIG 9.0e11    //as in Numer. recepies
#define RJ_TINY ( 2. * pow( DBL_MIN, ( 1./3. ) ) )
#define RJ_BIG  ( pow( DBL_MAX, ( 1./3. ) ) / 5. )
#define RJ_C1 ( 3.0/14.0 )
#define RJ_C2 ( 1.0/3.0 )
#define RJ_C3 ( 3.0/22.0 )
#define RJ_C4 ( 3.0/26.0 )
#define RJ_C5 ( 0.75 * RJ_C3 )
#define RJ_C6 ( 1.5 * RJ_C4 )
#define RJ_C7 ( 0.5 * RJ_C2 )
#define RJ_C8 ( RJ_C3 + RJ_C3 )

double legendreIntegrals::rj( double x, double y, double z, double p )
{
  double a = 0.0, alamb = 0.0, alpha = 0.0, ans = 0.0, ave = 0.0;
  double b = 0.0, beta = 0.0, delp = 0.0, delx = 0.0, dely = 0.0, delz = 0.0, ea = 0.0;
  double eb = 0.0, ec = 0.0, ed = 0.0, ee = 0.0, fac = 0.0, pt = 0.0, rcx = 0.0;
  double rho = 0.0, sqrtx = 0.0, sqrty = 0.0, sqrtz = 0.0, sum = 0.0, tau = 0.0;
  double xt = 0.0, yt = 0.0, zt = 0.0;
  /* Originally in Numerical recepies
     if ( MIN( MIN( x, y ), z ) < 0.0 ||
     MIN( MIN( x + y, x + z ), MIN( y + z, ABS( p ) ) ) < RJ_TINY ||
     MAX( MAX( x, y ), MAX( z, ABS( p ) ) ) > RJ_BIG ){
     _errorr( "legendreIntegrals::rj: invalid argument(s)\n" );
     }
  */
  //bug version from: http://www.nr.com/latest-known-bugs.html, last edit 02.09.2009
  if ( MIN( MIN( x, y ), z ) < 0.0 ||
       MIN( MIN( MIN( x + y, x + z ), y + z ), ABS( p ) ) < RJ_TINY ||
       MAX( MAX( MAX( x, y ), z ), ABS( p ) ) > RJ_BIG ){
    _errorr( "legendreIntegrals::rj: invalid argument(s)\n" );
  }
  
  sum = 0.0;
  fac = 1.0;
  
  if( p > 0.0 ){
    xt = x;
    yt = y;
    zt = z;
    pt = p;
  }
  else{
    xt = MIN( MIN( x, y ), z );
    zt = MAX( MAX( x, y ), z );
    yt = x + y + z - xt - zt;
    a = 1.0/( yt - p );
    b = a * ( zt - yt ) * ( yt - xt );
    pt = yt + b;
    rho = xt * zt/yt;
    tau = p * pt/yt;
    rcx = rc( rho, tau );
  }
  
  do{
    sqrtx = sqrt( xt );
    sqrty = sqrt( yt );
    sqrtz = sqrt( zt );
    alamb = sqrtx * ( sqrty + sqrtz ) + sqrty * sqrtz;
    alpha = SQR( pt * ( sqrtx + sqrty + sqrtz ) + sqrtx * sqrty * sqrtz );
    beta = pt * SQR( pt + alamb );
    sum += fac * rc( alpha, beta );
    fac = 0.25 * fac;
    xt = 0.25 * ( xt + alamb );
    yt = 0.25 *( yt + alamb );
    zt = 0.25 *( zt + alamb );
    pt = 0.25 *( pt + alamb );
    ave = 0.2 * ( xt + yt + zt + pt + pt );
    delx = ( ave - xt )/ave;
    dely = ( ave - yt )/ave;
    delz = ( ave - zt )/ave;
    delp = ( ave - pt )/ave;
  }
  /* Originally in Numerical recepies
     while( MAX( MAX( ABS( delx ), ABS( dely ) ), MAX( ABS( delz ), ABS( delp ) ) ) > RJ_ERRTOL );
  */
  //bug version from: http://www.nr.com/latest-known-bugs.html, last edit 02.09.2009
  while( MAX( MAX( MAX( ABS( delx ), ABS( dely ) ), ABS( delz ) ), ABS( delp ) ) > RJ_ERRTOL );
  
  ea = delx * ( dely + delz ) + dely * delz;
  eb = delx * dely * delz;
  ec = delp * delp;
  ed = ea - 3.0 * ec;
  ee = eb + 2.0 * delp * ( ea - ec );
  ans = 3.0 * sum + fac * ( 1.0 + ed * ( -RJ_C1 + RJ_C5 * ed - RJ_C6 * ee ) + eb * ( RJ_C7 + delp * ( -RJ_C8 + delp * RJ_C4 ) ) + delp * ea * ( RJ_C2 - delp * RJ_C3 ) -RJ_C2 * delp * ec )/( ave * sqrt( ave ) );
  
  if( p <= 0.0 ){
    ans = a * ( b * ans + 3.0 * ( rcx - rf( xt, yt, zt ) ) );
  }
  
  return( ans );
}//end of function: rj()

//********************************************************************************************************
// description: Solution of Carlson's degenerate elliptic integral - numerical recepies
// last edit:   24. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Computes Carlson's degenerate elliptic integral, RC(@x, @y). @x must be nonnegative and y must
// be nonzero. If y < 0, the Cauchy principal value is returned. TINY must be at least 5 times
// the machine underflow limit, BIG at most one fifth the machine maximum overflow limit.
//********************************************************************************************************

#define RC_ERRTOL 0.0012
//#define RC_TINY 1.69e-38    //as in Numer. recepies
//#define RC_SQRTNY 1.3e-19   //as in Numer. recepies
//#define RC_BIG 3.e37        //as in Numer. recepies
#define RC_TINY ( DBL_MIN * 5. )
#define RC_SQRTNY ( sqrt( RC_TINY ) )
#define RC_BIG ( DBL_MAX / 5. )
#define RC_TNBG ( RC_TINY * RC_BIG )
#define RC_COMP1 ( 2.236/RC_SQRTNY )
#define RC_COMP2 ( RC_TNBG * RC_TNBG/25.0 )
#define RC_THIRD ( 1.0/3.0)
#define RC_C1 0.3
#define RC_C2 ( 1.0/7.0 )
#define RC_C3 0.375
#define RC_C4 ( 9.0/22.0 )

double legendreIntegrals::rc( double x, double y )
{
  double alamb, ave, s, w, xt, yt;
  if ( x < 0.0 || y == 0.0 || ( x + ABS( y ) ) < RC_TINY ||
       ( x + ABS( y ) ) > RC_BIG || ( y <- RC_COMP1 && x > 0.0 && x < RC_COMP2 ) ){
    _errorr( "legendreIntegrals::rc: invalid argument(s)\n" );
  }
  
  if( y > 0.0 ){
    xt = x;
    yt = y;
    w = 1.0;
  }
  else{
    xt = x - y;
    yt = -y;
    w = sqrt( x )/sqrt( xt );
  }
  
  do{
    alamb = 2.0 * sqrt( xt ) * sqrt( yt ) + yt;
    xt = 0.25 * ( xt + alamb );
    yt = 0.25 * ( yt + alamb );
    ave = RC_THIRD * ( xt + yt + yt );
    s = ( yt - ave )/ave;
  }while( fabs( s ) > RC_ERRTOL );
  
  return( w * ( 1.0 + s * s * ( RC_C1 + s * ( RC_C2 + s * ( RC_C3 + s * RC_C4 ) ) ) )/sqrt( ave ) );
}//end of function: rc()


} // end of namespace mumech

/*end of file*/