muMECH  1.0
esei_Ellipsoid.cpp
Go to the documentation of this file.
1 //********************************************************************************************************
2 // code: ### ### ##### #### ## ##
3 // ## ## ######## ## ## ## ## ##
4 // ## ## ## ## ## ### ## ######
5 // ## ## ## ## ## ## ## ## ##
6 // ##### ## ## ##### #### ## ##
7 // ##
8 //
9 // name: eshelbySoluEllipticalIntegralsEllipsoid.cpp
10 // author(s): Jan Novak
11 // language: C, C++
12 // license: This program is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Lesser General Public License as published by
14 // the Free Software Foundation; either version 2 of the License, or
15 // (at your option) any later version.
16 //
17 // This program is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 // GNU Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public License
23 // along with this program; if not, write to the Free Software
24 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
25 //********************************************************************************************************
26 
27 #include "esei_Ellipsoid.h"
28 
29 #include "macros.h"
30 #include "eshelbySoluLambda.h"
31 #include "legendreIntegrals.h"
32 #include "elementaryFunctions.h"
33 #include "problem.h"
34 
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <math.h>
38 #include <cmath>
39 #include <iostream>
40 
41 
42 namespace mumech {
43 
44 //local constant definitions
45 #define _delta_lambda( lambda ) sqrt( ( SQR(I->a[0])+lambda )*( SQR(I->a[1])+lambda )*( SQR(I->a[2])+lambda ) )
46 
47 
48 // ********************************************************************************************************
49 // PUBLIC FUNCTIONS
50 // ********************************************************************************************************
51 
59 void eshelbySoluEllipticIntegralsEllipsoid :: giveEllipticIntegrals (double J[13], double lambda, bool intpoint)
60 {
61  const double *sort_a = I->a;
62 
63  double Theta = asin( sqrt( ( SQR( sort_a[0] ) - SQR( sort_a[2] ) ) / ( SQR( sort_a[0] ) + lambda ) ) );
64  double k = sqrt( ( SQR( sort_a[0] ) - SQR( sort_a[1] ) )/( SQR( sort_a[0] ) - SQR( sort_a[2] ) ) );
65  double mult = _4PI_ * sort_a[0] * sort_a[1] * sort_a[2];
66  double Dla = _delta_lambda( lambda );
67 
68 
69  //J[0] = this->eI( sort_a, k, Theta, mult );//I - this value is not necessary anywhere else
70  J[0] = 0.;
71 
72  J[1] = this->I1( sort_a, k, Theta, mult );//I1
73  J[3] = this->I3( sort_a, k, Theta, lambda, mult, intpoint );//I3
74  J[2] = this->I2( J[1], J[3], lambda, mult, Dla, intpoint );//I2
75 
76  J[5] = this->I12( sort_a, J[1], J[2] );//I12
77  J[6] = this->I13( sort_a, J[1], J[3] );//I13
78  J[4] = this->I11( sort_a, J[5], J[6], lambda, mult, Dla, intpoint );//I11
79 
80  J[7] = J[5];//I21 = I12
81  J[9] = this->I23( sort_a, J[2], J[3] );//I23
82  J[8] = this->I22( sort_a, J[7], J[9], lambda, mult, Dla, intpoint );//I22
83 
84  J[10] = J[6];//I31 = I13
85  J[11] = J[9];//I32 = I23
86  J[12] = this->I33( sort_a, J[11], J[10], lambda, mult, Dla, intpoint );//I33
87 
88 #ifdef TOM_DEBUG
89  std::cout << "J integrals" << std::endl;
90  for (int i = 0; i < 13; i++)
91  std::cout << "J[" << i << "] = " << J[i] << std::endl;
92  std::cout << std::endl;
93 #endif
94 }//end of function: giveEllipticIntegrals( )
95 
96 //********************************************************************************************************
97 // description: Function gives the derivatives of Ferers-Dysons integrals
98 // last edit: 11. 08. 2010
99 //-------------------------------------------------------------------------------------------------------
100 // note: @point - point record (data structure of a given point)
101 // @sort_a[3] - sorted dimensions of ellipsoidal semiaxes (a1>a2,>a3)
102 //********************************************************************************************************
103 void eshelbySoluEllipticIntegralsEllipsoid :: giveDerivativesOfEllipticIntegrals (Point *point, bool intpoint)
104 {
105  if (intpoint) _errorr("derivatives are 0, this function should not be called");
106 
107  double ndiff = I->ndiff_1;
108 
109  switch (I->P->give_diffType()){
110  case DT_COMPLEX :
111  {
112  eshelbySoluLambda Lambda;
113  cLambda CLa;
114 
115  CLa.cx1[0] = (point->loc_x[0],_DIFF_H_); CLa.cx1[1] = point->loc_x[1]; CLa.cx1[2] = point->loc_x[2];
116  CLa.cx2[0] = point->loc_x[0]; CLa.cx2[1] = (point->loc_x[1],_DIFF_H_); CLa.cx2[2] = point->loc_x[2];
117  CLa.cx3[0] = point->loc_x[0]; CLa.cx3[1] = point->loc_x[1]; CLa.cx3[2] = (point->loc_x[2],_DIFF_H_);
118 
119  CLa.cla1 = Lambda.giveLambda( I->a, CLa.cx1, IS_ELLIPSOID );
120  CLa.cla2 = Lambda.giveLambda( I->a, CLa.cx2, IS_ELLIPSOID );
121  CLa.cla3 = Lambda.giveLambda( I->a, CLa.cx3, IS_ELLIPSOID );
122 
123  CLa.cDla1 = _delta_lambda( CLa.cla1 );
124  CLa.cDla2 = _delta_lambda( CLa.cla2 );
125  CLa.cDla3 = _delta_lambda( CLa.cla3 );
126 
127  CLa.cdla_11 = Lambda.c_dLambda( I->a, CLa.cx1, CLa.cla1, _x1_ );
128  CLa.cdla_12 = Lambda.c_dLambda( I->a, CLa.cx2, CLa.cla2, _x1_ );
129  CLa.cdla_13 = Lambda.c_dLambda( I->a, CLa.cx3, CLa.cla3, _x1_ );
130 
131  CLa.cdla_22 = Lambda.c_dLambda( I->a, CLa.cx2, CLa.cla2, _x2_ );
132  CLa.cdla_23 = Lambda.c_dLambda( I->a, CLa.cx3, CLa.cla3, _x2_ );
133 
134  CLa.cdla_33 = Lambda.c_dLambda( I->a, CLa.cx3, CLa.cla3, _x3_ );
135 
136  CLa.mult = -2. * PI * I->a[0] * I->a[1] * I->a[2];
137 
138  //first derivatives Ii,i:
139  this->give_cdIi( point->dJi, point->ddJi, I->a, CLa );
140  //first derivatives Iij,i:
141  this->give_cdIij( point->dJij, point->ddJij, I->a, CLa );
142  //second derivatives Ii,ij:
143  this->give_cddIi( point->ddJi, point->ddJi );
144  //second derivatives Iij,ij:
145  this->give_cddIij( point->ddJij, point->ddJij );
146  return;
147  break;
148  }
149  case DT_NUMERICAL : // Real numerical derivatives
150  {
151  // NUMERICAL DERIVATIVES (using REAL ARGUMENTS)
152 
153  // Get perturbated points and their lambdas
154  double lambdas[NUM_PERTURB];
155  this->getPerturbatedLambdas(lambdas, point->loc_x);
156 
157  // Create and fill a two dimensional array containing values of J-integrals at perturbated points.
158  const int NUM_ELLIP_INT = 13;
159  double J_pert[NUM_PERTURB][NUM_ELLIP_INT]; // All perturbations of all elliptical integrals.
160  for (int i = 0; i < NUM_PERTURB; i++) // For each perturbated point
161  this->giveEllipticIntegrals(J_pert[i], lambdas[i], intpoint); // Get elliptical integrals for lambda of perturbated point.
162 
163  // Compute numerial derivatives of first and second order of J_i integrals.
164  for (int i = 0; i < 3; i++) // For each elliptic integral Ji
165  {
166  // See types.c, section "Derivatives of Ferers-Dysons' integrals".
167  // Ji_j = dJi[3*i+j-4], i,j=1..3.
168  point->dJi[3*i ] = (J_pert[_Xp_][i+1] - J_pert[_Xm_][i+1]) / (2.0 * ndiff);
169  point->dJi[3*i+1] = (J_pert[_Yp_][i+1] - J_pert[_Ym_][i+1]) / (2.0 * ndiff);
170  point->dJi[3*i+2] = (J_pert[_Zp_][i+1] - J_pert[_Zm_][i+1]) / (2.0 * ndiff);
171 
172  // See types.c, section "Derivatives of Ferers-Dysons' integrals".
173  // Ji_jk = ddJi[9*i+3*j+k-13], i,j,k=1..3.
174  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
175  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)
176  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);
177 
178  point->ddJi[9*i+3] = point->ddJi[9*i+1];
179  point->ddJi[9*i+4] = (J_pert[_Yp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Ym_][i+1]) / (ndiff * ndiff);
180  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);
181 
182  point->ddJi[9*i+6] = point->ddJi[9*i+2];
183  point->ddJi[9*i+7] = point->ddJi[9*i+5];
184  point->ddJi[9*i+8] = (J_pert[_Zp_][i+1] - 2.0 * J_pert[_CENTER_][i+1] + J_pert[_Zm_][i+1]) / (ndiff * ndiff);
185  }
186 
187  // Compute numerial derivatives of first and second order of J_ij integrals.
188  for (int i = 0; i < 9; i++) // For each elliptic integral Jij
189  {
190  // See types.c, section "Derivatives of Ferers-Dysons' integrals".
191  // Jij_k = dJij[9*k+3*i+j-13], i,j,k=1..3.
192  point->dJij[i] = (J_pert[_Xp_][i+4] - J_pert[_Xm_][i+4]) / (2.0 * ndiff);
193  point->dJij[i+9] = (J_pert[_Yp_][i+4] - J_pert[_Ym_][i+4]) / (2.0 * ndiff);
194  point->dJij[i+18] = (J_pert[_Zp_][i+4] - J_pert[_Zm_][i+4]) / (2.0 * ndiff);
195 
196  // See types.c, section "Derivatives of Ferers-Dysons' integrals".
197  // Jij_kl = ddJi[27*i+9*j+3*k+l-40] i,j,k,l=1..3.
198  point->ddJij[9*i] = (J_pert[_Xp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Xm_][i+4]) / (ndiff * ndiff);
199  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);
200  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);
201 
202  point->ddJij[9*i+3] = point->ddJij[9*i+1];
203  point->ddJij[9*i+4] = (J_pert[_Yp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Ym_][i+4]) / (ndiff * ndiff);
204  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);
205 
206  point->ddJij[9*i+6] = point->ddJij[9*i+2];
207  point->ddJij[9*i+7] = point->ddJij[9*i+5];
208  point->ddJij[9*i+8] = (J_pert[_Zp_][i+4] - 2.0 * J_pert[_CENTER_][i+4] + J_pert[_Zm_][i+4]) / (ndiff * ndiff);
209  }
210  break;
211  }
212  case DT_ANALITICAL :// Analytical derivatives
213  {
214  eshelbySoluLambda Lambda;
215  aLambda aLa;
216 
217  aLa.la = Lambda.giveLambda( I->a, point->loc_x, IS_ELLIPSOID );
218 
219  aLa.Dla = _delta_lambda( aLa.la );
220 
221  aLa.dla[_x1_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x1_ );
222  aLa.dla[_x2_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x2_ );
223  aLa.dla[_x3_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x3_ );
224 
225  aLa.ddla[_x1_*_x1_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x1_, _x1_ );
226  aLa.ddla[_x1_*_x2_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x1_, _x2_ );
227  aLa.ddla[_x1_*_x3_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x1_, _x3_ );
228 
229  aLa.ddla[_x2_*_x2_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x2_, _x2_ );
230  aLa.ddla[_x2_*_x3_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x2_, _x3_ );
231 
232  aLa.ddla[_x3_*_x3_] = Lambda.giveLambdaDerivative( I->a, point->loc_x, aLa.la, _x3_, _x3_ );
233 
234  aLa.mult = -2. * PI * I->a[0] * I->a[1] * I->a[2];
235 
236  //first derivatives Ii,i:
237  this->give_dIi( point->dJi, I->a, point->loc_x, aLa );
238  //first derivatives Iij,i:
239  this->give_dIij( point->dJij, I->a, point->loc_x, aLa );
240  //second derivatives Ii,ij:
241  this->give_ddIi( point->ddJi, I->a, point->loc_x, aLa );
242  //second derivatives Iij,ij:
243  this->give_ddIij( point->ddJij, I->a, point->loc_x, aLa );
244  break;
245  }
246  default:
247  _errorr( "Invalid type of differentiation.\n");
248  }
249 
250 #ifdef TOM_DEBUG
251  std::cout << "dJi integrals" << std::endl;
252  for (int i = 0; i < 9; i++)
253  std::cout << "dJi[" << i << "] = " << point->dJi[i] << std::endl;
254  std::cout << std::endl;
255  std::cout << "dJij integrals" << std::endl;
256  for (int i = 0; i < 27; i++)
257  std::cout << "dJij[" << i << "] = " << point->dJij[i] << std::endl;
258  std::cout << std::endl;
259  std::cout << "ddJi integrals" << std::endl;
260  for (int i = 0; i < 27; i++)
261  std::cout << "ddJi[" << i << "] = " << point->ddJi[i] << std::endl;
262  std::cout << std::endl;
263  std::cout << "ddJij integrals" << std::endl;
264  for (int i = 0; i < 81; i++)
265  std::cout << "ddJij[" << i << "] = " << point->ddJij[i] << std::endl;
266  std::cout << std::endl;
267 #endif
268 
269 }//end of function: giveDerivativeOfEllipticIntegrals( )
270 
271 
272 //********************************************************************************************************
273 // PROTECTED FUNCTIONS
274 //********************************************************************************************************
275 
276 //Ferers-Dysons integrals:
277 //------------------------
278 
279 //********************************************************************************************************
280 // description: Function gives the solution of the integral I
281 // last edit: 25. 08. 2010
282 //-------------------------------------------------------------------------------------------------------
283 // note: Solution adopted from Toshio Mura's book (1982)
284 // @a[3] - semiaxes, sorted as a1>a2>a3
285 // @k, Theta - limits of Legendre's elliptic integrals
286 // @mult = 4 * PI * a1 * a2 * a3
287 //********************************************************************************************************
288 double eshelbySoluEllipticIntegralsEllipsoid::eI( double a[3], double k, double Theta, double mult )
289 {
290  return ( mult * legendreIntegrals::ellf (Theta, k) / sqrt( SQR( a[0] ) - SQR( a[2] ) ) );
291 }
292 
293 //********************************************************************************************************
294 // description: Function gives the solution of the integral I_1
295 // last edit: 25. 08. 2010
296 //--------------------------------------------------------------------------------------------------------
297 // note: Solution adopted from Toshio Mura's book (1982).
298 // This solution is apparently wrong in Eq. 12.17, page 78.
299 // It holds for interior points only
300 //--------------------------------------------------------------------------------------------------------
301 // @a[3] - semiaxes, sorted as a1>a2>a3
302 // @k, Theta - limits of Legendre's elliptic integrals
303 // @mult = 4 * PI * a1 * a2 * a3
304 //********************************************************************************************************
305 double eshelbySoluEllipticIntegralsEllipsoid::I1(const double a[3], double k, double Theta, double mult )
306 {
307  // Mura book (12.17)
308  double a1 = a[0], a2 = a[1], a3 = a[2];
309  legendreIntegrals ellInt;
310 
311  // Mura book (11.17) or (12.17)
312  return( mult * ( ellInt.ellf( Theta, k ) - ellInt.elle( Theta, k ) )/
313  ( ( SQR( a1 ) - SQR( a2 ) ) * sqrt( SQR( a1 ) - SQR( a3 ) ) ) );
314 }//end of function: I1()
315 
316 //********************************************************************************************************
317 // description: Function gives the solution of the integral I_2
318 // last edit: 25. 08. 2010
319 //-------------------------------------------------------------------------------------------------------
320 // note: Solution adopted from Toshio Mura's book (1982)
321 // @a[3] - semiaxes, sorted as a1>a2>a3
322 // @lambda - lambda parameter according Mura's book, page 73
323 // @mult = 4 * PI * a1 * a2 * a3
324 // @Dla - Delta lambda, see the definition above
325 // intpoint - internal x external point
326 //********************************************************************************************************
327 double eshelbySoluEllipticIntegralsEllipsoid::I2( const double I1, double I3, double lambda,
328  double mult, double Dla, bool intpoint)
329 {
330  if (intpoint) return ( _4PI_ - I1 - I3 ); // internal points
331  else return ( mult / Dla - I1 - I3 ); // external points
332 }
333 
334 //********************************************************************************************************
335 // description: Function gives the solution of the integral I_3
336 // last edit: 25. 08. 2010
337 //-------------------------------------------------------------------------------------------------------
338 // note: Solution adopted from Toshio Mura's book (1982)
339 // @a[3] - semiaxes, sorted as a1>a2>a3
340 // @k, Theta - limits of Legendre's elliptic integrals
341 // @lambda - lambda parameter according Mura's book, page 73
342 // @mult = 4 * PI * a1 * a2 * a3
343 //********************************************************************************************************
344 double eshelbySoluEllipticIntegralsEllipsoid::I3( const double a[3], double k, double Theta, double lambda,
345  double mult, bool intpoint)
346 {
347  double a1 = a[0], a2 = a[1], a3 = a[2], value = 0.;
348  legendreIntegrals ellInt;
349 
350  if (intpoint) {//internal points
351  value = mult / ( ( SQR( a2 ) - SQR( a3 ) ) * sqrt( SQR( a1 ) - SQR( a3 ) ) ) *
352  ( a2 * sqrt( SQR(a1) - SQR(a3) )/( a1 * a3 ) - ellInt.elle( Theta, k ) );
353  }
354  else {//external points
355  value = mult / ( ( SQR( a2 ) - SQR( a3 ) ) * sqrt( SQR( a1 ) - SQR( a3 ) ) ) *
356  ( sqrt( SQR( a2 ) + lambda ) * sqrt( SQR( a1 ) - SQR( a3 ) ) * pow( ( SQR( a3 ) + lambda ), -.5 ) *
357  pow( ( SQR( a1 ) + lambda ), -.5 ) - ellInt.elle( Theta, k ) );
358  }
359 
360  return value;
361 }//end of function: I3()
362 
363 //********************************************************************************************************
364 // description: Function gives the solution of the integral I_11
365 // last edit: 25. 08. 2010
366 //-------------------------------------------------------------------------------------------------------
367 // note: Solution adopted from Toshio Mura's book (1982)
368 // @mult = 4 * PI * a1 * a2 * a3
369 // @Dla - Delta lambda, see the definition above
370 //********************************************************************************************************
371 double eshelbySoluEllipticIntegralsEllipsoid::I11(const double a[3], double I12, double I13, double lambda,
372  double mult, double Dla, bool intpoint)
373 {
374  if (intpoint) return ( _4PI_ / SQR( a[0] ) - I12 - I13 ) / 3.; // internal points
375  else return ( mult / ( ( SQR( a[0] ) + lambda ) * Dla ) - I12 - I13 ) / 3.; // external points
376 }
377 
378 //********************************************************************************************************
379 // description: Function gives the solution of the integral I_12
380 // last edit: 18. 08. 2009
381 //-------------------------------------------------------------------------------------------------------
382 // note: Solution adopted from Toshio Mura's book (1982)
383 //********************************************************************************************************
384 double eshelbySoluEllipticIntegralsEllipsoid::I12(const double a[3], double I1, double I2 )
385 {
386  return ( ( I2 - I1 )/( SQR( a[0] ) - SQR( a[1] ) ) );
387 }
388 
389 //********************************************************************************************************
390 // description: Function gives the solution of the integral I_13
391 // last edit: 18. 08. 2009
392 //-------------------------------------------------------------------------------------------------------
393 // note: Solution adopted from Toshio Mura's book (1982)
394 //********************************************************************************************************
395 double eshelbySoluEllipticIntegralsEllipsoid::I13(const double a[3], double I1, double I3 )
396 {
397  return( ( I3 - I1 )/( SQR( a[0] ) - SQR( a[2] ) ) );
398 }
399 
400 //********************************************************************************************************
401 // description: Function gives the solution of the integral I_22
402 // last edit: 25. 08. 2010
403 //-------------------------------------------------------------------------------------------------------
404 // note: Solution adopted from Toshio Mura's book (1982)
405 // @mult = 4 * PI * a1 * a2 * a3
406 // @Dla - Delta lambda, see the definition above
407 //********************************************************************************************************
408 double eshelbySoluEllipticIntegralsEllipsoid::I22(const double a[3], double I21, double I23, double lambda,
409  double mult, double Dla, bool intpoint )
410 {
411  if (intpoint) return ( _4PI_ / SQR( a[1] ) - I21 - I23 ) / 3.; // internal points
412  else return ( mult / ( ( SQR( a[1] ) + lambda ) * Dla ) - I21 - I23 ) / 3.; // external points
413 }
414 
415 //********************************************************************************************************
416 // description: Function gives the solution of the integral I_23
417 // last edit: 19. 08. 2009
418 //-------------------------------------------------------------------------------------------------------
419 // note: Solution adopted from Toshio Mura's book (1982)
420 //********************************************************************************************************
421 double eshelbySoluEllipticIntegralsEllipsoid::I23(const double a[3], double I2, double I3 )
422 {
423  return( ( I3 - I2 )/( SQR( a[1] ) - SQR( a[2] ) ) );
424 }
425 
426 //********************************************************************************************************
427 // description: Function gives the solution of the integral I_33
428 // last edit: 25. 08. 2010
429 //-------------------------------------------------------------------------------------------------------
430 // note: Solution adopted from Toshio Mura's book (1982)
431 // @mult = 4 * PI * a1 * a2 * a3
432 // @Dla - Delta lambda, see the definition above
433 //********************************************************************************************************
434 double eshelbySoluEllipticIntegralsEllipsoid::I33(const double a[3], double I32, double I31, double lambda,
435  double mult, double Dla, bool intpoint )
436 {
437  if (intpoint) return ( _4PI_ / SQR( a[2] ) - I31 - I32 ) / 3.; // internal points
438  else return ( mult / ( ( SQR( a[2] ) + lambda ) * Dla ) - I31 - I32 ) / 3.; // external points
439 }
440 
441 
442 //Analytical derivatives of Ferers-Dysons integrals:
443 //--------------------------------------------------
444 
445 //********************************************************************************************************
446 // description: Function gives the first or second derivative of a component of elliptical integrals
447 // last edit: 25. 08. 2010
448 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
449 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
450 // @x[3] - coordinates of a point
451 // @aLa - lambda parameter and other related quantities
452 // @component - component which derivative should be calculated
453 // @direction_1 - direction of the first derivative
454 // @direction_2 - direction of the seccond derivative: _x1_, _x2_, _x3_ or must be '_empty_'
455 // to return the first derivative according 'direction_1' only
456 //********************************************************************************************************
457 double eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegralDerivative(const double a[3], const double x[3],
458  aLambda aLa,
459  ellipticIntegralComponent component,
460  derivativeDirection direction_1,
461  derivativeDirection direction_2 )
462 {
463  double derivative = 0.;
464 
465  if( ( direction_1 == _x1_ || direction_1 == _x2_ || direction_1 == _x3_ ) && direction_2 == _empty_ ){
466  derivative = this->giveEllipticIntegralFirstDerivative( a, x, aLa, component, direction_1 );
467  }
468  else if( ( direction_1 == _x1_ || direction_1 == _x2_ || direction_1 == _x3_ ) &&
469  ( direction_2 == _x1_ || direction_2 == _x2_ || direction_2 == _x3_ ) ){
470  derivative = this->giveEllipticIntegralSeccondDerivative( a, x, aLa, component, direction_1, direction_2 );
471  }
472  else{
473  _errorr( "giveEllipticalIntegralDerivative: invalid derivative direction ");
474  }
475 
476  return derivative;
477 }//end of function: giveEllipticIntegralDerivative( )
478 
479 //********************************************************************************************************
480 // description: Function gives the first derivative of the Elliptic integral
481 // last edit: 25. 08. 2010
482 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
483 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
484 // @x[3] - coordinates of a point
485 // @aLa - lambda parameter and other related quantities
486 // @component - component which derivative should be calculated
487 // @direction - derivative direction: _x1_, _x2_, _x3_
488 //********************************************************************************************************
489 double eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegralFirstDerivative(const double a[3], const double x[3],
490  aLambda aLa,
491  ellipticIntegralComponent component,
492  derivativeDirection direction )
493 {
494  double derivative = 0.;
495 
496  if( component == _I_ ){
497  derivative = this->give_Ip_Derivative(aLa, direction );
498  }
499  else if( component == _I1_ || component == _I2_ || component == _I3_ ){
500  derivative = this->give_Iip_Derivative( a, aLa, component, direction );
501  }else if( component == _I11_ || component == _I12_ || component == _I13_ ||
502  component == _I21_ || component == _I22_ || component == _I23_ ||
503  component == _I31_ || component == _I32_ || component == _I33_ ){
504  derivative = this->give_Iijp_Derivative( a, aLa, component, direction );
505  }
506  else{
507  _errorr( "giveEllipticIntegralFirstDerivative: invalid integral component ");
508  }
509 
510  return derivative;
511 }//end of function: giveEllipticIntegralFirstDerivative()
512 
513 //********************************************************************************************************
514 // description: Function gives the second derivative of the Elliptic integral
515 // last edit: 25. 08. 2010
516 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
517 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
518 // @x[3] - coordinates of a point
519 // @aLa - lambda parameter and other related quantities
520 // @component - component which derivative should be calculated
521 // @direction_1 - direction of the first derivative: _x1_, _x2_, _x3_
522 // @direction_2 - direction of the second derivative: _x1_, _x2_, _x3_ or must be '_empty_'
523 // to return the first derivative according 'direction_1' only
524 //********************************************************************************************************
525 double eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegralSeccondDerivative(const double a[3], const double x[3],
526  aLambda aLa,
527  ellipticIntegralComponent component,
528  derivativeDirection direction_1,
529  derivativeDirection direction_2 )
530 {
531  double derivative = 0.;
532 
533  if( component == _I1_ || component == _I2_ || component == _I3_ ){
534  derivative = this->give_Iipq_Derivative( a, aLa, component, direction_1, direction_2 );
535  }else if( component == _I11_ || component == _I12_ || component == _I13_ ||
536  component == _I21_ || component == _I22_ || component == _I23_ ||
537  component == _I31_ || component == _I32_ || component == _I33_ ){
538  derivative = this->give_Iijpq_Derivative( a, aLa, component, direction_1, direction_2 );
539  }
540  else{
541  _errorr( "giveEllipticIntegralSeccondDerivative: invalid integral component ");
542  }
543 
544  return derivative;
545 }//end of function: giveEllipticIntegralSeccondDerivative()
546 
547 //********************************************************************************************************
548 // description: Function gives the first derivative I,p of Elliptic integral I according 'p' direction
549 // last edit: 25. 08. 2010
550 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
551 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
552 // @x[3] - coordinates of a point
553 // @aLa - lambda parameter and other related quantities
554 // @p - derivative direction: _x1_, _x2_, _x3_
555 //********************************************************************************************************
556 double eshelbySoluEllipticIntegralsEllipsoid::give_Ip_Derivative(aLambda aLa, derivativeDirection p)
557 {
558  //returned (integral) value calculation
559  return ( aLa.mult * aLa.dla[p] / aLa.Dla );
560 }
561 
562 //********************************************************************************************************
563 // description: Function gives the first derivative I_i,p of Elliptic integral I_i according 'p' direction
564 // last edit: 25. 08. 2010
565 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
566 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
567 // @x[3] - coordinates of a point
568 // @aLa - lambda parameter and other related quantities
569 // @component - component which derivative should be calculated
570 // @p - derivative direction: _x1_, _x2_, _x3_
571 //********************************************************************************************************
572 double eshelbySoluEllipticIntegralsEllipsoid::give_Iip_Derivative(const double a[3], aLambda aLa,
573  ellipticIntegralComponent component,
574  derivativeDirection p )
575 {
576  int i;
577 
578  switch( ( int ) component ){
579  case _I1_ : i = 0; break;
580  case _I2_ : i = 1; break;
581  case _I3_ : i = 2; break;
582  default:
583  _errorr( "give_Iip_Derivative: invalid integral component ");
584  }
585 
586  //returned (integral) value
587  return( aLa.mult * aLa.dla[p] / ( ( SQR( a[i] ) + aLa.la ) * aLa.Dla ) );
588 }
589 
590 //********************************************************************************************************
591 // description: Function gives the first derivative I_ij,p of Elliptic integral I_ij according 'p'
592 // direction
593 // last edit: 25. 08. 2010
594 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
595 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
596 // @x[3] - coordinates of a point
597 // @aLa - lambda parameter and other related quantities
598 // @component - component which derivative should be calculated
599 // @p - derivative direction: _x1_, _x2_, _x3_
600 //********************************************************************************************************
601 double eshelbySoluEllipticIntegralsEllipsoid::give_Iijp_Derivative(const double a[3], aLambda aLa,
602  ellipticIntegralComponent component,
603  derivativeDirection p )
604 {
605  int i, j;
606 
607  switch( ( int ) component ){
608  case _I11_ : i = 0; j = 0; break;
609  case _I12_ : i = 0; j = 1; break;
610  case _I13_ : i = 0; j = 2; break;
611  case _I21_ : i = 1; j = 0; break;
612  case _I22_ : i = 1; j = 1; break;
613  case _I23_ : i = 1; j = 2; break;
614  case _I31_ : i = 2; j = 0; break;
615  case _I32_ : i = 2; j = 1; break;
616  case _I33_ : i = 2; j = 2; break;
617  default:
618  _errorr( "give_Iijp_Derivative: invalid integral component ");
619  }
620 
621  //returned (integral) value
622  return( aLa.mult * aLa.dla[p] / ( ( SQR( a[i] ) + aLa.la ) * ( SQR( a[j] ) + aLa.la ) * aLa.Dla ) );
623 }
624 
625 //********************************************************************************************************
626 // description: Function gives the second derivative of the I_ij integral by 'p' and 'q' directions
627 // last edit: 25. 8. 2010
628 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
629 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
630 // @x[3] - coordinates of a point
631 // @aLa - lambda parameter and other related quantities
632 // @component - component which derivative should be calculated
633 // @p - first derivative direction: _x1_, _x2_, _x3_
634 // @q - second derivative direction: _x1_, _x2_, _x3_
635 //********************************************************************************************************
636 double eshelbySoluEllipticIntegralsEllipsoid::give_Iijpq_Derivative(const double a[3], aLambda aLa,
637  ellipticIntegralComponent component,
638  derivativeDirection p,
639  derivativeDirection q )
640 {
641  double value = 0.;
642  double dLa_p = aLa.dla[p]; //first derivative of lambda by 'p'
643  double dLa_q = aLa.dla[q]; //first derivative of lambda by 'q'
644  double dLa_pq = aLa.ddla[p * q]; //second derivative of lambda by 'pq'
645  double a1 = a[0], a2 = a[1], a3 = a[2];
646  int i, j;
647 
648  switch( component ){
649  case _I11_ : i = 0; j = 0; break;
650  case _I12_ : i = 0; j = 1; break;
651  case _I13_ : i = 0; j = 2; break;
652  case _I21_ : i = 1; j = 0; break;
653  case _I22_ : i = 1; j = 1; break;
654  case _I23_ : i = 1; j = 2; break;
655  case _I31_ : i = 2; j = 0; break;
656  case _I32_ : i = 2; j = 1; break;
657  case _I33_ : i = 2; j = 2; break;
658  default:
659  _errorr( "give_Iijp_Derivative: invalid integral component ");
660  }
661 
662  value = (a1*a2*a3*PI*(2*(SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la)*(SQR(a[i]) + aLa.la)*
663  dLa_q*dLa_p + 2*(SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la)*
664  (SQR(a[j]) + aLa.la)*dLa_q*dLa_p + (SQR(a[i]) + aLa.la)*(SQR(a[j]) + aLa.la)*
665  (SQR(a2)*SQR(a3) + SQR(a1)*(SQR(a2) + SQR(a3)) +
666  2*(SQR(a1) + SQR(a2) + SQR(a3))*aLa.la + 3*SQR(aLa.la))*dLa_q*dLa_p -
667  2*(SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la)*(SQR(a[i]) + aLa.la)*
668  (SQR(a[j]) + aLa.la)*dLa_pq))/
669  (pow((SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la),1.5)*SQR(SQR(a[i]) + aLa.la)*
670  SQR(SQR(a[j]) + aLa.la));
671 
672  return value;
673 }
674 
675 //********************************************************************************************************
676 // description: Function gives the second derivative of the I_i integral by 'p' and 'q' directions
677 // last edit: 25. 8. 2010
678 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
679 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
680 // @x[3] - coordinates of a point
681 // @aLa - lambda parameter and other related quantities
682 // @component - component which derivative should be calculated
683 // @p - first derivative direction: _x1_, _x2_, _x3_
684 // @q - second derivative direction: _x1_, _x2_, _x3_
685 //********************************************************************************************************
686 double eshelbySoluEllipticIntegralsEllipsoid::give_Iipq_Derivative(const double a[3], aLambda aLa,
687  ellipticIntegralComponent component,
688  derivativeDirection p,
689  derivativeDirection q )
690 {
691  double value = 0.;
692  double dLa_p = aLa.dla[p]; //first derivative of lambda by 'p'
693  double dLa_q = aLa.dla[q]; //first derivative of lambda by 'q'
694  double dLa_pq = aLa.ddla[p * q]; //second derivative of lambda by 'pq'
695  double a1 = a[0], a2 = a[1], a3 = a[2];
696  int i;
697 
698  switch( component ){
699  case _I1_ : i = 0; break;
700  case _I2_ : i = 1; break;
701  case _I3_ : i = 2; break;
702  default:
703  _errorr( "give_Iipq_Derivative: invalid integral component ");
704  }
705 
706  value = (a1*a2*a3*PI*(2*(SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la)*dLa_q*
707  dLa_p + (SQR(a[i]) + aLa.la)*
708  (SQR(a2)*SQR(a3) + SQR(a1)*(SQR(a2) + SQR(a3)) +
709  2*(SQR(a1) + SQR(a2) + SQR(a3))*aLa.la + 3*SQR(aLa.la))*dLa_q*dLa_p -
710  2*(SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la)*(SQR(a[i]) + aLa.la)*dLa_pq))/
711  (pow((SQR(a1) + aLa.la)*(SQR(a2) + aLa.la)*(SQR(a3) + aLa.la),1.5)*SQR(SQR(a[i]) + aLa.la));
712 
713  return value;
714 }
715 
716 //********************************************************************************************************
717 // description: Function gives all first derivatives of Ferers-Dysons' integrals Ii
718 // last edit: 25. 08. 2010
719 //-------------------------------------------------------------------------------------------------------
720 // note: @dJi[9] - table of elliptic integral derivatives to be calculated
721 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
722 // @x[3] - coordinates of a point
723 // @aLa - lambda parameter and other related quantities
724 //********************************************************************************************************
725 void eshelbySoluEllipticIntegralsEllipsoid::give_dIi( double dJi[9], const double sort_a[3], const double x[3],
726  aLambda aLa )
727 {
728  //dIi,1
729  dJi[0] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x1_, _empty_ );//I1,1
730  dJi[3] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x1_, _empty_ );//I2,1
731  dJi[6] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x1_, _empty_ );//I3,1
732  //dIi,2
733  dJi[1] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x2_, _empty_ );//I1,2
734  dJi[4] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x2_, _empty_ );//I2,2
735  dJi[7] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x2_, _empty_ );//I3,2
736  //dIi,3
737  dJi[2] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x3_, _empty_ );//I1,3
738  dJi[5] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x3_, _empty_ );//I2,3
739  dJi[8] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x3_, _empty_ );//I3,3
740 
741  return;
742 }
743 
744 //********************************************************************************************************
745 // description: Function gives all first derivatives of Ferers-Dysons' integrals Iij
746 // last edit: 25. 08. 2010
747 //-------------------------------------------------------------------------------------------------------
748 // note: @dJij[27] - table of elliptic integral derivatives to be calculated
749 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
750 // @x[3] - coordinates of a point
751 // @aLa - lambda parameter and other related quantities
752 //********************************************************************************************************
753 void eshelbySoluEllipticIntegralsEllipsoid::give_dIij( double dJij[27], const double sort_a[3], const double x[3],
754  aLambda aLa )
755 {
756  //dIij,1
757  dJij[0] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x1_, _empty_ );//I11,1
758  dJij[1] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x1_, _empty_ );//I12,1
759  dJij[2] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x1_, _empty_ );//I13,1
760 
761  dJij[3] = dJij[1];//I21,1
762  dJij[4] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x1_, _empty_ );//I22,1
763  dJij[5] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x1_, _empty_ );//I23,1
764 
765  dJij[6] = dJij[2];//I31,1
766  dJij[7] = dJij[5];//I32,1
767  dJij[8] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x1_, _empty_ );//I33,1
768 
769  //dIij,2
770  dJij[9] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x2_, _empty_ );//I11,2
771  dJij[10] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x2_, _empty_ );//I12,2
772  dJij[11] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x2_, _empty_ );//I13,2
773 
774  dJij[12] = dJij[10];//I21,2
775  dJij[13] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x2_, _empty_ );//I22,2
776  dJij[14] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x2_, _empty_ );//I23,2
777 
778  dJij[15] = dJij[11];//I31,2
779  dJij[16] = dJij[14];//I32,2
780  dJij[17] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x2_, _empty_ );//I33,2
781 
782  //dIij,3
783  dJij[18] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x3_, _empty_ );//I11,3
784  dJij[19] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x3_, _empty_ );//I12,3
785  dJij[20] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x3_, _empty_ );//I13,3
786 
787  dJij[21] = dJij[19];//I21,3
788  dJij[22] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x3_, _empty_ );//I22,3
789  dJij[23] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x3_, _empty_ );//I23,3
790 
791  dJij[24] = dJij[20];//I31,3
792  dJij[25] = dJij[23];//I32,3
793  dJij[26] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x3_, _empty_ );//I33,3
794 
795  return;
796 }//end of function: give_dIij( )
797 
798 //********************************************************************************************************
799 // description: Function gives all second derivatives of Ferers-Dysons' integrals Ii
800 // last edit: 25. 08. 2010
801 //-------------------------------------------------------------------------------------------------------
802 // note: @ddJi[27] - table of elliptic integral derivatives to be calculated
803 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
804 // @x[3] - coordinates of a point
805 // @aLa - lambda parameter and other related quantities
806 //********************************************************************************************************
807 void eshelbySoluEllipticIntegralsEllipsoid::give_ddIi( double ddJi[27], const double sort_a[3], const double x[3],
808  aLambda aLa )
809 {
810  //dI1,ij
811  ddJi[0] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x1_, _x1_ );//I1,11
812  ddJi[1] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x1_, _x2_ );//I1,12
813  ddJi[2] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x1_, _x3_ );//I1,13
814 
815  ddJi[3] = ddJi[1]; //I1,21
816  ddJi[4] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x2_, _x2_ );//I1,22
817  ddJi[5] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x2_, _x3_ );//I1,23
818 
819  ddJi[6] = ddJi[2]; //I1,31
820  ddJi[7] = ddJi[5]; //I1,32
821  ddJi[8] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I1_, _x3_, _x3_ );//I1,33
822 
823  //dI2,ij
824  ddJi[9] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x1_, _x1_ );//I2,11
825  ddJi[10] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x1_, _x2_ );//I2,12
826  ddJi[11] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x1_, _x3_ );//I2,13
827 
828  ddJi[12] = ddJi[10];//I2,21
829  ddJi[13] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x2_, _x2_ );//I2,22
830  ddJi[14] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x2_, _x3_ );//I2,23
831 
832  ddJi[15] = ddJi[11];//I2,31
833  ddJi[16] = ddJi[14];//I2,32
834  ddJi[17] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I2_, _x3_, _x3_ );//I2,33
835 
836  //dI3,ij
837  ddJi[18] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x1_, _x1_ );//I3,11
838  ddJi[19] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x1_, _x2_ );//I3,12
839  ddJi[20] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x1_, _x3_ );//I3,13
840 
841  ddJi[21] = ddJi[19];//I3,21
842  ddJi[22] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x2_, _x2_ );//I3,22
843  ddJi[23] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x2_, _x3_ );//I3,23
844 
845  ddJi[24] = ddJi[20];//I3,31
846  ddJi[25] = ddJi[23];//I3,32
847  ddJi[26] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I3_, _x3_, _x3_ );//I3,33
848 
849  return;
850 }//end of function: give_ddIi( )
851 
852 //********************************************************************************************************
853 // description: Function gives all second derivatives of Ferers-Dysons' integrals Iij
854 // last edit: 25. 08. 2010
855 //-------------------------------------------------------------------------------------------------------
856 // note: @ddJij[81] - table of elliptic integral derivatives to be calculated
857 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
858 // @x[3] - coordinates of a point
859 // @aLa - lambda parameter and other related quantities
860 //********************************************************************************************************
861 void eshelbySoluEllipticIntegralsEllipsoid::give_ddIij( double ddJij[27], const double sort_a[3], const double x[3],
862  aLambda aLa )
863 {
864  //ddI11,ij
865  ddJij[0] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x1_, _x1_ );//I11,11
866  ddJij[1] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x1_, _x2_ );//I11,12
867  ddJij[2] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x1_, _x3_ );//I11,13
868 
869  ddJij[3] = ddJij[1]; //I11,21
870  ddJij[4] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x2_, _x2_ );//I11,22
871  ddJij[5] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x2_, _x3_ );//I11,23
872 
873  ddJij[6] = ddJij[2]; //I11,31
874  ddJij[7] = ddJij[5]; //I11,32
875  ddJij[8] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I11_, _x3_, _x3_ );//I11,33
876 
877  //ddI12,ij
878  ddJij[9] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x1_, _x1_ );//I12,11
879  ddJij[10] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x1_, _x2_ );//I12,12
880  ddJij[11] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x1_, _x3_ );//I12,13
881 
882  ddJij[12] = ddJij[10]; //I12,21
883  ddJij[13] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x2_, _x2_ );//I12,22
884  ddJij[14] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x2_, _x3_ );//I12,23
885 
886  ddJij[15] = ddJij[11]; //I12,31
887  ddJij[16] = ddJij[14]; //I12,32
888  ddJij[17] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I12_, _x3_, _x3_ );//I12,33
889 
890  //ddI13,ij
891  ddJij[18] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x1_, _x1_ );//I13,11
892  ddJij[19] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x1_, _x2_ );//I13,12
893  ddJij[20] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x1_, _x3_ );//I13,13
894 
895  ddJij[21] = ddJij[19]; //I13,21
896  ddJij[22] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x2_, _x2_ );//I13,22
897  ddJij[23] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x2_, _x3_ );//I13,23
898 
899  ddJij[24] = ddJij[20]; //I13,31
900  ddJij[25] = ddJij[23]; //I13,32
901  ddJij[26] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I13_, _x3_, _x3_ );//I13,33
902 
903  //ddI21,ij
904  ddJij[27] = ddJij[9]; //I21,11
905  ddJij[28] = ddJij[10];//I21,12
906  ddJij[29] = ddJij[11];//I21,13
907 
908  ddJij[30] = ddJij[12];//I21,21
909  ddJij[31] = ddJij[13];//I21,22
910  ddJij[32] = ddJij[14];//I21,23
911 
912  ddJij[33] = ddJij[15];//I21,31
913  ddJij[34] = ddJij[16];//I21,32
914  ddJij[35] = ddJij[17];//I21,33
915 
916  //ddI22,ij
917  ddJij[36] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x1_, _x1_ );//I22,11
918  ddJij[37] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x1_, _x2_ );//I22,12
919  ddJij[38] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x1_, _x3_ );//I22,13
920 
921  ddJij[39] = ddJij[37]; //I22,21
922  ddJij[40] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x2_, _x2_ );//I22,22
923  ddJij[41] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x2_, _x3_ );//I22,23
924 
925  ddJij[42] = ddJij[38]; //I22,31
926  ddJij[43] = ddJij[41]; //I22,32
927  ddJij[44] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I22_, _x3_, _x3_ );//I22,33
928 
929  //ddI23,ij
930  ddJij[45] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x1_, _x1_ );//I23,11
931  ddJij[46] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x1_, _x2_ );//I23,12
932  ddJij[47] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x1_, _x3_ );//I23,13
933 
934  ddJij[48] = ddJij[46]; //I23,21
935  ddJij[49] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x2_, _x2_ );//I23,22
936  ddJij[50] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x2_, _x3_ );//I23,23
937 
938  ddJij[51] = ddJij[47]; //I23,31
939  ddJij[52] = ddJij[50]; //I23,32
940  ddJij[53] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I23_, _x3_, _x3_ );//I23,33
941 
942  //ddI31,ij
943  ddJij[54] = ddJij[18];//I31,11
944  ddJij[55] = ddJij[19];//I31,12
945  ddJij[56] = ddJij[20];//I31,13
946 
947  ddJij[57] = ddJij[21];//I31,21
948  ddJij[58] = ddJij[22];//I31,22
949  ddJij[59] = ddJij[23];//I31,23
950 
951  ddJij[60] = ddJij[24];//I31,31
952  ddJij[61] = ddJij[25];//I31,32
953  ddJij[62] = ddJij[26];//I31,33
954 
955  //ddI32,ij
956  ddJij[63] = ddJij[45];//I32,11
957  ddJij[64] = ddJij[46];//I32,12
958  ddJij[65] = ddJij[47];//I32,13
959 
960  ddJij[66] = ddJij[48];//I32,21
961  ddJij[67] = ddJij[49];//I32,22
962  ddJij[68] = ddJij[50];//I32,23
963 
964  ddJij[69] = ddJij[51];//I32,31
965  ddJij[70] = ddJij[52];//I32,32
966  ddJij[71] = ddJij[53];//I32,33
967 
968  //ddI33,ij
969  ddJij[72] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x1_, _x1_ );//I33,11
970  ddJij[73] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x1_, _x2_ );//I33,12
971  ddJij[74] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x1_, _x3_ );//I33,13
972 
973  ddJij[75] = ddJij[73]; //I33,21
974  ddJij[76] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x2_, _x2_ );//I33,22
975  ddJij[77] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x2_, _x3_ );//I33,23
976 
977  ddJij[78] = ddJij[74]; //I33,31
978  ddJij[79] = ddJij[77]; //I33,32
979  ddJij[80] = this->giveEllipticIntegralDerivative( sort_a, x, aLa, _I33_, _x3_, _x3_ );//I33,33
980 
981  return;
982 }//end of function: give_ddIij( )
983 
984 //Complex numerical derivatives of Ferers-Dysons integrals:
985 //---------------------------------------------------------
986 
987 //********************************************************************************************************
988 // description: Function gives the first complex derivative of Ii_i in 'i' direction
989 // last edit: 02. 08. 2010
990 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
991 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
992 // @la - lambda parameter (not a Lame constant !!!)
993 // @Dla - Delta lambda
994 // @cdla_i - first complex derivative of lambda in 'i' direction
995 // @mult - is the multiplier calculated as: -2. * PI * a[0] * a[1] * a[2]
996 // @component - component which derivative should be calculated
997 //********************************************************************************************************
998 std::complex<double> eshelbySoluEllipticIntegralsEllipsoid::cIi_i(const double a[3], std::complex<double> la,
999  std::complex<double> Dla, std::complex<double> cdla_i,
1000  double mult,
1001  ellipticIntegralComponent component )
1002 {
1003  int comp = component - 1;
1004 
1005  return( mult * cdla_i / ( ( SQR( a[comp] ) + la ) * Dla ) );
1006 }
1007 
1008 //********************************************************************************************************
1009 // description: Function gives the first complex derivative of Iij_i in 'i' direction
1010 // last edit: 04. 08. 2010
1011 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
1012 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
1013 // @la - lambda parameter (not a Lame constant !!!)
1014 // @Dla - Delta lambda
1015 // @cdla_i - first complex derivative of lambda in 'i' direction
1016 // @mult - is the multiplier calculated as: -2. * PI * a[0] * a[1] * a[2]
1017 // @component - component which derivative should be calculated
1018 //********************************************************************************************************
1019 std::complex<double> eshelbySoluEllipticIntegralsEllipsoid::cIij_i(const double a[3], std::complex<double> la,
1020  std::complex<double> Dla, std::complex<double> cdla_i,
1021  double mult,
1022  ellipticIntegralComponent component )
1023 {
1024  int i = 0, j = 0;
1025 
1026  switch( ( int ) component ){
1027  case _I11_ :
1028  /*default i = 0; j = 0;*/
1029  break;
1030  case _I12_ :
1031  /*default i = 0;*/ j = 1;
1032  break;
1033  case _I13_ :
1034  /*default i = 0;*/ j = 2;
1035  break;
1036  case _I21_ :
1037  i = 1; /*default j = 0;*/
1038  break;
1039  case _I22_ :
1040  i = 1; j = 1;
1041  break;
1042  case _I23_ :
1043  i = 1; j = 2;
1044  break;
1045  case _I31_ :
1046  i = 2; /*default j = 0;*/
1047  break;
1048  case _I32_ :
1049  i = 2; j = 1;
1050  break;
1051  case _I33_ :
1052  i = 2; j = 2;
1053  break;
1054  default:
1055  _errorr( "cIij_i: invalid integral component ");
1056  }//end of switch: component
1057 
1058  return( mult * cdla_i / ( ( SQR( a[i] ) + la ) * ( SQR( a[j] ) + la ) * Dla ) );
1059 }//end of function: cIij_i( )
1060 
1061 //********************************************************************************************************
1062 // description: Function gives all first derivatives of Ferers-Dysons' integrals Ii by using
1063 // the complex differentiation
1064 // last edit: 03. 08. 2010
1065 //-------------------------------------------------------------------------------------------------------
1066 // note: @RedJi[9] - real part of i-th elliptic integral derivative
1067 // @ImdJi[27] - imaginary part of i-th elliptic integral derivative
1068 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
1069 // @x[3] - coordinates of a point
1070 // @Cla - lambda parameter related quantities
1071 //********************************************************************************************************
1072 void eshelbySoluEllipticIntegralsEllipsoid::give_cdIi( double RedJi[9], double ImdJi[27],
1073  const double sort_a[3], cLambda CLa )
1074 {
1075  std::complex<double> aux = 0.;
1076 
1077  //dIi,1
1078  aux = this->cIi_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I1_ );//cI1,1
1079  RedJi[0] = real( aux );//I1,1
1080  ImdJi[0] = imag( aux );//I1,11*_DIFF_H_
1081 
1082  aux = this->cIi_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I2_ );//cI2,1
1083  RedJi[3] = real( aux );//I2,1
1084  ImdJi[9] = imag( aux );//I2,11*_DIFF_H_
1085 
1086  aux = this->cIi_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I3_ );//cI3,1
1087  RedJi[6] = real( aux );//I3,1
1088  ImdJi[18] = imag( aux );//I3,11*_DIFF_H_
1089 
1090  ImdJi[1] = imag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I1_ ) );//I1,12*_DIFF_H_
1091  ImdJi[10] = imag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I2_ ) );//I2,12*_DIFF_H_
1092  ImdJi[19] = imag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I3_ ) );//I3,12*_DIFF_H_
1093 
1094  ImdJi[2] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I1_ ) );//I1,13*_DIFF_H_
1095  ImdJi[11] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I2_ ) );//I2,13*_DIFF_H_
1096  ImdJi[20] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I3_ ) );//I3,13*_DIFF_H_
1097 
1098  //dIi,2
1099  aux = this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I1_ );//cI1,2
1100  RedJi[1] = real( aux );//I1,2
1101  ImdJi[4] = imag( aux );//I1,22*_DIFF_H_
1102 
1103  aux = this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I2_ );//cI2,2
1104  RedJi[4] = real( aux );//I2,2
1105  ImdJi[13] = imag( aux );//I2,22*_DIFF_H_
1106 
1107  aux = this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I3_ );//cI3,2
1108  RedJi[7] = real( aux );//I3,2
1109  ImdJi[22] = imag( aux );//I3,22*_DIFF_H_
1110 
1111  ImdJi[3] = ImdJi[1];//I1,21*_DIFF_H_
1112  ImdJi[12] = ImdJi[10];//I2,21*_DIFF_H_
1113  ImdJi[21] = ImdJi[19];//I3,21*_DIFF_H_
1114 
1115  ImdJi[5] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I1_ ) );//I1,23*_DIFF_H_
1116  ImdJi[14] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I2_ ) );//I2,23*_DIFF_H_
1117  ImdJi[23] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I3_ ) );//I3,23*_DIFF_H_
1118 
1119  //dIi,3
1120  aux = this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I1_ );//cI1,3
1121  RedJi[2] = real( aux );//I1,3
1122  ImdJi[8] = imag( aux );//I1,33*_DIFF_H_
1123 
1124  aux = this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I2_ );//cI2,3
1125  RedJi[5] = real( aux );//I2,3
1126  ImdJi[17] = imag( aux );//I2,33*_DIFF_H_
1127 
1128  aux = this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I3_ );//cI3,3
1129  RedJi[8] = real( aux );//I3,3
1130  ImdJi[26] = imag( aux );//I3,33*_DIFF_H_
1131 
1132  ImdJi[6] = ImdJi[2];//I1,31*_DIFF_H_
1133  ImdJi[15] = ImdJi[11];//I2,31*_DIFF_H_
1134  ImdJi[24] = ImdJi[20];//I3,31*_DIFF_H_
1135 
1136  ImdJi[7] = ImdJi[5];//I1,32*_DIFF_H_
1137  ImdJi[16] = ImdJi[14];//I2,32*_DIFF_H_
1138  ImdJi[25] = ImdJi[23];//I3,32*_DIFF_H_
1139 
1140  return;
1141 }//end of function: give_cdIi( )
1142 
1143 //********************************************************************************************************
1144 // description: Function gives all first derivatives of Ferers-Dysons' integrals Iij by using
1145 // the complex differentiation
1146 // last edit: 09. 08. 2010
1147 //-------------------------------------------------------------------------------------------------------
1148 // note: @dJij[27] - table of elliptic integral derivatives to be calculated
1149 // @cdJij[27] - auxiliary table of all elliptic integral derivatives (needful for Jij_ij)
1150 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
1151 // @Cla - lambda parameter related quantities
1152 //********************************************************************************************************
1153 void eshelbySoluEllipticIntegralsEllipsoid::give_cdIij( double RedJij[27], double ImdJij[81],
1154  const double sort_a[3], cLambda CLa )
1155 {
1156  std::complex<double> aux = 0.;
1157 
1158  //dIij,1
1159  aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I11_ );//cI11,1
1160  RedJij[0] = real( aux );//I11,1
1161  ImdJij[0] = imag( aux );//I11,11*_DIFF_H_
1162 
1163  aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I12_ );//cI12,1
1164  RedJij[1] = real( aux );//I12,1
1165  ImdJij[9] = imag( aux );//I12,11*_DIFF_H_
1166 
1167  aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I13_ );//cI13,1
1168  RedJij[2] = real( aux );//I13,1
1169  ImdJij[18] = imag( aux );//I13,11*_DIFF_H_
1170 
1171  ImdJij[1] = imag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I11_ ) );//I11,12*_DIFF_H_
1172  ImdJij[10] = imag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I12_ ) );//I12,12*_DIFF_H_
1173  ImdJij[19] = imag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I13_ ) );//I13,12*_DIFF_H_
1174 
1175  ImdJij[2] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I11_ ) );//I11,13*_DIFF_H_
1176  ImdJij[11] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I12_ ) );//I12,13*_DIFF_H_
1177  ImdJij[20] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I13_ ) );//I13,13*_DIFF_H_
1178 
1179  RedJij[3] = RedJij[1];//I21,1
1180  ImdJij[27] = ImdJij[9];//I21,11*_DIFF_H_
1181 
1182  aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I22_ );//cI22,1
1183  RedJij[4] = real( aux );//I22,1
1184  ImdJij[36] = imag( aux );//I22,11*_DIFF_H_
1185 
1186  ImdJij[37] = imag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I22_ ) );//I22,12*_DIFF_H_
1187  ImdJij[38] = imag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I22_ ) );//I22,13*_DIFF_H_
1188 
1189  aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I23_ );//cI23,1
1190  RedJij[5] = real( aux );//I23,1
1191  ImdJij[45] = imag( aux );//I23,11*_DIFF_H_
1192 
1193  ImdJij[46] = imag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I23_ ) );//I23,12*_DIFF_H_
1194  ImdJij[47] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I23_ ) );//I23,13*_DIFF_H_
1195 
1196  RedJij[6] = RedJij[2];//I31,1
1197  ImdJij[54] = ImdJij[18];//I31,11*_DIFF_H_
1198 
1199  RedJij[7] = RedJij[5];//I32,1
1200  ImdJij[63] = ImdJij[45];//I32,11*_DIFF_H_
1201 
1202  aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I33_ );//cI33,1
1203  RedJij[8] = real( aux );//I33,1
1204  ImdJij[72] = imag( aux );//I33,11*_DIFF_H_
1205 
1206  ImdJij[37] = imag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I22_ ) );//I22,12*_DIFF_H_
1207  ImdJij[38] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I22_ ) );//I22,13*_DIFF_H_
1208 
1209  //dIij,2
1210  aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I11_ );//cI11,2
1211  RedJij[9] = real( aux );//I11,2
1212  ImdJij[4] = imag( aux );//I11,22*_DIFF_H_
1213 
1214  aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I12_ );//cI12,2
1215  RedJij[10] = real( aux );//I12,2
1216  ImdJij[13] = imag( aux );//I12,22*_DIFF_H_
1217 
1218  aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I13_ );//cI13,2
1219  RedJij[11] = real( aux );//I13,2
1220  ImdJij[22] = imag( aux );//I13,22*_DIFF_H_
1221 
1222  ImdJij[5] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I11_ ) );//I11,23*_DIFF_H_
1223  ImdJij[14] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I12_ ) );//I12,23*_DIFF_H_
1224  ImdJij[23] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I13_ ) );//I13,23*_DIFF_H_
1225 
1226  ImdJij[8] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I11_ ) );//I11,33*_DIFF_H_
1227  ImdJij[17] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I12_ ) );//I12,33*_DIFF_H_
1228  ImdJij[26] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I13_ ) );//I13,33*_DIFF_H_
1229 
1230  RedJij[12] = RedJij[10];//I21,2
1231  ImdJij[31] = ImdJij[13];//I21,22*_DIFF_H_
1232 
1233  aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I22_ );//cI22,2
1234  RedJij[13] = real( aux );//I22,2
1235  ImdJij[40] = imag( aux );//I22,22*_DIFF_H_
1236 
1237  ImdJij[41] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I22_ ) );//I22,23*_DIFF_H_
1238  ImdJij[44] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I22_ ) );//I22,33*_DIFF_H_
1239 
1240  aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I23_ );//cI23,2
1241  RedJij[14] = real( aux );//I23,2
1242  ImdJij[49] = imag( aux );//I23,22*_DIFF_H_
1243 
1244  ImdJij[50] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I23_ ) );//I23,23*_DIFF_H_
1245 
1246  RedJij[15] = RedJij[11];//I31,2
1247  ImdJij[58] = ImdJij[22];//I31,22*_DIFF_H_
1248 
1249  RedJij[16] = RedJij[14];//I32,2
1250  ImdJij[67] = ImdJij[49];//I32,22*_DIFF_H_
1251 
1252  aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I33_ );//cI33,2
1253  RedJij[17] = real( aux );//I33,2
1254  ImdJij[76] = imag( aux );//I33,22*_DIFF_H_
1255 
1256  ImdJij[77] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I33_ ) );//I33,23*_DIFF_H_
1257 
1258  //dIij,3
1259  aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I11_ );//cI11,3
1260  RedJij[18] = real( aux );//I11,3
1261  ImdJij[8] = imag( aux );//I11,33*_DIFF_H_
1262 
1263  aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I12_ );//cI12,3
1264  RedJij[19] = real( aux );//I12,3
1265  ImdJij[17] = imag( aux );//I12,33*_DIFF_H_
1266 
1267  aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I13_ );//cI13,3
1268  RedJij[20] = real( aux );//I13,3
1269  ImdJij[26] = imag( aux );//I13,33*_DIFF_H_
1270 
1271  RedJij[21] = RedJij[19];//I21,3
1272  ImdJij[35] = ImdJij[17];//I21,33*_DIFF_H_
1273 
1274  aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I22_ );//cI22,3
1275  RedJij[22] = real( aux );//I22,3
1276  ImdJij[44] = imag( aux );//I22,33*_DIFF_H_
1277 
1278  aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I23_ );//cI23,3
1279  RedJij[23] = real( aux );//I23,3
1280  ImdJij[53] = imag( aux );//I23,33*_DIFF_H_
1281 
1282  RedJij[24] = RedJij[20];//I31,3
1283  ImdJij[62] = ImdJij[26];//I31,33*_DIFF_H_
1284 
1285  RedJij[25] = RedJij[23];//I32,3
1286  ImdJij[71] = ImdJij[53];//I32,33*_DIFF_H_
1287 
1288  aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I33_ );//cI33,3
1289  RedJij[26] = real( aux );//I33,3
1290  ImdJij[80] = imag( aux );//I33,33*_DIFF_H_
1291 
1292  ImdJij[73] = imag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I33_ ) );//I33,12*_DIFF_H_
1293  ImdJij[74] = imag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I33_ ) );//I33,13*_DIFF_H_
1294 
1295  //initializing of remaining components
1296  ImdJij[3] = ImdJij[1];//I11,21*_DIFF_H_
1297  ImdJij[6] = ImdJij[2];//I11,31*_DIFF_H_
1298  ImdJij[7] = ImdJij[5];//I11,32*_DIFF_H_
1299  ImdJij[12] = ImdJij[10];//I12,21*_DIFF_H_
1300  ImdJij[15] = ImdJij[11];//I12,31*_DIFF_H_
1301  ImdJij[16] = ImdJij[14];//I12,32*_DIFF_H_
1302  ImdJij[21] = ImdJij[19];//I13,21*_DIFF_H_
1303  ImdJij[24] = ImdJij[20];//I13,31*_DIFF_H_
1304  ImdJij[25] = ImdJij[23];//I13,32*_DIFF_H_
1305  ImdJij[28] = ImdJij[10];//I21,12*_DIFF_H_
1306  ImdJij[29] = ImdJij[11];//I21,13*_DIFF_H_
1307  ImdJij[30] = ImdJij[12];//I21,21*_DIFF_H_
1308  ImdJij[32] = ImdJij[14];//I21,23*_DIFF_H_
1309  ImdJij[33] = ImdJij[15];//I21,31*_DIFF_H_
1310  ImdJij[34] = ImdJij[16];//I21,32*_DIFF_H_
1311  ImdJij[39] = ImdJij[37];//I22,21*_DIFF_H_
1312  ImdJij[42] = ImdJij[38];//I22,31*_DIFF_H_
1313  ImdJij[43] = ImdJij[41];//I22,32*_DIFF_H_
1314  ImdJij[48] = ImdJij[46];//I23,21*_DIFF_H_
1315  ImdJij[51] = ImdJij[47];//I23,31*_DIFF_H_
1316  ImdJij[52] = ImdJij[50];//I23,32*_DIFF_H_
1317  ImdJij[55] = ImdJij[19];//I31,12*_DIFF_H_
1318  ImdJij[56] = ImdJij[20];//I31,13*_DIFF_H_
1319  ImdJij[57] = ImdJij[21];//I31,21*_DIFF_H_
1320  ImdJij[59] = ImdJij[23];//I31,23*_DIFF_H_
1321  ImdJij[60] = ImdJij[24];//I31,31*_DIFF_H_
1322  ImdJij[61] = ImdJij[25];//I31,32*_DIFF_H_
1323  ImdJij[64] = ImdJij[46];//I32,12*_DIFF_H_
1324  ImdJij[65] = ImdJij[47];//I32,13*_DIFF_H_
1325  ImdJij[66] = ImdJij[48];//I32,21*_DIFF_H_
1326  ImdJij[68] = ImdJij[50];//I32,23*_DIFF_H_
1327  ImdJij[69] = ImdJij[51];//I32,31*_DIFF_H_
1328  ImdJij[70] = ImdJij[52];//I32,32*_DIFF_H_
1329  ImdJij[75] = ImdJij[73];//I33,21*_DIFF_H_
1330  ImdJij[78] = ImdJij[74];//I33,31*_DIFF_H_
1331  ImdJij[79] = ImdJij[77];//I33,32*_DIFF_H_
1332 
1333  return;
1334 }//end of function: give_cdIij( )
1335 
1336 //********************************************************************************************************
1337 // description: Function gives all second derivatives of Ferers-Dysons' integrals Ii by using
1338 // the complex differentiation
1339 // last edit: 04. 08. 2010
1340 //-------------------------------------------------------------------------------------------------------
1341 // note: @ddJi[27] - table of second derivatives of elliptic integrals to be calculated
1342 // @ImdJi[27] - table of the imaginary parts of the first complex derivatives of
1343 // elliptic integrals Ii
1344 //********************************************************************************************************
1345 void eshelbySoluEllipticIntegralsEllipsoid::give_cddIi( double ddJi[27], double ImdJi[27] )
1346 {
1347  for( int i = 0; i < 27; i++ ){
1348  ddJi[i] = ImdJi[i]/_DIFF_H_;
1349  }
1350 
1351  return;
1352 }
1353 
1354 //********************************************************************************************************
1355 // description: Function gives all second derivatives of Ferers-Dysons' integrals Iij by using
1356 // the complex differentiation
1357 // last edit: 05. 08. 2010
1358 //-------------------------------------------------------------------------------------------------------
1359 // note: @ddJij[81] - table of second derivatives of elliptic integrals to be calculated
1360 // @ImdJij[81] - table of the imaginary parts of the first complex derivatives of
1361 // elliptic integrals Iij
1362 //********************************************************************************************************
1363 void eshelbySoluEllipticIntegralsEllipsoid::give_cddIij( double ddJij[81], double ImdJij[81] )
1364 {
1365  for( int i = 0; i < 81; i++ ){
1366  ddJij[i] = ImdJij[i]/_DIFF_H_;
1367  }
1368 
1369  return;
1370 }
1371 
1372 
1373 //Numerical derivatives of Ferers-Dysons integrals:
1374 //---------------------------------------------------------
1375 
1376 //********************************************************************************************************
1377 // description: Function gives the first numerical derivative of Ii_i in 'i' direction
1378 // last edit: 26. 08. 2010
1379 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
1380 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
1381 // @la - lambda parameter (not a Lame constant !!!)
1382 // @Dla - Delta lambda
1383 // @dla_i - first derivative of lambda in 'i' direction
1384 // @mult - is the multiplier calculated as: -2. * PI * a[0] * a[1] * a[2]
1385 // @component - component which derivative should be calculated
1386 //********************************************************************************************************
1387 double eshelbySoluEllipticIntegralsEllipsoid::nIi_i(const double a[3], double la, double Dla, double dla_i,
1388  double mult, ellipticIntegralComponent component )
1389 {
1390  int comp = component - 1;
1391 
1392  return( mult * dla_i / ( ( SQR( a[comp] ) + la ) * Dla ) );
1393 }
1394 
1395 //********************************************************************************************************
1396 // description: Function gives the first numerical derivative of Iij_i in 'i' direction
1397 // last edit: 27. 08. 2010
1398 // note: Solution adopted from Toshio Mura's book (1982), page(s) 73
1399 // @a[3] - dimensions of an ellipsoid (a1, a2, a3 semi-axies dimensions)
1400 // @la - lambda parameter (not a Lame constant !!!)
1401 // @Dla - Delta lambda
1402 // @dla_i - first derivative of lambda in 'i' direction
1403 // @mult - is the multiplier calculated as: -2. * PI * a[0] * a[1] * a[2]
1404 // @component - component which derivative should be calculated
1405 //********************************************************************************************************
1406 double eshelbySoluEllipticIntegralsEllipsoid::nIij_i(const double a[3], double la, double Dla, double dla_i,
1407  double mult, ellipticIntegralComponent component )
1408 {
1409  int i = 0, j = 0;
1410 
1411  switch( ( int ) component ){
1412  case _I11_ :
1413  /*default i = 0; j = 0;*/
1414  break;
1415  case _I12_ :
1416  /*default i = 0;*/ j = 1;
1417  break;
1418  case _I13_ :
1419  /*default i = 0;*/ j = 2;
1420  break;
1421  case _I21_ :
1422  i = 1; /*default j = 0;*/
1423  break;
1424  case _I22_ :
1425  i = 1; j = 1;
1426  break;
1427  case _I23_ :
1428  i = 1; j = 2;
1429  break;
1430  case _I31_ :
1431  i = 2; /*default j = 0;*/
1432  break;
1433  case _I32_ :
1434  i = 2; j = 1;
1435  break;
1436  case _I33_ :
1437  i = 2; j = 2;
1438  break;
1439  default:
1440  _errorr( "nIij_i: invalid integral component ");
1441  }//end of switch: component
1442 
1443  return( mult * dla_i / ( ( SQR( a[i] ) + la ) * ( SQR( a[j] ) + la ) * Dla ) );
1444 }//end of function: nIij_i( )
1445 
1446 //********************************************************************************************************
1447 // description: Function gives all first derivatives of Ferers-Dysons' integrals Ii by using
1448 // numerical differentiation
1449 // last edit: 27. 08. 2010
1450 //-------------------------------------------------------------------------------------------------------
1451 // note: @dJi[9] - i-th elliptic integral derivative
1452 // @ddJih[27] - second eliptic integral derivative * h
1453 // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
1454 // @x[3] - coordinates of a point
1455 // @nla - lambda parameter related quantities
1456 //********************************************************************************************************
1457 void eshelbySoluEllipticIntegralsEllipsoid::give_ndIi( double dJi[9], double ddJih[27],
1458  const double sort_a[3], nLambda nLa )
1459 {
1460  //dIi,1
1461  dJi[0] = this->nIi_i( sort_a, nLa.la, nLa.Dla, nLa.dla, nLa.mult, _I1_ );//I1,1
1462  ddJih[0] = this->nIi_i( sort_a, nLa.la1, nLa.Dla1, nLa.dla_11, nLa.mult, _I1_ ) - dJi[0];//I1,11*h
1463 
1464  dJi[3] = this->nIi_i( sort_a, nLa.la, nLa.Dla, nLa.dla, nLa.mult, _I2_ );//I2,1
1465  ddJih[9] = this->nIi_i( sort_a, nLa.la1, nLa.Dla1, nLa.dla_11, nLa.mult, _I2_ ) - dJi[3];//I2,11*h
1466 
1467  dJi[6] = this->nIi_i( sort_a, nLa.la, nLa.Dla, nLa.dla, nLa.mult, _I3_ );//I3,1
1468  ddJih[18] = this->nIi_i( sort_a, nLa.la1, nLa.Dla1, nLa.dla_11, nLa.mult, _I3_ ) - dJi[6];//I3,11*h
1469 
1470 
1471  // ddJih[1] = this->nIi_i( sort_a, nLa.la2, nLa.Dla2, nLa.dla_12, nLa.mult, _I1_ ) - dJi[3];//I1,12*h
1472  //cimag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I1_ ) );//I1,12*h
1473  // ImdJi[10] = cimag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I2_ ) );//I2,12*_DIFF_H_
1474  // ImdJi[19] = cimag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I3_ ) );//I3,12*_DIFF_H_
1475 
1476  // ImdJi[2] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I1_ ) );//I1,13*_DIFF_H_
1477  // ImdJi[11] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I2_ ) );//I2,13*_DIFF_H_
1478  // ImdJi[20] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I3_ ) );//I3,13*_DIFF_H_
1479 
1480  // //dIi,2
1481  // aux = this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I1_ );//cI1,2
1482  // RedJi[1] = creal( aux );//I1,2
1483  // ImdJi[4] = cimag( aux );//I1,22*_DIFF_H_
1484 
1485  // aux = this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I2_ );//cI2,2
1486  // RedJi[4] = creal( aux );//I2,2
1487  // ImdJi[13] = cimag( aux );//I2,22*_DIFF_H_
1488 
1489  // aux = this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I3_ );//cI3,2
1490  // RedJi[7] = creal( aux );//I3,2
1491  // ImdJi[22] = cimag( aux );//I3,22*_DIFF_H_
1492 
1493  // ImdJi[3] = ImdJi[1];//I1,21*_DIFF_H_
1494  // ImdJi[12] = ImdJi[10];//I2,21*_DIFF_H_
1495  // ImdJi[21] = ImdJi[19];//I3,21*_DIFF_H_
1496 
1497  // ImdJi[5] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I1_ ) );//I1,23*_DIFF_H_
1498  // ImdJi[14] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I2_ ) );//I2,23*_DIFF_H_
1499  // ImdJi[23] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I3_ ) );//I3,23*_DIFF_H_
1500 
1501  // //dIi,3
1502  // aux = this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I1_ );//cI1,3
1503  // RedJi[2] = creal( aux );//I1,3
1504  // ImdJi[8] = cimag( aux );//I1,33*_DIFF_H_
1505 
1506  // aux = this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I2_ );//cI2,3
1507  // RedJi[5] = creal( aux );//I2,3
1508  // ImdJi[17] = cimag( aux );//I2,33*_DIFF_H_
1509 
1510  // aux = this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I3_ );//cI3,3
1511  // RedJi[8] = creal( aux );//I3,3
1512  // ImdJi[26] = cimag( aux );//I3,33*_DIFF_H_
1513 
1514  // ImdJi[6] = ImdJi[2];//I1,31*_DIFF_H_
1515  // ImdJi[15] = ImdJi[11];//I2,31*_DIFF_H_
1516  // ImdJi[24] = ImdJi[20];//I3,31*_DIFF_H_
1517 
1518  // ImdJi[7] = ImdJi[5];//I1,32*_DIFF_H_
1519  // ImdJi[16] = ImdJi[14];//I2,32*_DIFF_H_
1520  // ImdJi[25] = ImdJi[23];//I3,32*_DIFF_H_
1521 
1522  return;
1523 }//end of function: give_ndIi( )
1524 
1525 //********************************************************************************************************
1526 // // description: Function gives all first derivatives of Ferers-Dysons' integrals Iij by using
1527 // // the complex differentiation
1528 // // last edit: 09. 08. 2010
1529 // //-------------------------------------------------------------------------------------------------------
1530 // // note: @dJij[27] - table of elliptic integral derivatives to be calculated
1531 // // @cdJij[27] - auxiliary table of all elliptic integral derivatives (needful for Jij_ij)
1532 // // @sort_a[3] - sorted semiaxes' dimensoins of an ellipsoid (a1>a2>a3)
1533 // // @Cla - lambda parameter related quantities
1534 // //********************************************************************************************************
1535 // void eshelbySoluEllipticIntegralsEllipsoid::give_ndIij( double RedJij[27], double ImdJij[81],
1536 // const double sort_a[3], cLambda CLa )
1537 // {
1538 // complex double aux = 0.;
1539 
1540 // //dIij,1
1541 // aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I11_ );//cI11,1
1542 // RedJij[0] = creal( aux );//I11,1
1543 // ImdJij[0] = cimag( aux );//I11,11*_DIFF_H_
1544 
1545 // aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I12_ );//cI12,1
1546 // RedJij[1] = creal( aux );//I12,1
1547 // ImdJij[9] = cimag( aux );//I12,11*_DIFF_H_
1548 
1549 // aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I13_ );//cI13,1
1550 // RedJij[2] = creal( aux );//I13,1
1551 // ImdJij[18] = cimag( aux );//I13,11*_DIFF_H_
1552 
1553 // ImdJij[1] = cimag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I11_ ) );//I11,12*_DIFF_H_
1554 // ImdJij[10] = cimag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I12_ ) );//I12,12*_DIFF_H_
1555 // ImdJij[19] = cimag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I13_ ) );//I13,12*_DIFF_H_
1556 
1557 // ImdJij[2] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I11_ ) );//I11,13*_DIFF_H_
1558 // ImdJij[11] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I12_ ) );//I12,13*_DIFF_H_
1559 // ImdJij[20] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I13_ ) );//I13,13*_DIFF_H_
1560 
1561 // RedJij[3] = RedJij[1];//I21,1
1562 // ImdJij[27] = ImdJij[9];//I21,11*_DIFF_H_
1563 
1564 // aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I22_ );//cI22,1
1565 // RedJij[4] = creal( aux );//I22,1
1566 // ImdJij[36] = cimag( aux );//I22,11*_DIFF_H_
1567 
1568 // ImdJij[37] = cimag( this->cIi_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I22_ ) );//I22,12*_DIFF_H_
1569 // ImdJij[38] = cimag( this->cIi_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I22_ ) );//I22,13*_DIFF_H_
1570 
1571 // aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I23_ );//cI23,1
1572 // RedJij[5] = creal( aux );//I23,1
1573 // ImdJij[45] = cimag( aux );//I23,11*_DIFF_H_
1574 
1575 // ImdJij[46] = cimag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I23_ ) );//I23,12*_DIFF_H_
1576 // ImdJij[47] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I23_ ) );//I23,13*_DIFF_H_
1577 
1578 // RedJij[6] = RedJij[2];//I31,1
1579 // ImdJij[54] = ImdJij[18];//I31,11*_DIFF_H_
1580 
1581 // RedJij[7] = RedJij[5];//I32,1
1582 // ImdJij[63] = ImdJij[45];//I32,11*_DIFF_H_
1583 
1584 // aux = this->cIij_i( sort_a, CLa.cla1, CLa.cDla1, CLa.cdla_11, CLa.mult, _I33_ );//cI33,1
1585 // RedJij[8] = creal( aux );//I33,1
1586 // ImdJij[72] = cimag( aux );//I33,11*_DIFF_H_
1587 
1588 // ImdJij[37] = cimag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I22_ ) );//I22,12*_DIFF_H_
1589 // ImdJij[38] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I22_ ) );//I22,13*_DIFF_H_
1590 
1591 // //dIij,2
1592 // aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I11_ );//cI11,2
1593 // RedJij[9] = creal( aux );//I11,2
1594 // ImdJij[4] = cimag( aux );//I11,22*_DIFF_H_
1595 
1596 // aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I12_ );//cI12,2
1597 // RedJij[10] = creal( aux );//I12,2
1598 // ImdJij[13] = cimag( aux );//I12,22*_DIFF_H_
1599 
1600 // aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I13_ );//cI13,2
1601 // RedJij[11] = creal( aux );//I13,2
1602 // ImdJij[22] = cimag( aux );//I13,22*_DIFF_H_
1603 
1604 // ImdJij[5] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I11_ ) );//I11,23*_DIFF_H_
1605 // ImdJij[14] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I12_ ) );//I12,23*_DIFF_H_
1606 // ImdJij[23] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I13_ ) );//I13,23*_DIFF_H_
1607 
1608 // ImdJij[8] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I11_ ) );//I11,33*_DIFF_H_
1609 // ImdJij[17] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I12_ ) );//I12,33*_DIFF_H_
1610 // ImdJij[26] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I13_ ) );//I13,33*_DIFF_H_
1611 
1612 // RedJij[12] = RedJij[10];//I21,2
1613 // ImdJij[31] = ImdJij[13];//I21,22*_DIFF_H_
1614 
1615 // aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I22_ );//cI22,2
1616 // RedJij[13] = creal( aux );//I22,2
1617 // ImdJij[40] = cimag( aux );//I22,22*_DIFF_H_
1618 
1619 // ImdJij[41] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I22_ ) );//I22,23*_DIFF_H_
1620 // ImdJij[44] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I22_ ) );//I22,33*_DIFF_H_
1621 
1622 // aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I23_ );//cI23,2
1623 // RedJij[14] = creal( aux );//I23,2
1624 // ImdJij[49] = cimag( aux );//I23,22*_DIFF_H_
1625 
1626 // ImdJij[50] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I23_ ) );//I23,23*_DIFF_H_
1627 
1628 // RedJij[15] = RedJij[11];//I31,2
1629 // ImdJij[58] = ImdJij[22];//I31,22*_DIFF_H_
1630 
1631 // RedJij[16] = RedJij[14];//I32,2
1632 // ImdJij[67] = ImdJij[49];//I32,22*_DIFF_H_
1633 
1634 // aux = this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_22, CLa.mult, _I33_ );//cI33,2
1635 // RedJij[17] = creal( aux );//I33,2
1636 // ImdJij[76] = cimag( aux );//I33,22*_DIFF_H_
1637 
1638 // ImdJij[77] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_23, CLa.mult, _I33_ ) );//I33,23*_DIFF_H_
1639 
1640 // //dIij,3
1641 // aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I11_ );//cI11,3
1642 // RedJij[18] = creal( aux );//I11,3
1643 // ImdJij[8] = cimag( aux );//I11,33*_DIFF_H_
1644 
1645 // aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I12_ );//cI12,3
1646 // RedJij[19] = creal( aux );//I12,3
1647 // ImdJij[17] = cimag( aux );//I12,33*_DIFF_H_
1648 
1649 // aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I13_ );//cI13,3
1650 // RedJij[20] = creal( aux );//I13,3
1651 // ImdJij[26] = cimag( aux );//I13,33*_DIFF_H_
1652 
1653 // RedJij[21] = RedJij[19];//I21,3
1654 // ImdJij[35] = ImdJij[17];//I21,33*_DIFF_H_
1655 
1656 // aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I22_ );//cI22,3
1657 // RedJij[22] = creal( aux );//I22,3
1658 // ImdJij[44] = cimag( aux );//I22,33*_DIFF_H_
1659 
1660 // aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I23_ );//cI23,3
1661 // RedJij[23] = creal( aux );//I23,3
1662 // ImdJij[53] = cimag( aux );//I23,33*_DIFF_H_
1663 
1664 // RedJij[24] = RedJij[20];//I31,3
1665 // ImdJij[62] = ImdJij[26];//I31,33*_DIFF_H_
1666 
1667 // RedJij[25] = RedJij[23];//I32,3
1668 // ImdJij[71] = ImdJij[53];//I32,33*_DIFF_H_
1669 
1670 // aux = this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_33, CLa.mult, _I33_ );//cI33,3
1671 // RedJij[26] = creal( aux );//I33,3
1672 // ImdJij[80] = cimag( aux );//I33,33*_DIFF_H_
1673 
1674 // ImdJij[73] = cimag( this->cIij_i( sort_a, CLa.cla2, CLa.cDla2, CLa.cdla_12, CLa.mult, _I33_ ) );//I33,12*_DIFF_H_
1675 // ImdJij[74] = cimag( this->cIij_i( sort_a, CLa.cla3, CLa.cDla3, CLa.cdla_13, CLa.mult, _I33_ ) );//I33,13*_DIFF_H_
1676 
1677 // //initializing of remaining components
1678 // ImdJij[3] = ImdJij[1];//I11,21*_DIFF_H_
1679 // ImdJij[6] = ImdJij[2];//I11,31*_DIFF_H_
1680 // ImdJij[7] = ImdJij[5];//I11,32*_DIFF_H_
1681 // ImdJij[12] = ImdJij[10];//I12,21*_DIFF_H_
1682 // ImdJij[15] = ImdJij[11];//I12,31*_DIFF_H_
1683 // ImdJij[16] = ImdJij[14];//I12,32*_DIFF_H_
1684 // ImdJij[21] = ImdJij[19];//I13,21*_DIFF_H_
1685 // ImdJij[24] = ImdJij[20];//I13,31*_DIFF_H_
1686 // ImdJij[25] = ImdJij[23];//I13,32*_DIFF_H_
1687 // ImdJij[28] = ImdJij[10];//I21,12*_DIFF_H_
1688 // ImdJij[29] = ImdJij[11];//I21,13*_DIFF_H_
1689 // ImdJij[30] = ImdJij[12];//I21,21*_DIFF_H_
1690 // ImdJij[32] = ImdJij[14];//I21,23*_DIFF_H_
1691 // ImdJij[33] = ImdJij[15];//I21,31*_DIFF_H_
1692 // ImdJij[34] = ImdJij[16];//I21,32*_DIFF_H_
1693 // ImdJij[39] = ImdJij[37];//I22,21*_DIFF_H_
1694 // ImdJij[42] = ImdJij[38];//I22,31*_DIFF_H_
1695 // ImdJij[43] = ImdJij[41];//I22,32*_DIFF_H_
1696 // ImdJij[48] = ImdJij[46];//I23,21*_DIFF_H_
1697 // ImdJij[51] = ImdJij[47];//I23,31*_DIFF_H_
1698 // ImdJij[52] = ImdJij[50];//I23,32*_DIFF_H_
1699 // ImdJij[55] = ImdJij[19];//I31,12*_DIFF_H_
1700 // ImdJij[56] = ImdJij[20];//I31,13*_DIFF_H_
1701 // ImdJij[57] = ImdJij[21];//I31,21*_DIFF_H_
1702 // ImdJij[59] = ImdJij[23];//I31,23*_DIFF_H_
1703 // ImdJij[60] = ImdJij[24];//I31,31*_DIFF_H_
1704 // ImdJij[61] = ImdJij[25];//I31,32*_DIFF_H_
1705 // ImdJij[64] = ImdJij[46];//I32,12*_DIFF_H_
1706 // ImdJij[65] = ImdJij[47];//I32,13*_DIFF_H_
1707 // ImdJij[66] = ImdJij[48];//I32,21*_DIFF_H_
1708 // ImdJij[68] = ImdJij[50];//I32,23*_DIFF_H_
1709 // ImdJij[69] = ImdJij[51];//I32,31*_DIFF_H_
1710 // ImdJij[70] = ImdJij[52];//I32,32*_DIFF_H_
1711 // ImdJij[75] = ImdJij[73];//I33,21*_DIFF_H_
1712 // ImdJij[78] = ImdJij[74];//I33,31*_DIFF_H_
1713 // ImdJij[79] = ImdJij[77];//I33,32*_DIFF_H_
1714 
1715 // return;
1716 // }//end of function: give_ndIij( )
1717 
1718 // Lambda computed according to Mura p. 86.
1719 double eshelbySoluEllipticIntegralsEllipsoid::getLambda(const double a[3], double x1, double x2, double x3)
1720 {
1721  double A = x1 * x1;
1722  double B = x2 * x2;
1723  double C = x3 * x3;
1724  double D = a[0]*a[0];
1725  double E = a[1]*a[1];
1726  double F = a[2]*a[2];
1727 
1728  // Coeficients of cubic equation. Obtained
1729  double aa = 1.0;
1730  double bb = D + E + F - A - B - C;
1731  double cc = D*E + D*F + E*F - A*E - A*F - B*D - B*F - C*D - C*E;
1732  double dd = D*E*F - A*E*F - B*D*F - C*D*E;
1733 
1734  // Discriminant of cubic equation. See http://en.wikipedia.org/wiki/Cubic_function#Roots_of_a_cubic_function
1735  double Discriminant = 18.0*aa*bb*cc*dd - 4.0*pow(bb, 3.0)*dd + SQR(bb)*SQR(cc) - 4.0*aa*pow(cc, 3.0) - 27.0*SQR(aa)*SQR(dd);
1736  if (Discriminant <= 0.0)
1737  {
1738  _errorr("Negative of zero discriminant in getLambda().");
1739  }
1740 
1741  // Write the depressed form of cubid equation. See http://en.wikipedia.org/wiki/Cubic_function#Reduction_to_a_depressed_cubic
1742  double pp = (3.0 * aa * cc - SQR(bb)) / 3.0 / SQR(aa);
1743  double qq = (2.0 * pow(bb, 3.0) - 9.0 * aa * bb * cc + 27.0 * SQR(aa) * dd) / 27.0 / pow(aa, 3.0);
1744 
1745  if (pp > 0.0)
1746  {
1747  _errorr("Coeficient pp has to be positive in getLambda().");
1748  }
1749 
1750  // Roots of depresed cubic equations. See http://en.wikipedia.org/wiki/Cubic_function#Three_real_roots
1751  double t0 = 2.0 * sqrt(-pp / 3.0) * cos(1.0 / 3.0 *acos(3.0 * qq / 2.0 / pp * sqrt(-3.0 / pp)) - 0.0 * 2.0 * PI / 3.0);
1752  double t1 = 2.0 * sqrt(-pp / 3.0) * cos(1.0 / 3.0 *acos(3.0 * qq / 2.0 / pp * sqrt(-3.0 / pp)) - 1.0 * 2.0 * PI / 3.0);
1753  double t2 = 2.0 * sqrt(-pp / 3.0) * cos(1.0 / 3.0 *acos(3.0 * qq / 2.0 / pp * sqrt(-3.0 / pp)) - 2.0 * 2.0 * PI / 3.0);
1754 
1755  // Roots of cubic equation in general (original) form.
1756  double root0 = t0 - bb / 3.0;
1757  double root1 = t1 - bb / 3.0;
1758  double root2 = t2 - bb / 3.0;
1759 
1760  // Return the highest of the roots.
1761  return MAX(MAX(root0, root1), root2);
1762 }
1763 
1764 } // end of namespace mumech
1765 
1766 /*end of file*/
mumech::cLambda::cdla_13
std::complex< double > cdla_13
Definition: types.h:688
mumech::eshelbySoluEllipticIntegralsEllipsoid::eI
double eI(double a[3], double k, double Theta, double mult)
Definition: esei_Ellipsoid.cpp:288
MAX
#define MAX(a, b)
Definition: macros.h:101
mumech::cLambda::cdla_22
std::complex< double > cdla_22
Definition: types.h:689
mumech::nLambda::Dla
double Dla
Definition: types.h:707
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_ndIi
void give_ndIi(double dJi[9], double ddJih[27], const double sort_a[3], nLambda nLa)
Definition: esei_Ellipsoid.cpp:1457
mumech::Point::loc_x
double loc_x[3]
Local coordinates of a point.
Definition: matrix.h:137
mumech::eshelbySoluEllipticIntegralsEllipsoid::I33
double I33(const double a[3], double I32, double I31, double lambda, double mult, double Dla, bool intpoint)
Definition: esei_Ellipsoid.cpp:434
esei_Ellipsoid.h
Class eshelbySoluEllipticIntegralsEllipsoid.
mumech::cLambda::cx2
std::complex< double > cx2[3]
Definition: types.h:678
mumech::cLambda::cDla2
std::complex< double > cDla2
Definition: types.h:684
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_dIi
void give_dIi(double dJi[9], const double sort_a[3], const double x[3], aLambda aLa)
Definition: esei_Ellipsoid.cpp:725
mumech::eshelbySoluEllipticIntegralsEllipsoid::I2
double I2(double I1, double I3, double lambda, double mult, double Dla, bool intpoint)
Definition: esei_Ellipsoid.cpp:327
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_Ip_Derivative
double give_Ip_Derivative(aLambda aLa, derivativeDirection p)
Definition: esei_Ellipsoid.cpp:556
mumech::legendreIntegrals
Class legendreIntegrals provides functions calculating values of elliptical integrals.
Definition: legendreIntegrals.h:38
mumech::cLambda::cdla_12
std::complex< double > cdla_12
Definition: types.h:687
legendreIntegrals.h
Class legendreIntegrals.
mumech::Inclusion::P
const Problem * P
problem description
Definition: inclusion.h:69
PI
#define PI
Definition: macros.h:71
mumech::eshelbySoluEllipticIntegrals::I
const InclusionRecord3D * I
Definition: esei.h:45
mumech::cLambda::cDla3
std::complex< double > cDla3
Definition: types.h:685
mumech::nLambda::la
double la
Definition: types.h:706
mumech::eshelbySoluEllipticIntegralsEllipsoid::I12
double I12(const double a[3], double I1, double I2)
Definition: esei_Ellipsoid.cpp:384
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_Iip_Derivative
double give_Iip_Derivative(const double a[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection p)
Definition: esei_Ellipsoid.cpp:572
mumech::eshelbySoluEllipticIntegralsEllipsoid::I3
double I3(const double a[3], double k, double Theta, double lambda, double mult, bool intpoint)
Definition: esei_Ellipsoid.cpp:344
mumech::eshelbySoluEllipticIntegralsEllipsoid::cIij_i
std::complex< double > cIij_i(const double a[3], std::complex< double > la, std::complex< double > Dla, std::complex< double > cdla_i, double mult, ellipticIntegralComponent component)
Definition: esei_Ellipsoid.cpp:1019
mumech::legendreIntegrals::ellf
static double ellf(double phi, double ak)
Function calculating the Legendre&#39;s integral of the first kind - numerical recepies.
Definition: legendreIntegrals.cpp:49
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_Iijp_Derivative
double give_Iijp_Derivative(const double a[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection p)
Definition: esei_Ellipsoid.cpp:601
NUM_PERTURB
#define NUM_PERTURB
Definition: types.h:732
mumech::aLambda::la
double la
Definition: types.h:697
mumech::derivativeDirection
derivativeDirection
Definition: types.h:208
mumech::eshelbySoluEllipticIntegralsEllipsoid::nIi_i
double nIi_i(const double a[3], double la, double Dla, double dla_i, double mult, ellipticIntegralComponent component)
Definition: esei_Ellipsoid.cpp:1387
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_cddIi
void give_cddIi(double ddJi[27], double ImdJi[27])
Definition: esei_Ellipsoid.cpp:1345
mumech::eshelbySoluLambda::giveLambdaDerivative
double giveLambdaDerivative(const double a[3], const double x[3], double lambda, derivativeDirection direction)
Definition: eshelbySoluLambda.cpp:167
mumech::ellipticIntegralComponent
ellipticIntegralComponent
Definition: types.h:198
macros.h
The header file of usefull macros.
mumech::eshelbySoluEllipticIntegrals::getPerturbatedLambdas
void getPerturbatedLambdas(double *lambdas, const double loc_x[3])
Helper function.
Definition: esei.cpp:51
elementaryFunctions.h
Collection of the functions of basic manipulations, some of them can be used instead of their counter...
mumech::eshelbySoluEllipticIntegralsEllipsoid::I1
double I1(const double a[3], double k, double Theta, double mult)
Definition: esei_Ellipsoid.cpp:305
mumech::Point::ddJij
double ddJij[81]
Second derivatives of elliptic integral Iij.
Definition: matrix.h:145
mumech::Problem::give_diffType
DiffTypes give_diffType(void) const
Give type of ...
Definition: problem.h:336
mumech::Point
Single Point data structure - contribution from Single inclusion.
Definition: matrix.h:133
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_cdIi
void give_cdIi(double RedJi[9], double ImdJi[27], const double sort_a[3], cLambda CLa)
Definition: esei_Ellipsoid.cpp:1072
mumech::cLambda::cx1
std::complex< double > cx1[3]
Definition: types.h:677
mumech::eshelbySoluEllipticIntegralsEllipsoid::I11
double I11(const double a[3], double I12, double I13, double lambda, double mult, double Dla, bool intpoint)
Definition: esei_Ellipsoid.cpp:371
mumech::eshelbySoluEllipticIntegralsEllipsoid::I23
double I23(const double a[3], double I2, double I3)
Definition: esei_Ellipsoid.cpp:421
mumech::Inclusion::ndiff_1
double ndiff_1
derivative step for the first derivations
Definition: inclusion.h:108
mumech::eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegralFirstDerivative
double giveEllipticIntegralFirstDerivative(const double a[3], const double x[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection direction)
Definition: esei_Ellipsoid.cpp:489
mumech::eshelbySoluLambda::giveLambda
static double giveLambda(const double a[3], const double x[3], InclusionGeometry shape)
Definition: eshelbySoluLambda.cpp:70
mumech::cLambda::cdla_33
std::complex< double > cdla_33
Definition: types.h:691
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_ddIi
void give_ddIi(double ddJi[27], const double sort_a[3], const double x[3], aLambda aLa)
Definition: esei_Ellipsoid.cpp:807
mumech::nLambda::dla
double dla
Definition: types.h:708
mumech::nLambda::Dla1
double Dla1
Definition: types.h:715
mumech::eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegralSeccondDerivative
double giveEllipticIntegralSeccondDerivative(const double a[3], const double x[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection direction_1, derivativeDirection direction_2)
Definition: esei_Ellipsoid.cpp:525
_DIFF_H_
#define _DIFF_H_
Definition: types.h:666
mumech::eshelbySoluEllipticIntegralsEllipsoid::giveDerivativesOfEllipticIntegrals
void giveDerivativesOfEllipticIntegrals(Point *point, bool intpoint)
Function gives the values of Ferers-Dyson&#39;s elliptic integral derivatives of the inclusion this->I...
Definition: esei_Ellipsoid.cpp:103
SQR
#define SQR(a)
Definition: macros.h:97
mumech::Point::dJi
double dJi[9]
First derivatives of elliptic integral Ii.
Definition: matrix.h:142
mumech::Inclusion::a
double * a
Inclusion semiaxes&#39; dimensions in global arrangement.
Definition: inclusion.h:76
mumech::Point::ddJi
double ddJi[27]
Second derivatives of elliptic integral Ii.
Definition: matrix.h:144
mumech::cLambda::cla2
std::complex< double > cla2
Definition: types.h:681
mumech::eshelbySoluLambda
Class of functions calculating the values of &#39;lambda&#39; parameter and its derivatives.
Definition: eshelbySoluLambda.h:36
_errorr
#define _errorr(_1)
Definition: gelib.h:160
mumech::nLambda::la1
double la1
Definition: types.h:712
mumech::eshelbySoluEllipticIntegralsEllipsoid::I22
double I22(const double a[3], double I21, double I23, double lambda, double mult, double Dla, bool intpoint)
Definition: esei_Ellipsoid.cpp:408
_delta_lambda
#define _delta_lambda(lambda)
Definition: esei_Ellipsoid.cpp:45
mumech::nLambda::dla_11
double dla_11
Definition: types.h:718
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_cddIij
void give_cddIij(double ddJij[81], double ImdJij[81])
Definition: esei_Ellipsoid.cpp:1363
mumech::aLambda::ddla
double ddla[10]
Definition: types.h:700
mumech::nLambda::mult
double mult
Definition: types.h:728
mumech::Point::dJij
double dJij[27]
First derivatives of elliptic integral Iij.
Definition: matrix.h:143
mumech::eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegralDerivative
double giveEllipticIntegralDerivative(const double a[3], const double x[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection direction_1, derivativeDirection direction_2)
Definition: esei_Ellipsoid.cpp:457
mumech::eshelbySoluLambda::c_dLambda
std::complex< double > c_dLambda(const double a[3], std::complex< double > x[3], std::complex< double > lambda, derivativeDirection direction)
Definition: eshelbySoluLambda.cpp:286
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_cdIij
void give_cdIij(double RedJij[27], double ImdJij[81], const double sort_a[3], cLambda CLa)
Definition: esei_Ellipsoid.cpp:1153
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_ddIij
void give_ddIij(double ddJi[27], const double sort_a[3], const double x[3], aLambda aLa)
Definition: esei_Ellipsoid.cpp:861
mumech::aLambda::Dla
double Dla
Definition: types.h:698
mumech::eshelbySoluEllipticIntegralsEllipsoid::getLambda
double getLambda(const double a[3], double x1, double x2, double x3)
Returns lambda for a given point (x1, x2, x3)
Definition: esei_Ellipsoid.cpp:1719
mumech::eshelbySoluEllipticIntegralsEllipsoid::cIi_i
std::complex< double > cIi_i(const double a[3], std::complex< double > la, std::complex< double > Dla, std::complex< double > cdla_i, double mult, ellipticIntegralComponent component)
Definition: esei_Ellipsoid.cpp:998
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_Iijpq_Derivative
double give_Iijpq_Derivative(const double a[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection p, derivativeDirection q)
Definition: esei_Ellipsoid.cpp:636
mumech::aLambda::mult
double mult
Definition: types.h:701
mumech::cLambda::cx3
std::complex< double > cx3[3]
Definition: types.h:679
mumech::legendreIntegrals::elle
static double elle(double phi, double ak)
Function calculating the Legendre&#39;s integral of the seccond kind - numerical recepies.
Definition: legendreIntegrals.cpp:63
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_Iipq_Derivative
double give_Iipq_Derivative(const double a[3], aLambda aLa, ellipticIntegralComponent component, derivativeDirection p, derivativeDirection q)
Definition: esei_Ellipsoid.cpp:686
mumech::aLambda::dla
double dla[4]
Definition: types.h:699
mumech::eshelbySoluEllipticIntegralsEllipsoid::I13
double I13(const double a[3], double I1, double I3)
Definition: esei_Ellipsoid.cpp:395
_4PI_
#define _4PI_
Definition: types.h:50
mumech::cLambda::cdla_23
std::complex< double > cdla_23
Definition: types.h:690
mumech::eshelbySoluEllipticIntegralsEllipsoid::nIij_i
double nIij_i(const double a[3], double la, double Dla, double dla_i, double mult, ellipticIntegralComponent component)
Definition: esei_Ellipsoid.cpp:1406
mumech::cLambda::cla1
std::complex< double > cla1
Definition: types.h:680
mumech::eshelbySoluEllipticIntegralsEllipsoid::give_dIij
void give_dIij(double dJij[27], const double sort_a[3], const double x[3], aLambda aLa)
Definition: esei_Ellipsoid.cpp:753
mumech::cLambda::cla3
std::complex< double > cla3
Definition: types.h:682
mumech::cLambda::cdla_11
std::complex< double > cdla_11
Definition: types.h:686
problem.h
Class Problem.
mumech::cLambda::cDla1
std::complex< double > cDla1
Definition: types.h:683
mumech::cLambda::mult
double mult
Definition: types.h:692
mumech::eshelbySoluEllipticIntegralsEllipsoid::giveEllipticIntegrals
void giveEllipticIntegrals(double J[13], double lambda, bool intpoint)
Function gives the values of Ferers-Dyson&#39;s elliptic integrals of the inclusion this->I.
Definition: esei_Ellipsoid.cpp:59
This reference manual is part of the muMECH documentation. Generated on Wed Feb 24 2016 22:50:07 for muMECH by doxygen 1.8.11.