<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Poroplastic Data Fitting
Step1:  Summary
In the cells to follow, the following material parameters were found
$$\begin{align}
B_0 &= 14617807286.8\
B_1 &= 40384983097.2\
B_2 &= 385649437.858\
P_0 & = −164761936.257 \
P_1 & = 3.20119273834e−10\
P_2 & = 7.39166987894e−18\
P_3 & = 0.0983914345654\
G_1 & = 9647335534.93 \
G_2 & = 2.3838775292e−09 \
G_3 & = −7.40942609805e−07\
\end{align}$$
Read in the Data
Read in the hydrostatic data and compute derived values.
Step2:  Hydrostatic Response
Elastic Unloading Curve
Plot the pressure vs. volume strain curve and determine the section in which elastic unloading occurs
Step3:  It appears that the unloading occurs at data point 101 and continues until the end of the data.  This curve will be used to fit the bulk modulus parameters.  Below, scipy is used to optimize the parameters to the curve.
Step4:  Poro response
With the bulk response determined, find the porosity parameters
Step5:  Shear Response
Step6:  The axial stress versus axial strain plot shows that the response is linear, meaning that the elastic modulus is constant.
Step7:  The shear modulus can now be determined
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from numpy import *
from bokeh import *
from bokeh.plotting import *
output_notebook()
from matmodlab2 import *
from pandas import read_excel
from scipy.optimize import leastsq
diff = lambda x: np.ediff1d(x, to_begin=0.)
trace = lambda x, s='SIG': x[s+'11'] + x[s+'22'] + x[s+'33']
RTJ2 = lambda x: sqrt(((x['SIG11']-x['SIG22'])**2 + 
                       (x['SIG22']-x['SIG33'])**2 + 
                       (x['SIG33']-x['SIG22'])**2)/6.)
Explanation: Poroplastic Data Fitting
End of explanation
df = read_excel('porodata.xlsx', sheetname='hydrostatic')
df['EV'] = trace(df, 'STRAIN')
df['I1'] = trace(df, 'SIG')
df['dEV'] = diff(df['EV'])
df['dI1'] = diff(df['I1'])
Explanation: Summary
In the cells to follow, the following material parameters were found
$$\begin{align}
B_0 &= 14617807286.8\
B_1 &= 40384983097.2\
B_2 &= 385649437.858\
P_0 & = −164761936.257 \
P_1 & = 3.20119273834e−10\
P_2 & = 7.39166987894e−18\
P_3 & = 0.0983914345654\
G_1 & = 9647335534.93 \
G_2 & = 2.3838775292e−09 \
G_3 & = −7.40942609805e−07\
\end{align}$$
Read in the Data
Read in the hydrostatic data and compute derived values.
End of explanation
plot = figure(x_axis_label='Volume Strain', y_axis_label='Pressure')
plot.circle(-df['EV'], -df['I1']/3.)
plot.text(-df['EV'], -df['I1']/3.,
    text=range(len(df)),text_color="#333333",
    text_align="left", text_font_size="5pt")
show(plot)
Explanation: Hydrostatic Response
Elastic Unloading Curve
Plot the pressure vs. volume strain curve and determine the section in which elastic unloading occurs
End of explanation
kfun = lambda B0, B1, B2, I1: B0 + B1 * exp(-B2 / abs(I1))
def kmm_bulk(x, fac, I1, K):
    B0, B1, B2 = x * fac
    return K - kfun(B0, B1, B2, I1)
imax = 101
df1 = df.iloc[imax:].copy()
K = np.array(df1['dI1'] / 3. / df1['dEV'])
b0 = np.array((K[-1], K[0] - K[-1], 1e9))
fac = 1e9
B, icov = leastsq(kmm_bulk, b0/fac, args=(fac, df1['I1'], K))
B0, B1, B2 = B * fac
B0, B1, B2
plot = figure(x_axis_label='Bulk Modulus', y_axis_label='Pressure')
plot.circle(-df1['I1']/3., K)
plot.line(-df['I1']/3., kfun(B0, B1, B2, df['I1']), color='red')
show(plot)
Explanation: It appears that the unloading occurs at data point 101 and continues until the end of the data.  This curve will be used to fit the bulk modulus parameters.  Below, scipy is used to optimize the parameters to the curve.
End of explanation
df['EP'] = df['I1'] / 3. / kfun(B0, B1, B2, df['I1']) - df['EV']
p3 = max(df['EP'])
df['PORO'] = p3 - df['EP']
plot = figure(x_axis_label='Plastic Strain', y_axis_label='Pressure')
plot.circle(df['EP'], -df['I1']/3.)
show(plot)
plot = figure(x_axis_label='Pressure', y_axis_label='PORO')
df2 = df.iloc[:imax].copy()
plot.circle(-df2['I1']/3., df2['PORO'])
show(plot)
def pfun(P0, P1, P2, P3, I1):
    xi = -I1 / 3. + P0
    return P3 * exp(-(P1 + P2 * xi) * xi)
    
def kmm_poro(x, fac, I1, P):
    p0, p1, p2, p3 = asarray(x) * fac
    return P - pfun(p0, p1, p2, p3, I1)
p0 = (1, 1, 1, p3)
fac = np.array([1e8, 1e-10, 1e-18, 1])
p, icov = leastsq(kmm_poro, p0, args=(fac, df2['I1'], df2['PORO']))
P0, P1, P2, P3 = p * fac
P0, P1, P2, P3
plot = figure(x_axis_label='Pressure', y_axis_label='PORO')
plot.circle(-df2['I1']/3., df2['PORO'], legend='Data')
plot.line(-df2['I1']/3., pfun(P0, P1, P2, P3, df2['I1']), color='red', legend='Fit')
show(plot)
Explanation: Poro response
With the bulk response determined, find the porosity parameters
End of explanation
keys = (2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 22.5, 30.0)
colors = ('red', 'blue', 'orange', 'purple', 
          'green', 'black', 'magenta', 'teal', 'cyan')
df2 = {}
p = figure(x_axis_label='I1', y_axis_label='Sqrt[J2]')
p1 = figure(x_axis_label='Axial Strain', y_axis_label='Axial Stress')
for (i, key) in enumerate(keys):
    key = 'txc p={0:.01f}MPa'.format(key)
    x = read_excel('porodata.xlsx', sheetname=key)
    x['I1'] = trace(x, 'SIG')
    x['RTJ2'] = RTJ2(x)
    df2[key] = x
    p.circle(-df2[key]['I1'], df2[key]['RTJ2'], legend=key[4:], color=colors[i])
    
    # determine where hydrostatic preload ends
    j = nonzero(x['SIG11'] - x['SIG22'])[0]
    E0, S0 = df2[key]['STRAIN11'][j[0]], df2[key]['SIG11'][j[0]]
    p1.circle(-df2[key]['STRAIN11'][j]+E0, -df2[key]['SIG11'][j]+S0,
              legend=key[4:], color=colors[i])
p.legend.orientation = 'horizontal'
show(p1)
show(p)
Explanation: Shear Response
End of explanation
key = 'txc p=2.5MPa'
j = nonzero(df2[key]['SIG11'] - df2[key]['SIG22'])[0]
df3 = df2[key].iloc[j].copy()
E0, S0 = df3['STRAIN11'].iloc[0], df3['SIG11'].iloc[0]
EF, SF = df3['STRAIN11'].iloc[-1], df3['SIG11'].iloc[-1]
E = (SF - S0) / (EF - E0)
print('{0:E}'.format(E))
Explanation: The axial stress versus axial strain plot shows that the response is linear, meaning that the elastic modulus is constant.
End of explanation
G = lambda I1: 3 * kfun(B0, B1, B2, I1) * E / (9 * kfun(B0, B1, B2, I1) - E)
gfun = lambda g0, g1, g2, rtj2: g0 * (1 - g1 * exp(-g2 * rtj2)) / (1 - g1)
def kmm_shear(x, fac, rtj2, G):
    g0, g1, g2 = asarray(x) * fac
    return G - gfun(g0, g1, g2, rtj2)
g = asarray(G(df3['I1']))
g0 = (g[0], .0001, 0)
fac = 1.
g, icov = leastsq(kmm_shear, g0, args=(fac, RTJ2(df3), g))
G0, G1, G2 = g * fac
G0, G1, G2
p2 = figure(x_axis_label='Sqrt[J2]', y_axis_label='Shear Modulus')
p2.circle(RTJ2(df3), G(df3['I1']))
p2.line(RTJ2(df3), gfun(G0, G1, G2, RTJ2(df3)), color='red')
show(p2)
Explanation: The shear modulus can now be determined
End of explanation
<END_TASK>