"""Computation of the global geometric factor"""
import numpy as np
from ...matlab_funcs import lscov, legendre
[docs]def coeff2ggf(coeff, poisson_ratio=.45):
"""Compute the global geometric factor from stress coefficients
The radial displacements of an elastic sphere can be expressed in
terms of Legendre polynomials (see :cite:`Lure1964` equation 6.2.9)
whose coefficients are computed from the Legendre decomposition of
the radial stress.
Notes
-----
- For a :math:`\\sigma_0 \cos^n(\\theta)` stress profile,
the GGF already includes the peak stress according to:
.. math::
\\text{GGF} = \\sigma_0 F_\\text{G}.
- This is a conversion of the Matlab script GGF.m to Python. The
code solves a linear system of equations to determine all
Legendre coefficients. The new implementation in
:func:`ggf.legendre2ggf` uses the direct solution and
thus should be preferred.
"""
## Parameters used in the program
Theta1 =0
# Nu =input('enter the value of nu (Poisson s ratio) \n')
Nu = poisson_ratio
##
SigmaMat= coeff
n=SigmaMat.size-1
SigmaMatamended = SigmaMat/2 # amendement due to the correction of the paper by Ananthakrishnan
Sigmatot=np.zeros((2*n+2,1))
for k in range(n):
Sigmatot[k]=SigmaMatamended[k]
#Sigmatot[k+n+1] = 0
Mat=np.zeros((2*n+2,2*n+2))
r=1# r is set to 1 for the calculation.
N=n
for d in range(n+1):
Mat[d,d]=(r**d)*(d + 1)*(d**2-d-2-2*Nu)#
Mat[d,d+N+1]=(r**(d - 2))*d*(d - 1)#
for d in range(2,n):
Mat[d+N,d+1]=(r**(d+1))*(-1+2*(d+1)+(d+1)**2+2*Nu)*(-(d+1))*(d+2)/(1+2*(d+1)) # i pour le coefficient et j pour la ligne de sigma pour le coefficient P(d)
Mat[d+N,d-1]=(r**(d-1))*(-1+2*(d-1)+((d-1)**2)+2*Nu)*(d)*(d-1)/(1+2*(d-1))
Mat[d+N,d+N+1+1]=(r**(d-1))*(d)*(-(d+1))*(d+2)/(1+2*(d+1))
Mat[d+N,d+N]=(r**(d-3))*(d-2)*(d)*(d-1)/(1+2*(d-1))
d=0
Mat[d+N,d+1]=(r**(d+1))*(-1+2*(d+1)+(d+1)**2+2*Nu)*(-(d+1))*(d+2)/(1+2*(d+1))
Mat[d+N,d+N+1+1]=(r**(d-1))*(d)*(-(d+1))*(d+2)/(1+2*(d+1))
d=1
Mat[d+N,d+1]=(r**(d+1))*(-1+2*(d+1)+((d+1)**2)+2*Nu)*(-(d+1))*(d+2)/(1+2*(d+1))
Mat[d+N,d+N+1+1]=(r**(d-1))*(d)*(-(d+1))*(d+2)/(1+2*(d+1))
d=n
Mat[d+N,d-1]=(r**(d-1))*(-1+2*(d-1)+((d-1)**2)+2*Nu)*(d)*(d-1)/(1+2*(d-1))
Mat[d+N,d+N]=(r**(d-3))*(d-2)*(d)*(d-1)/(1+2*(d-1))
d=n+1
Mat[d+N,d-1]=(r**(d-1))*(-1+2*(d-1)+((d-1)**2)+2*Nu)*(d)*(d-1)/(1+2*(d-1))
Mat[d+N,d+N]=(r**(d-3))*(d-2)*(d)*(d-1)/(1+2*(d-1))
Mat1=Mat
sol=lscov(Mat1,Sigmatot)
PL = np.zeros(n+1, dtype=float)
w0 = np.zeros(n+1, dtype=float)
## Calculation of the Global Geometrical Factor including sigma0
ww=0
for k in range(0,n):
PLA=legendre(k,np.cos(Theta1))
PL[k]=np.real_if_close(PLA[0][0]) # legendre polynomial with m=0!
w0[k]=((sol[k])*(r**(k+1))*(k+1 )*(k-2 +4*Nu) +(sol[n+k+1])*r**(k-1)*k)*PL[k]
ww=ww+w0[k]
GF=ww
return GF