esei_OblateSpheroid.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                         ###  ###  #####  ####  ##  ##
//                       ##  ##  ########  ##    ##  ## ##  ##
//                       ##  ##  ## ## ##  ###   ##     ######
//                       ##  ##  ##    ##  ##    ##  ## ##  ##
//                       #####   ##    ##  #####  ####  ##  ##
//                       ##
//
// name:           eshelbySoluEllipticIntegralsOblateSpheroid.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 "eshelbySoluLambda.h"
#include "legendreIntegrals.h"
#include "elementaryFunctions.h"
#include "esei_OblateSpheroid.h"
#include "inclusion.h"


namespace mumech {


// --------------------------------------------------------------------------
// PUBLIC FUNCTIONS

// Function gives the values of Ferers-Dysons integrals
// Last edit:   11-Sep-2013
void eshelbySoluEllipticIntegralsOblateSpheroid::giveEllipticIntegrals(double J[13], double lambda, bool intpoint)
{
#ifdef DEBUG
  // Oblate spheroid: a1 = a2 > a3
  if (I->a[0] != I->a[1])  _errorr( "Values of semiaxes do not fit the condition for oblate spheroid (a1 = a2). " );
  if (I->a[0] <= I->a[2])  _errorr( "Values of semiaxes do not fit the condition for oblate spheroid (a1 > a3)." );
#endif
  
  double a1 = I->a[0];
  double a3 = I->a[2];
  
  J[0] = 0.;
  
  J[1] = this->I1(a1, a3, lambda, intpoint); // I1
  J[2] = this->I2(J[1]); // I2
  J[3] = this->I3(a1, a3, J[1], lambda, intpoint); // 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(a1, a3, J[6], lambda, intpoint); // I12
  J[4] = this->I11(J[5]); // I11

  J[7] = J[5]; // I21
  J[8] = this->I22(J[5]); // I22
  J[9] = this->I23(J[6]); // I23

  J[10] = J[6]; // I31
  J[11] = J[9]; // I32
  J[12] = this->I33(a1, a3, J[6], lambda, intpoint); // I33

  //#ifdef MY_DEBUG
  //    std::cout << "Ii, Iij = ";
  //    for (int i = 0; i < 13; i++)
  //      std::cout << J[i] << ", ";
  //    std::cout << std::endl;
  //#endif
}

// Function gives the derivatives of Ferers-Dysons integrals
// Last edit:   11-Sep-2013
void eshelbySoluEllipticIntegralsOblateSpheroid::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);
    }
  
}

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

// Function gives the solution of the integral I
// last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsOblateSpheroid::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 eshelbySoluEllipticIntegralsOblateSpheroid::I1(double a1, double a3, double lambda, bool intpoint)
{
  if(intpoint) // Internal points
    {
      // Mura (11.28)
      return 2.0 * PI * SQR(a1) * a3 / pow(SQR(a1) - SQR(a3), 3.0 / 2.0) * (acos(a3 / a1) - a3 / a1 * sqrt(1.0 - SQR(a3 / a1)));
    }
  else // External points
    {
      // Mura (12.21)
      double b = get_b(a1, a3, lambda);
      double d = get_d(a1, a3, lambda);
      return 2.0 * PI * SQR(a1) * a3 / pow(SQR(a1) - SQR(a3), 3.0 / 2.0) * (acos(b) - b * d);
    }
}

// Function gives the solution of the integral I_2
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsOblateSpheroid::I2(double I1)
{
  // Mura (11.28), (12.21)
  return I1;
}

// Function gives the solution of the integral I_3
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsOblateSpheroid::I3(double a1, double a3, double I1, double lambda, bool intpoint)
{
  if(intpoint) // Internal points
    {
      // Mura (11.28)
      return _4PI_ - 2.0 * I1;
    }
  else // External points
    {
      // Mura (12.21)
      double b = get_b(a1, a3, lambda);
      double d = get_d(a1, a3, lambda);
      return _4PI_ * SQR(a1) * a3 * (d / b - acos(b)) / pow(SQR(a1) - SQR(a3), 3.0 / 2.0);
    }
}

// Function gives the solution of the integral I_11
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsOblateSpheroid::I11(double I12)
{
  // Mura (11.28), also applies to external field.
  return I12;
}

// Function gives the solution of the integral I_12
// Last edit:   11-Sep-2013
double eshelbySoluEllipticIntegralsOblateSpheroid::I12(double a1, double a3, double I13, double lambda, bool intpoint)
{
  if(intpoint) // Internal points
    {
      // Mura (11.28)
      return PI / SQR(a1) - 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_ * SQR(a1) * a3 / (SQR(a1) + lambda) / delta_lambda - I13);
    }
}

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

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

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

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

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

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

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

// Lambda computed according to Mura p. 86. The property that a1 == a2 is taken into account to simplify the computation.
double eshelbySoluEllipticIntegralsOblateSpheroid::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(x2)) * SQR(a3) - SQR(x3) * SQR(a1);

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


} // end of namespace mumech

/*end of file*/