esei_Sphere.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           eshelbySoluEllipticIntegralsSphere.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 <cmath>
#include <iostream>

#include "macros.h"
#include "esei_Sphere.h"
#include "elementaryFunctions.h"
#include "esei.h"
#include "problem.h"


namespace mumech {

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

  
void eshelbySoluEllipticIntegralsSphere::giveEllipticIntegrals(double J[13], double lambda, bool intpoint)
{
  double a0 = I->a[0];
  
  J[0] = 0.;
  if (intpoint)
    {
      J[1] = J[2] = J[3] = 4.0 * PI / 3.0;
      J[4] = J[5] = J[6] = J[7] = J[8] = J[9] = J[10] = J[11] = J[12] = 4.0 * PI / 5.0 / SQR(a0);
    }
  else
    {
      J[1] = J[2] = J[3] = 4.0 * PI * a0*a0*a0 / 3.0 / pow(a0*a0 + lambda, 3.0 / 2.0);
      J[4] = J[5] = J[6] = J[7] = J[8] = J[9] = J[10] = J[11] = J[12] = 4.0 * PI * a0*a0*a0 / 5.0 / pow(a0*a0 + lambda, 5.0 / 2.0);
    }
}//end of function: giveEllipticIntegrals( )


// ****************************************************************************************
// start: Alternative, more time demanding, evaluation of integrals
// ****************************************************************************************

// termitovo: PROC JE TADY TATO FCE, DRIVE BYLA POUZIVANA, ALE JE TO ASI NEJAKY EXPERIMENT, NE?
//
// //********************************************************************************************************
// // description: Function gives the values of all elliptic integrals
// // last edit:   08. 06. 2010
// //-------------------------------------------------------------------------------------------------------
// // note:        @J[13]     - elliptic integral values
// //              @a[3]      - dimensions of sphere radii (a1=a2=a3)
// //              @lambda    - lambda parameter according Mura's book, page 73
// //              @loc_x[3]  - local coordinates of a point, in this case are the same as the global ones
// //********************************************************************************************************
// void eshelbySoluEllipticIntegralsSphere::giveEllipticIntegrals( double J[13], const double a[3], double lambda,
//  const double loc_x[3] )
// {
//  double aa = SQR( a[0] ), mult = 4. * PI * aa * a[0], xx = 0., aux_x = 0.,
//    n_x[5] = { 0., 0., 0., 0., 0. };
// 
//  if( lambda > MACHINE_EPS ){
//    xx = SQR( loc_x[0] ) + SQR( loc_x[1] ) + SQR( loc_x[2] ); aux_x = sqrt( xx );
//    n_x[0] = aux_x; n_x[1] = n_x[0] * xx; n_x[2] = n_x[1] * xx; n_x[3] = n_x[2] * xx; n_x[4] = n_x[3] * xx;
//  }else ;
// 
//  //Ferers-Dysons integrals:
//  //------------------------
//  this->give_Ii_and_Iij( J, mult, n_x, lambda, aa );
// 
//  return;
// }//end of function: giveEllipticIntegrals( )
// 
// //********************************************************************************************************
// // description: Function gives the values of all elliptic Ferers-Dysons integrals
// // last edit:   08. 05. 2010
// //-------------------------------------------------------------------------------------------------------
// // note:        @J[13]     - table of elliptic integral values to be calculated
// //              @mult      - 4 * pi() * a^3
// //              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
// //              @lambda    - lambda parameter
// //              @aa        - square power of a[i] radius
// //********************************************************************************************************
// void eshelbySoluEllipticIntegralsSphere::give_Ii_and_Iij( double J[13], double mult, double norm_x[5],
//                            double lambda, double aa , bool intpoint )
// {
//   //J[0] = this->eI( mult, norm_x, a, lambda ); //this value is not necessary anywhere else
//   J[0] = 0.;
//   J[1] = J[2] = J[3] = this->Ii( mult, norm_x, lambda, intpoint);
//   J[4] = J[5] = J[6] = J[7] = J[8] = J[9] = J[10] = J[11] = J[12] = this->Iij( mult, norm_x, lambda, aa );
// 
//   return;
// }//end of function: give_Ii_and_Iij( )
// 
// //********************************************************************************************************
// // description: Function gives the solution of the integral I_1
// // last edit:   05. 05. 2010
// //-------------------------------------------------------------------------------------------------------
// // note:        I1 = I2 = I3
// //              @mult      - 4 * pi() * a^3
// //              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
// //              @a         - sphere radii
// //              @lambda    - lambda parameter
// //********************************************************************************************************
// double eshelbySoluEllipticIntegralsSphere::Ii( double mult, double norm_x[5], double lambda, bool intpoint )
// {
//   return ( intpoint  ?  4. * PI / 3.  :  mult/( 3. * norm_x[1] ) );
// }
// 
// //********************************************************************************************************
// // description: Function gives the solution of the integral I_11
// // last edit:   05. 05. 2010
// //-------------------------------------------------------------------------------------------------------
// // note:        I11 = I12 = I13 = I22 = I23 = I33
// //              @mult      - 4 * pi() * a^3
// //              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
// //              @a[3]      - sphere radii
// //              @lambda    - lambda parameter
// //              @aa        - square power of a[i] radius
// //********************************************************************************************************
// double eshelbySoluEllipticIntegralsSphere::Iij( double mult, double norm_x[5], double lambda, double aa, bool intpoint)
// {
//   return ( intpoint  ?  4. * PI / ( 5. * aa )  :  mult/( 5. * norm_x[2] ) );
// }
// //********************************************************************************************************
// // description: Function gives the solution of the integral I
// // last edit:   05. 05. 2010
// //-------------------------------------------------------------------------------------------------------
// // note:        @mult      - 4 * pi() * a^3
// //              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
// //              @a[3]      - sphere radii
// //              @lambda    - lambda parameter
// //********************************************************************************************************
// double eshelbySoluEllipticIntegralsSphere::eI( double mult, double norm_x[5], const double a[3], double lambda, bool intpoint)
// {
//   return ( intpoint  ? mult/a[0] : mult/norm_x[0] );
// }
// 
// ****************************************************************************************
// finish: Alternative, more time demanding, evaluation of integrals
// ****************************************************************************************



//********************************************************************************************************
// description: Function gives the values of all elliptic integrals' derivatives
// last edit:   08. 06. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @point     - point record (data structure of a given point)
//              @a[3]      - dimensions of sphere semiaxes (radii) (a1=a2=a3)
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::giveDerivativesOfEllipticIntegrals(Point *point, bool intpoint)
{
  if (intpoint) _errorr("derivatives are 0, this function should not be called");
  
  double ndiff = I->ndiff_1;

  double a0 = I->a[0];
    
  double aa = SQR( a0 ), mult = 4. * PI * aa * a0, xx = 0., aux_x = 0.,
    n_x[5] = { 0., 0., 0., 0., 0. };
  
  // termitovo: toto je i v alternativni (asi pomaly) verzi vypoctu integralu, takze se to muze vykostit i odtud mozna
  //            a hlavne by se tady lambda vubec nemello pocitat, stejne jako v jinych tvarech
  //if( point->la > MACHINE_EPS )
  //{
      xx = SQR( point->loc_x[0] ) + SQR( point->loc_x[1] ) + SQR( point->loc_x[2] ); aux_x = sqrt( xx );
      n_x[0] = aux_x; n_x[1] = n_x[0] * xx; n_x[2] = n_x[1] * xx; n_x[3] = n_x[2] * xx; n_x[4] = n_x[3] * xx;
      //  }
      //else
      //errol; // test, sem by to nemelo dojit, a point->la by nemel byt lambda ???
  
  
  //point->position == _EXT_POINT_
  switch (I->P->give_diffType()){
  case DT_COMPLEX : {
    std::complex<double> cdJi[13], cdJij[27];
      
    //first derivatives Ii,i:
    this->give_cdIi( cdJi, mult, point->loc_x );
    //first derivatives Iij,i:
    this->give_cdIij( cdJij, mult, point->loc_x );
    //second derivatives Ii,ij:
    this->give_cddIi( point->ddJi, cdJi );
    //second derivatives Iij,ij:
    this->give_cddIij( point->ddJij, cdJij );
      
    //complex-to-real conversion
    for( int i = 0; i < 13; i++ ) point->dJi[i] = real( cdJi[i] );
    for( int i = 0; i < 27; i++ ) point->dJij[i] = real( cdJij[i] );
    break;
  }
  case DT_NUMERICAL : { // Real numerical derivatives
    // 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);
      }
    break;
  }
  case DT_ANALITICAL :// Analytical derivatives
    {
      //first derivatives Ii,i:
      this->give_dIi( point->dJi, mult, n_x, point->loc_x );
      //first derivatives Iij,i:
      this->give_dIij( point->dJij, mult, n_x, point->loc_x );
      //second derivatives Ii,ij:
      this->give_ddIi( point->ddJi, mult, n_x, point->loc_x );
      //second derivatives Iij,ij:
      this->give_ddIij( point->ddJij, mult, n_x, point->loc_x );
      break;
    }
  default:
    _errorr( "Invalid type of differentiation.\n");
  }
}//end of function: giveDerivativesOfEllipticIntegrals( )
  
//********************************************************************************************************
//                                            PROTECTED FUNCTIONS
//********************************************************************************************************

//Ferers-Dysons integrals:
//------------------------




//Analytical derivatives of Ferers-Dysons integrals:
//--------------------------------------------------

//********************************************************************************************************
// description: Function gives the first derivative of the integral Ii
// last edit:   05. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        I1_i = I2_i = I3_i
//              @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//              @direction - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
double eshelbySoluEllipticIntegralsSphere::Ii_i( double mult, double norm_x[5], double x[3],
                         derivativeDirection direction )
{
  int dir = direction - 1;

  return( -mult * x[dir]/norm_x[2] );
}//end of function: Ii_i()

//********************************************************************************************************
// description: Function gives the second derivative of the integral Ii
// last edit:   05. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        I1_ij = I2_ij = I3_ij
//              @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//              @direction - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
double eshelbySoluEllipticIntegralsSphere::Ii_ij( double mult, double norm_x[5], double x[3],
                          derivativeDirection direction_1,
                          derivativeDirection direction_2 )
{
  int dir_1 = direction_1 - 1, dir_2 = direction_2 - 1;
  double aux = 5. * mult * x[dir_1] * x[dir_2]/norm_x[3];

  return ( dir_1 != dir_2 )? aux : ( aux - mult/norm_x[2] );
}//end of function: Ii_ij()

//********************************************************************************************************
// description: Function gives the first derivative of the integral Iij
// last edit:   05. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//              @direction - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
double eshelbySoluEllipticIntegralsSphere::Iij_i( double mult, double norm_x[5], double x[3],
                          derivativeDirection direction )
{
  int dir = direction - 1;

  return( -mult * x[dir]/norm_x[3] );
}//end of function: Iij_i()

//********************************************************************************************************
// description: Function gives the second derivative of the integral Iij
// last edit:   08. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//              @direction - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
double eshelbySoluEllipticIntegralsSphere::Iij_ij( double mult, double norm_x[5], double x[3],
                           derivativeDirection direction_1,
                           derivativeDirection direction_2 )
{
  int dir_1 = direction_1 - 1, dir_2 = direction_2 - 1;
  double aux = 7. * mult * x[dir_1] * x[dir_2]/norm_x[4];

  return ( dir_1 != dir_2 )? aux : ( aux - mult/norm_x[3] );
}//end of function: Ii_ij()


//********************************************************************************************************
// description: Function gives all first derivatives of Ferers-Dysons' integrals Ii
// last edit:   08. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @dJi[9]    - table of elliptic integral derivatives to be calculated
//              @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_dIi( double dJi[13], double mult, double norm_x[5],
                           double x[3] )
{
  //dIi,1
  dJi[0] = dJi[3] = dJi[6] = this->Ii_i( mult, norm_x, x, _x1_ ); //I1,1; I2,1; I3,1
  //dIi,2
  dJi[1] = dJi[4] = dJi[7] = this->Ii_i( mult, norm_x, x, _x2_ ); //I1,2; I2,2; I3,2
  //dIi,3
  dJi[2] = dJi[5] = dJi[8] = this->Ii_i( mult, norm_x, x, _x3_ ); //I1,3; I2,3; I3,3

  return;
}//end of function: give_dIi( )

//********************************************************************************************************
// description: Function gives all first derivatives of Ferers-Dysons' integrals Iij
// last edit:   08. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @dJij[27]  - table of elliptic integral derivatives to be calculated
//              @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_dIij( double dJij[27], double mult, double norm_x[5],
                            double x[3] )
{
  //dIij,1
  dJij[0]  = dJij[1] = dJij[2] = this->Iij_i( mult, norm_x, x, _x1_ ); //I11,1; I12,1; I13,1
  dJij[3]  = dJij[4] = dJij[5] = dJij[0]; //I21,1; I22,1; I23,1
  dJij[6]  = dJij[7] = dJij[8] = dJij[0]; //I31,1; I32,1; I33,1
  //dIij,2
  dJij[9]  = dJij[10] = dJij[11] = this->Iij_i( mult, norm_x, x, _x2_ ); //I11,2; I12,2; I13,2
  dJij[12] = dJij[13] = dJij[14] = dJij[9]; //I21,2; I22,2; I23,2
  dJij[15] = dJij[16] = dJij[17] = dJij[9]; //I31,2; I32,2; I33,2
  //dIij,3
  dJij[18] = dJij[19] = dJij[20] = this->Iij_i( mult, norm_x, x, _x3_ ); //I11,3; I12,3; I13,3
  dJij[21] = dJij[22] = dJij[23] = dJij[18]; //I21,3; I22,3; I23,3
  dJij[24] = dJij[25] = dJij[26] = dJij[18]; //I31,3; I32,3; I33,3

  return;
}//end of function: give_dIij( )

//********************************************************************************************************
// description: Function gives all second derivatives of Ferers-Dysons' integrals Ii
// last edit:   08. 05. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @ddJi[27]  - table of elliptic integral derivatives to be calculated
//              @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_ddIi( double ddJi[27], double mult, double norm_x[5],
                            double x[3] )
{
  //dI1,ij
  ddJi[0] = this->Ii_ij( mult, norm_x, x, _x1_, _x1_ ); //I1,11
  ddJi[1] = this->Ii_ij( mult, norm_x, x, _x1_, _x2_ ); //I1,12
  ddJi[2] = this->Ii_ij( mult, norm_x, x, _x1_, _x3_ ); //I1,13

  ddJi[3] = ddJi[1]; //I1,21
  ddJi[4] = this->Ii_ij( mult, norm_x, x, _x2_, _x2_ ); //I1,22
  ddJi[5] = this->Ii_ij( mult, norm_x, x, _x2_, _x3_ ); //I1,23

  ddJi[6] = ddJi[2]; //I1,31
  ddJi[7] = ddJi[5]; //I1,32
  ddJi[8] = this->Ii_ij( mult, norm_x, x, _x3_, _x3_ ); //I1,33

  //dI1,ij = dI2,ij = dI3,ij
  ddJi[9]  = ddJi[18] = ddJi[0]; //I2,11; I3,11
  ddJi[10] = ddJi[19] = ddJi[1]; //I2,12; I3,12
  ddJi[11] = ddJi[20] = ddJi[2]; //I2,13; I3,13

  ddJi[12] = ddJi[21] = ddJi[3]; //I2,21; I3,21
  ddJi[13] = ddJi[22] = ddJi[4]; //I2,22; I3,22
  ddJi[14] = ddJi[23] = ddJi[5]; //I2,23; I3,23

  ddJi[15] = ddJi[24] = ddJi[6]; //I2,31; I3,31
  ddJi[16] = ddJi[25] = ddJi[7]; //I2,32; I3,32
  ddJi[17] = ddJi[26] = ddJi[8]; //I2,33; I3,33

  return;
}//end of function: give_ddIi( )

//********************************************************************************************************
// description: Function gives all second derivatives of Ferers-Dysons' integrals Iij
// last edit:   15. 06. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @ddJij[81] - table of elliptic integral derivatives to be calculated
//              @mult      - 4 * pi() * a^3
//              @norm_x[5] - length of vector x and its powers norm_x = {x, x^3, x^5, x^7, x^9}
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_ddIij( double ddJij[81], double mult, double norm_x[5],
                             double x[3] )
{
  int j = 0;

  //dI11,ij
  ddJij[0] = this->Iij_ij( mult, norm_x, x, _x1_, _x1_ ); //I11,11
  ddJij[1] = this->Iij_ij( mult, norm_x, x, _x1_, _x2_ ); //I11,12
  ddJij[2] = this->Iij_ij( mult, norm_x, x, _x1_, _x3_ ); //I11,13

  ddJij[3] = ddJij[1]; //I11,21
  ddJij[4] = this->Iij_ij( mult, norm_x, x, _x2_, _x2_ ); //I11,22
  ddJij[5] = this->Iij_ij( mult, norm_x, x, _x2_, _x3_ ); //I11,23

  ddJij[6] = ddJij[2]; //I11,31
  ddJij[7] = ddJij[5]; //I11,32
  ddJij[8] = this->Iij_ij( mult, norm_x, x, _x3_, _x3_ ); //I11,33

  //dI11,ij = dI12,ij = dI13,ij = dI21,ij = dI22,ij = dI23,ij = dI31,ij = dI32,ij = dI33,ij
  for( int i = 1; i < 9; i++ ){
    j = i * 9;
    ddJij[j] = ddJij[0];
    ddJij[j + 1] = ddJij[1];
    ddJij[j + 2] = ddJij[2];

    ddJij[j + 3] = ddJij[3];
    ddJij[j + 4] = ddJij[4];
    ddJij[j + 5] = ddJij[5];

    ddJij[j + 6] = ddJij[6];
    ddJij[j + 7] = ddJij[7];
    ddJij[j + 8] = ddJij[8];
  }//end of loop 'for: i'

  return;
}//end of function: give_ddIij( )

//Complex numerical derivatives of Ferers-Dysons integrals:
//---------------------------------------------------------

//********************************************************************************************************
// description: Function gives the first derivative of the integral Ii by using
//              the complex differentiation
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        I1_i = I2_i = I3_i
//              @mult        - 4 * pi() * a^3
//              @x[3]        - coordinates of a point
//              @direction   - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
std::complex<double> eshelbySoluEllipticIntegralsSphere::cIi_i( double mult, std::complex<double> x[3],
                                derivativeDirection direction )
{
  int dir = direction - 1;
  std::complex<double> cn2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]);

  return( -mult * x[dir]/pow( cn2, 2.5 ) );
}//end of function: cIi_i()

//********************************************************************************************************
// description: Function gives the first derivative of the integral Iij by using
//              the complex differentiation
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @mult      - 4 * pi() * a^3
//              @x[3]      - coordinates of a point
//              @direction - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
std::complex<double> eshelbySoluEllipticIntegralsSphere::cIij_i( double mult, std::complex<double> x[3],
                                 derivativeDirection direction )
{
  int dir = direction - 1;
  std::complex<double> cn2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]);

  return( -mult * x[dir]/pow( cn2, 3.5 ) );
}//end of function: cIij_i()

//********************************************************************************************************
// description: Function gives all first derivatives of Ferers-Dysons' integrals Ii by using
//              the complex differentiation
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @dJi[9]    - table of elliptic integral derivatives to be calculated
//              @mult      - 4 * pi() * a^3
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_cdIi( std::complex<double> dJi[13], double mult, double x[3] )
{
  std::complex<double> cx[3] = { x[0], x[1], x[2] };

  this->update_cx( cx, _x1_ );
  dJi[0] = this->cIi_i( mult, cx, _x1_ ); //I1,1

  this->update_cx( cx, _x2_ );
  dJi[3] = this->cIi_i( mult, cx, _x1_ ); //I2,1
  dJi[1] = dJi[4] = this->cIi_i( mult, cx, _x2_ ); //I1,2; I2,2

  this->update_cx( cx, _x3_ );
  dJi[6] = this->cIi_i( mult, cx, _x1_ ); //I3,1
  dJi[7] = this->cIi_i( mult, cx, _x2_ ); //I3,2
  dJi[2] = dJi[5] = dJi[8] = this->cIi_i( mult, cx, _x3_ ); //I1,3; I2,3; I3,3

  return;
}//end of function: give_cdIi( )

//********************************************************************************************************
// description: Function gives all first derivatives of Ferers-Dysons' integrals Iij by using
//              the complex differentiation
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @dJij[27]  - table of elliptic integral derivatives to be calculated
//              @mult      - 4 * pi() * a^3
//              @x[3]      - coordinates of a point
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_cdIij( std::complex<double> dJij[27], double mult, double x[3] )
{
  std::complex<double> cx[3] = { x[0], x[1], x[2] };

  this->update_cx( cx, _x1_ );
  dJij[0] = this->cIij_i( mult, cx, _x1_ );
  this->update_cx( cx, _x2_ );
  dJij[1] = this->cIij_i( mult, cx, _x1_ );
  dJij[9] = this->cIij_i( mult, cx, _x2_ );
  this->update_cx( cx, _x3_ );
  dJij[2] = this->cIij_i( mult, cx, _x1_ ); //I11,1; I12,1; I13,1
  dJij[10] = this->cIij_i( mult, cx, _x2_ );
  dJij[18] = this->cIij_i( mult, cx, _x3_ );

  //dIij,1
  dJij[3] = dJij[4] = dJij[5] = dJij[0]; //I21,1; I22,1; I23,1
  dJij[6] = dJij[7] = dJij[8] = dJij[0]; //I31,1; I32,1; I33,1
  //dIij,2
  dJij[11] = dJij[9]; //I11,2; I12,2; I13,2
  dJij[12] = dJij[13] = dJij[14] = dJij[9]; //I21,2; I22,2; I23,2
  dJij[15] = dJij[16] = dJij[17] = dJij[9]; //I31,2; I32,2; I33,2
  //dIij,3
  dJij[19] = dJij[20] = dJij[18]; //I11,3; I12,3; I13,3
  dJij[21] = dJij[22] = dJij[23] = dJij[18]; //I21,3; I22,3; I23,3
  dJij[24] = dJij[25] = dJij[26] = dJij[18]; //I31,3; I32,3; I33,3

  return;
}//end of function: give_cdIij( )

//********************************************************************************************************
// description: Function updates vectors cx according the derivative direction
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @cx[3]     - complex coordinates of a point 'x' (initialized to its real counterparts)
//              @direction - derivative direction (_x1_, _x2_, _x3_)
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::update_cx( std::complex<double> cx[3], derivativeDirection direction )
{
  if( direction == _x1_ ){
    cx[0].imag(_DIFF_H_ );
    cx[1].imag(0);
    cx[2].imag(0);
  }
  else if( direction == _x2_ ){
    cx[0].imag(0);
    cx[1].imag( _DIFF_H_);
    cx[2].imag(0);
  }
  else if( direction == _x3_ ){
    cx[0].imag(0);
    cx[1].imag(0);
    cx[2].imag(_DIFF_H_);
  }
  else{
    _errorr( "update_cx: Invalid derivative direction");
  }

  return ;
}//end of function: update_cx( )

//********************************************************************************************************
// description: Function gives all second derivatives of Ferers-Dysons' integrals Ii by using
//              the complex differentiation
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @ddJi[27]  - table of second derivatives of elliptic integrals to be calculated
//              @dJi[9]    - table of the first complex derivatives of elliptic integrals
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_cddIi( double ddJi[27], std::complex<double> dJi[13] )
{
  //dI1,ij
  ddJi[0] = std::imag( dJi[0] ) / _DIFF_H_; //I1,11
  ddJi[1] = std::imag( dJi[3] ) / _DIFF_H_; //I1,12
  ddJi[2] = std::imag( dJi[6] ) / _DIFF_H_; //I1,13

  ddJi[3] = ddJi[1]; //I1,21
  ddJi[4] = std::imag( dJi[4] ) / _DIFF_H_; //I1,22
  ddJi[5] = std::imag( dJi[7] ) / _DIFF_H_; //I1,23

  ddJi[6] = ddJi[2]; //I1,31
  ddJi[7] = ddJi[5]; //I1,32
  ddJi[8] = std::imag( dJi[8] ) / _DIFF_H_; //I1,33

  //dI1,ij = dI2,ij = dI3,ij
  ddJi[9]  = ddJi[18] = ddJi[0]; //I2,11; I3,11
  ddJi[10] = ddJi[19] = ddJi[1]; //I2,12; I3,12
  ddJi[11] = ddJi[20] = ddJi[2]; //I2,13; I3,13

  ddJi[12] = ddJi[21] = ddJi[3]; //I2,21; I3,21
  ddJi[13] = ddJi[22] = ddJi[4]; //I2,22; I3,22
  ddJi[14] = ddJi[23] = ddJi[5]; //I2,23; I3,23

  ddJi[15] = ddJi[24] = ddJi[6]; //I2,31; I3,31
  ddJi[16] = ddJi[25] = ddJi[7]; //I2,32; I3,32
  ddJi[17] = ddJi[26] = ddJi[8]; //I2,33; I3,33

  return;
}//end of function: give_cddIi( )

//********************************************************************************************************
// description: Function gives all second derivatives of Ferers-Dysons' integrals Iij by using
//              the complex differentiation
// last edit:   08. 07. 2010
//-------------------------------------------------------------------------------------------------------
// note:        @ddJij[81]  - table of second derivatives of elliptic integrals to be calculated
//              @dJij[27]   - table of the first complex derivatives of elliptic integrals
//********************************************************************************************************
void eshelbySoluEllipticIntegralsSphere::give_cddIij( double ddJij[81], std::complex<double> dJij[27] )
{
  int j = 0;

  //dI11,ij
  ddJij[0] = std::imag( dJij[0] ) / _DIFF_H_; //I11,11
  ddJij[1] = std::imag( dJij[1] ) / _DIFF_H_; //I11,12
  ddJij[2] = std::imag( dJij[2] ) / _DIFF_H_; //I11,13

  ddJij[3] = ddJij[1]; //I11,21
  ddJij[4] = std::imag( dJij[9] ) / _DIFF_H_; //I11,22
  ddJij[5] = std::imag( dJij[10] ) / _DIFF_H_; //I11,23

  ddJij[6] = ddJij[2]; //I11,31
  ddJij[7] = ddJij[5]; //I11,32
  ddJij[8] = std::imag( dJij[18] ) / _DIFF_H_;; //I11,33

  //dI11,ij = dI12,ij = dI13,ij = dI21,ij = dI22,ij = dI23,ij = dI31,ij = dI32,ij = dI33,ij
  for( int i = 1; i < 9; i++ ){
    j = i * 9;
    ddJij[j] = ddJij[0];
    ddJij[j + 1] = ddJij[1];
    ddJij[j + 2] = ddJij[2];

    ddJij[j + 3] = ddJij[3];
    ddJij[j + 4] = ddJij[4];
    ddJij[j + 5] = ddJij[5];

    ddJij[j + 6] = ddJij[6];
    ddJij[j + 7] = ddJij[7];
    ddJij[j + 8] = ddJij[8];
  }//end of loop 'for: i'

  return;
}//end of function: give_cddIij( )
// Lambda computed according to Mura p. 86.
double eshelbySoluEllipticIntegralsSphere::getLambda(const double a[3], double x1, double x2, double x3)
{
  // The cubic equation for lambda simplifies to
  return SQR(x1) + SQR(x2) + SQR(x3) - SQR(a[0]);
}


} // end of namespace mumech

/*end of file*/