esei_ProlateSpheroid.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluEllipticIntegralsProlateSpheroid.cpp
// author(s):      Jan Novak, Tomas Janda
// 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 <cmath>
#include <iostream>

#include "macros.h"
#include "esei_ProlateSpheroid.h"
#include "eshelbySoluLambda.h"
#include "legendreIntegrals.h"
#include "elementaryFunctions.h"
#include "inclusion.h"


namespace mumech {
 
// Function gives the values of Ferers-Dysons integrals
// Last edit:   11-Sep-2013
void eshelbySoluEllipticIntegralsProlateSpheroid::giveEllipticIntegrals(double J[13], double lambda, bool intpoint)
{
#ifdef DEBUG
  if (! (I->a[0] > I->a [1] && I->a[1] == I->a[2] && I->a[2] > 0.0))  _errorr("Semiaxes are not OK.");
#endif
  
  double a1 = I->a[0];
  double a3 = I->a[2];
  
  J[0] = 0.0;
  
  J[2] = this->I2(a1, a3, lambda, intpoint); // I2
  J[1] = this->I1(a1, a3, J[2], lambda, intpoint); // I1
  J[3] = this->I3(J[2]); // I3
  
  // (Order changed due to sequential evaluation I13 -> I12 -> I11)
  J[6] = this->I13(a1, a3, J[1], J[3]); // I13
  J[5] = this->I12(J[6]); // I12
  J[4] = this->I11(a1, a3, J[6], lambda, intpoint); // I11
  
  double I33 = this->I33(a1, a3, J[6], lambda, intpoint);
  
  J[7] = J[5]; // I21
  J[8] = this->I22(I33); // I22
  J[9] = this->I23(I33); // I23
  
  J[10] = J[6]; // I31
  J[11] = J[9]; // I32
  J[12] = I33; // I33
}

// Function gives the derivatives of Ferers-Dysons integrals
// Last edit:   11-Sep-2013
void eshelbySoluEllipticIntegralsProlateSpheroid::giveDerivativesOfEllipticIntegrals(Point *point, bool intpoint)
{
  if (intpoint) _errorr("derivatives are 0, this function should not be called");
  
  double ndiff = I->ndiff_1;

  // NUMERICAL DERIVATIVES (using REAL ARGUMENTS)
  
  // Get perturbated points and their lambdas
  double lambdas[NUM_PERTURB];
  this->getPerturbatedLambdas(lambdas, point->loc_x);
  
  // Create and fill a two dimensional array containing values of J-integrals at perturbated points.
  const int NUM_ELLIP_INT = 13;
  double J_pert[NUM_PERTURB][NUM_ELLIP_INT]; // All perturbations of all elliptical integrals.
  for (int i = 0; i < NUM_PERTURB; i++) // For each perturbated point
    this->giveEllipticIntegrals(J_pert[i], lambdas[i], intpoint); // Get elliptical integrals for lambda of perturbated point.
  
  // Compute numerial derivatives of first and second order of J_i integrals.
  for (int i = 0; i < 3; i++) // For each elliptic integral Ji
    {
      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Ji_j = dJi[3*i+j-4], i,j=1..3.
      point->dJi[3*i] = (J_pert[_Xp_][i+1] - J_pert[_Xm_][i+1]) / (2.0 * ndiff);
      point->dJi[3*i+1] = (J_pert[_Yp_][i+1] - J_pert[_Ym_][i+1]) / (2.0 * ndiff);
      point->dJi[3*i+2] = (J_pert[_Zp_][i+1] - J_pert[_Zm_][i+1]) / (2.0 * ndiff);
      
      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Ji_jk = ddJi[9*i+3*j+k-13], i,j,k=1..3.
      point->ddJi[9*i] =   (J_pert[_Xp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Xm_][i+1]) / (ndiff * ndiff); // d^2f(x,y,z)/dx^2 = (f(x+h,y,z) - 2*f(x,y,z) + f(x-h,y,z)) / h^2
      point->ddJi[9*i+1] = (J_pert[_XpYp_][i+1] - J_pert[_XpYm_][i+1] - J_pert[_XmYp_][i+1] + J_pert[_XmYm_][i+1]) / (4.0*ndiff*ndiff); // d^2f(x,y,z)/dxdy = (f(x+h,y+h,z) - f(x+h,y-h,z) - f(x-h,y+h,z) + f(x-h,y-h,z)) / (4*h^2)
      point->ddJi[9*i+2] = (J_pert[_XpZp_][i+1] - J_pert[_XpZm_][i+1] - J_pert[_XmZp_][i+1] + J_pert[_XmZm_][i+1]) / (4.0*ndiff*ndiff);
      
      point->ddJi[9*i+3] = point->ddJi[9*i+1];
      point->ddJi[9*i+4] = (J_pert[_Yp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Ym_][i+1]) / (ndiff * ndiff);
      point->ddJi[9*i+5] = (J_pert[_YpZp_][i+1] - J_pert[_YpZm_][i+1] - J_pert[_YmZp_][i+1] + J_pert[_YmZm_][i+1]) / (4.0*ndiff*ndiff);
      
      point->ddJi[9*i+6] = point->ddJi[9*i+2];
      point->ddJi[9*i+7] = point->ddJi[9*i+5];
      point->ddJi[9*i+8] = (J_pert[_Zp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Zm_][i+1]) / (ndiff * ndiff);
    }

  // Compute numerial derivatives of first and second order of J_ij integrals.
  for (int i = 0; i < 9; i++) // For each elliptic integral Jij
    {
      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Jij_k = dJij[9*k+3*i+j-13], i,j,k=1..3.
      point->dJij[i] =  (J_pert[_Xp_][i+4] - J_pert[_Xm_][i+4]) / (2.0 * ndiff);
      point->dJij[i+9] = (J_pert[_Yp_][i+4] - J_pert[_Ym_][i+4]) / (2.0 * ndiff);
      point->dJij[i+18] = (J_pert[_Zp_][i+4] - J_pert[_Zm_][i+4]) / (2.0 * ndiff);

      // See types.c, section "Derivatives of Ferers-Dysons' integrals".
      // Jij_kl = ddJi[27*i+9*j+3*k+l-40] i,j,k,l=1..3.
      point->ddJij[9*i] =   (J_pert[_Xp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Xm_][i+4]) / (ndiff * ndiff);
      point->ddJij[9*i+1] = (J_pert[_XpYp_][i+4] - J_pert[_XpYm_][i+4] - J_pert[_XmYp_][i+4] + J_pert[_XmYm_][i+4]) / (4.0*ndiff*ndiff);
      point->ddJij[9*i+2] = (J_pert[_XpZp_][i+4] - J_pert[_XpZm_][i+4] - J_pert[_XmZp_][i+4] + J_pert[_XmZm_][i+4]) / (4.0*ndiff*ndiff);

      point->ddJij[9*i+3] = point->ddJij[9*i+1];
      point->ddJij[9*i+4] = (J_pert[_Yp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Ym_][i+4]) / (ndiff * ndiff);
      point->ddJij[9*i+5] = (J_pert[_YpZp_][i+4] - J_pert[_YpZm_][i+4] - J_pert[_YmZp_][i+4] + J_pert[_YmZm_][i+4]) / (4.0*ndiff*ndiff);

      point->ddJij[9*i+6] = point->ddJij[9*i+2];
      point->ddJij[9*i+7] = point->ddJij[9*i+5];
      point->ddJij[9*i+8] = (J_pert[_Zp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Zm_][i+4]) / (ndiff * ndiff);
    }
  //}
  //else _errorr("point position is not defined");
}


// --------------------------------------------------------------------------
// PROTECTED FUNCTIONS

// Function gives the solution of the integral I
// last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::eI( double a[3], double k, double Theta, double mult )
{
  double a1 = a[0], a3 = a[2];
  legendreIntegrals ellInt;
  
  return( mult * ellInt.ellf( Theta, k ) / sqrt( SQR( a1 ) - SQR( a3 ) ) );
}

// Function gives the solution of the integral I_1
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I1(double a1, double a3, double I2, double lambda, bool intpoint)
{
  if (intpoint) // Internal points
    {
      // Mura (11.29)
      return _4PI_ - 2.0 * I2;
    }
  else // External points
    {
      // Mura (12.22)
      double b = get_b(a1, a3, lambda);
      double d = get_d(a1, a3, lambda);
      return _4PI_ * a1 * SQR(a3) * (acosh(b) - d / b) / pow(SQR(a1) - SQR(a3), 3.0 / 2.0);
    }
}

// Function gives the solution of the integral I_2
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I2(double a1, double a3, double lambda, bool intpoint)
{
  if (intpoint) // Internal points
    {
      // Mura (11.29)
      return 2.0 * PI * a1 * SQR(a3) / pow(SQR(a1) - SQR(a3), 3.0 / 2.0) * (a1 / a3 * sqrt(SQR(a1 / a3) - 1.0) - acosh(a1 / a3));
    }
  else // External points
    {
      // Mura (12.22)
      double b = get_b(a1, a3, lambda);
      double d = get_d(a1, a3, lambda);
      return 2.0 * PI * a1 * SQR(a3) / pow(SQR(a1) - SQR(a3), 3.0 / 2.0) * (b * d - acosh(b));
    }
}

// Function gives the solution of the integral I_3
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I3(double I2)
{
  return I2;
}

// Function gives the solution of the integral I_11
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I11(double a1, double a3, double I13, double lambda, bool intpoint)
{
  if (intpoint) // Internal points
    {
      // Mura (11.29)
      return 1.0 / 3.0 * (_4PI_ / SQR(a1) - 2.0 * I13);
    }
  else // External points
    {
      // Mura (12.18)
      double delta_lambda = get_delta_lambda(a1, a3, lambda);
      return 1.0 / 3.0 * (_4PI_ * a1 * SQR(a3) / (SQR(a1) + lambda) / delta_lambda - 2.0 * I13);
    }
}

// Function gives the solution of the integral I_12
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I12(double I13)
{
  // Mura (11.29), also applies to external field.
  return I13;
}

// Function gives the solution of the integral I_13
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I13( double a1, double a3, double I1, double I3 )
{
  // Mura (11.29), (12.18)
  return (I1 - I3)/(SQR(a3) - SQR(a1));
}

// Function gives the solution of the integral I_22
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I22(double I33)
{
  // Mura (11.29), also applies to external field.
  return I33;
}

// Function gives the solution of the integral I_23
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I23(double I33)
{
  return I33;
}

// Function gives the solution of the integral I_33
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsProlateSpheroid::I33(double a1, double a3, double I13, double lambda, bool intpoint)
{
  if (intpoint) // Internal points
    {
      // Mura (11.29)
      return PI / SQR(a3) - I13 / 4.0;
    }
  else // External points
    {
      double delta_lambda = get_delta_lambda(a1, a3, lambda);
      // Derived from third line in (12.18) and I11 = I22 = I12
      return 1.0 / 4.0 * (_4PI_ * a1 * SQR(a3) / (SQR(a3) + lambda) / delta_lambda - I13);
    }
}

// Auxiliary value b, (Mura, p.94)
double eshelbySoluEllipticIntegralsProlateSpheroid::get_b(double a1, double a3, double lambda)
{
  return sqrt((SQR(a1) + lambda) / (SQR(a3) + lambda));
}

// Auxiliary value d, (Mura, p.94)
double eshelbySoluEllipticIntegralsProlateSpheroid::get_d(double a1, double a3, double lambda)
{
  return sqrt((SQR(a1) - SQR(a3)) / (SQR(a3) + lambda));
}

// Auxiliary value Delta_lambda, (Mura, p.93)
double eshelbySoluEllipticIntegralsProlateSpheroid::get_delta_lambda(double a1, double a3, double lambda)
{
  return sqrt((SQR(a1) + lambda) * (SQR(a3) + lambda) * (SQR(a3) + lambda)); // Prolate spheroid: a2 == a3
}

// Lambda computed according to Mura p. 86. The property that a2 == a3 is taken into account to simplify the computation.
double eshelbySoluEllipticIntegralsProlateSpheroid::getLambda(const double a[3], double x1, double x2, double x3)
{
  double a1 = a[0];
  double a3 = a[2];

  // Coeficients of quadratic equation
  // a = 1
  double b = SQR(a1) + SQR(a3) - SQR(x1) - SQR(x2) - SQR(x3);
  double c = SQR(a1) * SQR(a3) - SQR(x1) * SQR(a3) - (SQR(x2) + SQR(x3)) * SQR(a1);

  return (-b + sqrt(SQR(b) - 4.0 * c)) / 2.0;
}



} // end of namespace mumech

/*end of file*/