Last update: August 2018, added extended elastic impedance code and various helper functions.
Connolly introduced the concept of elastic impedance ($EI$) in 1999 as an elegant solution to tie the interpretation of far angle stacks to well data (near-offsets can be tied to well data through simple acoustic impedance).
He rearranged the Aki-Richards equations so that the output is not reflectivity but something like a proper log, which is strikingly similar in shape to acoustic impedance ($AI$), until something interesting happens (like a hydrocarbon-bearing sand) and then elastic impedance shows a different response from $AI$.
The formula looks like this:
$$ EI = V_p^a V_s^b \rho^c $$where:
$$ a = 1 + \tan^2 \theta $$$$ b = -8 k \sin^2 \theta $$$$ c = 1 - 4 k \sin^2 \theta $$$$ k = \left( \frac{V_s}{V_p} \right)^2 $$The $k$ factor is usually kept constant in a given area.
Connolly also presents an alternative form which "can be used for generating synthetics" (his own words, check the Appendix in Connolly, 1999):
$$ EI = V_p \rho^* $$$$ \rho^\ast = V_p^{(\tan^2 \theta)} V_s^b \rho^c ) $$In practice we have $EI$ expressed as the product of velocity and a pseudo-density ($ \rho^* $) that is angle-dependant. So any kind of 1D modeling can be used to generate angle-dependant synthetic seismograms, and the usual calibration of $ V_p $ through time-depth tables ensures the correct positioning in time.
Whitcome in 2002 introduced the normalized $EI$ function to allow comparisons between $EI$ logs calculated at different angles:
$$ EI^{NORM} = \rho_0 V_{p0} \, \left[ \left( \frac{V_p}{V_{p0}}\right)^a \left( \frac{V_s}{V_{s0}} \right)^b \left( \frac{\rho}{\rho_0} \right)^c\right]$$Where the scaling factors $V_{p0}$ $V_{s0}$ $\rho_0$ are averages of the input elastic logs computed over the range of interest or the values at the top of the target zone.
I discuss some practical aspects of using elastic impedance in the context of generating synthetic seismograms and also show how I reproduced a couple of figures from "Seismic Amplitude" by Simm & Bacon in this notebook here: "The relationship between reflectivity and elastic impedance", or how to reproduce figure 5.62 from Seismic Amplitude by Simm & Bacon (2014).
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.signal import detrend
from ipywidgets import interact, fixed
import ipywidgets as widgets
%matplotlib inline
# comment out the following if you're not on a Mac with HiDPI display
%config InlineBackend.figure_format = 'retina'
def calc_ei(vp, vs, rho, alpha, scal=None, k=0.25):
alpha = np.radians(alpha)
a = 1 + (np.tan(alpha)) ** 2
b = -8 * k * ((np.sin(alpha)) ** 2)
c = 1 - 4 * k * ((np.sin(alpha)) ** 2)
if scal is None:
ei = vp**a * vs**b * rho** c
else:
vp0, vs0, rho0 = scal[0], scal[1], scal[2]
ei = vp0*rho0 * ( (vp/vp0) ** a * (vs/vs0) ** b * (rho/rho0) ** c)
return ei
w2=pd.read_csv('qsiwell2_frm.csv', index_col=0)
z1, z2 = 2100, 2150
ff = (w2.index>=z1) & (w2.index<=z2)
scal=w2.VP[ff].mean(), w2.VS[ff].mean(), w2.RHO[ff].mean()
ei = calc_ei(w2.VP, w2.VS, w2.RHO,30)
ei_norm = calc_ei(w2.VP, w2.VS, w2.RHO,30,scal)
f, ax = plt.subplots(nrows=1,ncols=3,sharey=True,figsize=(8,6))
ax[0].plot(gr,z,'k', label='GR')
ax[1].plot(ei,z,'r', label='EI 30')
ax[2].plot(ei_norm,z,'r', label='EI 30 norm')
ax[2].plot(ip,z,'k', label='AI')
ax[0].set_ylim(2250,2100)
ax[0].set_xlabel('[API]')
ax[1].set_xlabel('[m/s*g/cc]')
ax[2].set_xlabel('[m/s*g/cc]')
ax[0].set_ylabel('Depth [m]')
ax[1].locator_params(nbins=5, axis='x')
for aa in ax:
aa.grid(which='both')
aa.legend()
ax[0].set_title('Lithology')
ax[1].set_title('EI')
ax[2].set_title('AI and Normalized EI')
f.suptitle('Comparison Acoustic vs Elastic Impedance at 30 degrees', fontsize=16);
Just as an additional (and potentially useful) piece of information, here's how you would specify the constants discussed above ($k$ plus the scaling factors $V_{p0}$ $V_{s0}$ $\rho_0$) in the well-tie module ("Well Seismic Fusion") in Halliburton's Decision Space:
from IPython.display import Image
Image(filename='DecisionSpace_welltie_EI_Calc.png', width=400)
I haven't found further details on this but it looks like that you have indeed to specify the Vp/Vs ratio and not the $k$ constant, so given $Vp/Vs = 2$ then $ k = (Vs/Vp)^2 = 0.25 $.
The normalization constants are those that I call scaling factors.
def calc_eei(vp, vs, rho, chi, scal=None, k=0.25):
'''
EI (C) aadm 2018
extended elastic impedance
'''
chi = np.radians(chi)
p = np.cos(chi) + np.sin(chi)
q = -8 * k * np.sin(chi)
r = np.cos(chi) - 4 * k * np.sin(chi)
if scal is None:
vp0, vs0, rho0 = np.nanmean(vp), np.nanmean(vs), np.nanmean(rho)
else:
vp0, vs0, rho0 = scal[0], scal[1], scal[2]
eei = vp0*rho0 * ( (vp/vp0) ** p * (vs/vs0) ** q * (rho/rho0) ** r)
return eei
def eei_scaling(wells,names,tops,start,window,brine=False):
tmp0=np.zeros((len(wells),3))
tmp1=np.zeros((len(wells)))
if brine:
idvp,idvs,idrho,fluidtype=['VP_FRMB','VS_FRMB','RHO_FRMB','brine']
else:
idvp,idvs,idrho,fluidtype=['VP','VS','RHO','insitu']
for i,ww in enumerate(wells):
z1 = tops[names[i]][start]
z2 = tops[names[i]][start]+window
ww = ww.loc[(ww.index>=z1) & (ww.index<=z2)]
k = (np.nanmean(ww[idvs]) / np.nanmean(ww[idvp])) ** 2 # avg over interval of interest
scal=np.nanmean(np.column_stack((ww[idvp], ww[idvs], ww[idrho])), axis=0)
tmp0[i,:]=scal
tmp1[i]=k
print('{:.<12s} vp0={:<6.0f} vs0={:<6.0f} rho0={:<6.2f} vp/vs={:<6.2f} k={:<6.2f} ({:s})'.format(names[i],scal[0],scal[1],scal[2],(1/k)**.5,k,fluidtype))
vp0,vs0,rho0=tmp0.mean(axis=0)
k0=tmp1.mean()
print('{:.<12s} vp0={:<6.0f} vs0={:<6.0f} rho0={:<6.2f} vp/vs={:<6.2f} k={:<6.2f} ({:s})'.format('[AVG]',vp0,vs0,rho0,(1/k0)**.5,k0,fluidtype))
return tmp0.mean(axis=0),tmp1.mean()
def eei_diagnostic(ww,ztop,zbot,scal=None,k=None,ln=False,chi_min=None,chi_max=None,brine=False,targetlog='GR',removetrend=False,name=None):
label_eei = 'EEI'
if ln:
label_eei = 'ln('+label_eei+')'
if brine:
idvp,idvs,idrho=['VP_FRMB','VS_FRMB','RHO_FRMB']
else:
idvp,idvs,idrho=['VP','VS','RHO']
ww = ww.loc[(ww.index>=ztop) & (ww.index<=zbot)]
ww = ww.interpolate(limit_direction='both')
if chi_min is None:
chi_min, chi_max = -90,90
chi=np.linspace(chi_min,chi_max)
corr_eei = np.zeros(chi.size)
eei_array = np.zeros((ww.index.size,chi.size))
target=ww[targetlog]
if k is None:
k = (np.nanmean(ww[idvs]) / np.nanmean(ww[idvp]))**2
if scal is None:
scal = np.nanmean(np.column_stack((ww[idvp], ww[idvs], ww[idrho])), axis=0)
for n,X in enumerate(chi):
eei = calc_eei(ww[idvp],ww[idvs],ww[idrho],X,scal=scal,k=k)
if ln:
eei = np.log(eei)
if removetrend:
eei = pd.Series(data=detrend(eei), index=ww.index)
target = pd.Series(data=detrend(target), index=ww.index)
corr_eei[n]=target.corr(eei)
eei_array[:,n]=eei.values
maxcorr_eei=np.max(np.abs(corr_eei))
maxcorr_eei_id=np.argmax(np.abs(corr_eei))
print('max corr. EEI-{:s}={:.2f} @ X={:.0f}'.format(targetlog,maxcorr_eei,chi[maxcorr_eei_id]))
fig = plt.subplots(figsize=(10,5))
ax0 = plt.subplot2grid((1,5), (0,0), colspan=1)
ax1 = plt.subplot2grid((1,5), (0,1), colspan=1)
ax2 = plt.subplot2grid((1,5), (0,2), colspan=3)
axtarget = ax0.twiny()
axtarget.plot(target, ww.index, lw=1, color='.4')
axtarget.set_xlabel(targetlog, color='0.4')
axtarget.tick_params(axis='x', colors='0.4')
axip = ax1.twiny()
axip.plot(ww[idvp]*ww[idrho], ww.index, lw=1, color='black')
axip.set_xlabel('Ip')
ax2.plot(chi,corr_eei, ls='-', color='black')
ax2.set_xlabel('CHI angle')
ax2.set_ylabel('cross-correlation')
ax2.grid()
ax2.set_ylim(-1,1)
ax2.set_xlim(chi.min(),chi.max())
ax2.yaxis.set_label_position('right'), ax2.yaxis.tick_right()
if removetrend:
ax2.set_title('(de-trended)')
ax0.set_ylabel('DEPTH (m MD)')
ax1.set_yticklabels([])
for aa in [ax0,ax1]:
aa.plot(eei_array[:,maxcorr_eei_id], ww.index, lw=1, color='red')
aa.set_xlabel('{:s}{:.0f}'.format(label_eei,chi[maxcorr_eei_id]), color='red')
aa.tick_params(axis='x', colors='red')
aa.set_ylim(zbot,ztop)
aa.grid()
if name is not None:
plt.suptitle(name, weight='bold')
return eei_array,corr_eei
def eei_optimalchi(wells,names,tops,start,window,scal,k,ln=False,chi_min=None,chi_max=None,brine=False,litholog='VSH',fluidlog='SWE',removetrend=False):
label_eei = 'EEI'
if ln:
label_eei = 'ln('+label_eei+')'
if removetrend:
label_eei = 'de-trended '+label_eei
if brine:
idvp,idvs,idrho=['VP_FRMB','VS_FRMB','RHO_FRMB']
else:
idvp,idvs,idrho=['VP','VS','RHO']
if chi_min is None:
chi_min, chi_max = -90,90
chi=np.linspace(chi_min,chi_max)
corr_litho, corr_fluid = (np.zeros(chi.size) for _ in range(2))
fig,ax=plt.subplots(nrows=1,ncols=2,figsize=(12,6))
for i,ww in enumerate(wells):
z1 = tops[names[i]][start]
z2 = tops[names[i]][start]+window
ww = ww.loc[(ww.index>=z1) & (ww.index<=z2)]
ww = ww.interpolate(limit_direction='both')
litho,fluid=ww[litholog],ww[fluidlog]
for n,X in enumerate(chi):
eei = calc_eei(ww[idvp],ww[idvs],ww[idrho],X,scal=scal,k=k)
if ln:
eei = np.log(eei)
if removetrend:
eei = pd.Series(data=detrend(eei), index=ww.index)
litho = pd.Series(data=detrend(litho), index=ww.index)
fluid = pd.Series(data=detrend(fluid), index=ww.index)
corr_litho[n]=litho.corr(eei)
corr_fluid[n]=fluid.corr(eei)
maxcorr_litho = np.max(np.abs(corr_litho))
maxcorr_litho_id=np.argmax(np.abs(corr_litho))
maxcorr_fluid = np.max(np.abs(corr_fluid))
maxcorr_fluid_id=np.argmax(np.abs(corr_fluid))
print('{:s} Litho: max corr. EEI={:.2f} @ X={:.0f}'.format(names[i],maxcorr_litho,chi[maxcorr_litho_id]))
print('{:s} Fluid: max corr. EEI={:.2f} @ X={:.0f}'.format(names[i],maxcorr_fluid,chi[maxcorr_fluid_id]))
ax[0].plot(chi,corr_litho, label=names[i])
ax[1].plot(chi,corr_fluid, label=names[i])
ax[0].set_title('{}-{}'.format(label_eei,litholog))
ax[1].set_title('{}-{}'.format(label_eei,fluidlog))
for aa in ax:
aa.grid(b=True)
aa.set_xlabel('CHI angle'), aa.set_ylabel('cross-correlation')
aa.legend()
aa.set_ylim(-1,1)
plt.tight_layout()
def eei_scanlogs(wells,names,tops,start,window,scal,k,ln=False,chi_min=None,chi_max=None,brine=False,litholog='VSH',fluidlog='SWE',removetrend=False):
label_eei = 'EEI'
if ln:
label_eei = 'ln('+label_eei+')'
if brine:
idvp,idvs,idrho=['VP_FRMB','VS_FRMB','RHO_FRMB']
else:
idvp,idvs,idrho=['VP','VS','RHO']
if chi_min is None:
chi_min, chi_max = -90,90
chi=np.linspace(chi_min,chi_max,10)
corr_litho, corr_fluid = (np.zeros(chi.size) for _ in range(2))
for i,ww in enumerate(wells):
z1 = tops[names[i]][start]
z2 = tops[names[i]][start]+window
ww = ww.loc[(ww.index>=z1) & (ww.index<=z2)]
ww = ww.interpolate(limit_direction='both')
eei_array = np.zeros((ww.index.size,chi.size))
litho,fluid=ww[litholog],ww[fluidlog]
for n,X in enumerate(chi):
eei = calc_eei(ww[idvp],ww[idvs],ww[idrho],X,scal=scal,k=k)
if ln:
eei = np.log(eei)
if removetrend:
eei = pd.Series(data=detrend(eei), index=ww.index)
litho = pd.Series(data=detrend(litho), index=ww.index)
fluid = pd.Series(data=detrend(fluid), index=ww.index)
eei_array[:,n] = eei.values
f,ax=plt.subplots(nrows=2,ncols=10,sharey=True,figsize=(12,8))
for hh in range(10):
ax[0,hh].set_xticklabels([])
ax[1,hh].set_xlabel('{:s}{:.0f}'.format(label_eei,chi[hh]), color='red')
axli = ax[0,hh].twiny()
axli.plot(litho, ww.index, lw=1, color='0.4')
axli.set_xlabel(litholog, color='0.4')
axli.tick_params(axis='x', colors='0.4')
axli.set_ylim(z2,z1)
axfl = ax[1,hh].twiny()
axfl.plot(fluid, ww.index, lw=1, color='blue')
axfl.set_xlabel(fluidlog, color='blue')
axfl.tick_params(axis='x', colors='blue')
axfl.set_ylim(z2,z1)
for n in range(2):
ax[n,hh].plot(eei_array[:,hh], ww.index, lw=1, color='red')
ax[n,hh].tick_params(axis='x', colors='red')
ax[0,0].set_ylim(z2,z1)
ax[0,0].set_ylabel('DEPTH (m MD)')
ax[1,0].set_ylabel('DEPTH (m MD)')
if removetrend:
f.suptitle('{} (de-trended)'.format(names[i]))
else:
f.suptitle('{}'.format(names[i]))
plt.subplots_adjust(hspace=.2,wspace=0.2,left=.05,right=.98)
I will define some arbitrary markers to guide the interpretation (actually top_sand
is the top of the Heimdal reservoir, as gathered from Quantitative Seismic Interpretation by Avseth et al., 2005).
wells= [ w2]
names= ['WELL_2']
tmp1={'top_shale': 2078, 'top_sand': 2153, 'owc': 2183, 'base_sand': 2250}
markers=pd.DataFrame.from_dict(tmp1, orient='index', columns=['WELL_2'])
In case of multiple wells I would add the markers for the other wells in this way and have a single Dataframe with all my stratigraphic markers.
tmp2={'top_shale': 3000, 'top_sand': 3100, 'owc': 3130, 'base_sand': 3200}
dummy2=pd.DataFrame.from_dict(tmp2, orient='index', columns=['WELL_X'])
pd.concat([markers,dummy2],axis=1)
WELL_2 | WELL_X | |
---|---|---|
top_shale | 2078 | 3000 |
top_sand | 2155 | 3100 |
owc | 2188 | 3130 |
base_sand | 2250 | 3200 |
scal0,k0=eei_scaling(wells,names,markers,'top_sand',200,brine=False)
WELL_2...... vp0=2956 vs0=1382 rho0=2.19 vp/vs=2.14 k=0.22 (insitu) [AVG]....... vp0=2956 vs0=1382 rho0=2.19 vp/vs=2.14 k=0.22 (insitu)
eei_optimalchi(wells,names,markers,'top_sand',200,scal0,k0,brine=False,removetrend=False)
WELL_2 Litho: max corr. EEI=0.43 @ X=57 WELL_2 Fluid: max corr. EEI=0.60 @ X=17
Thomas & al. (2013) suggest to use the natural logarithm of EEI and perform the correlation analysis on the detrended target and logarithm-EEI curves.
From the quick tests I have run I think that detrending is effective only when using particularly long windows (how long? perhaps >100m?).
Also, the articles states to de-trend the target logs too, specifically mentioning that petrophsical logs "will contain some trend". That may be the case if using a normal acquired log as target like GR for example, but I'm not so sure about water saturation and porosity logs calculated with some more accuracy.
The following is an example of scanning for optimal $\chi$ angles using larger windows, using original, de-trended logs and ln(EEI):
print('input: EEI, trend removal=No')
eei_optimalchi(wells,names,markers,'top_sand',100,scal0,k0,chi_min=-180,chi_max=180,ln=False,brine=False,removetrend=False)
print('input: EEI, trend removal=Yes')
eei_optimalchi(wells,names,markers,'top_sand',100,scal0,k0,chi_min=-180,chi_max=180,ln=False,brine=False,removetrend=True)
print('input: ln(EEI), trend removal=Yes')
eei_optimalchi(wells,names,markers,'top_sand',100,scal0,k0,chi_min=-180,chi_max=180,ln=True,brine=False,removetrend=True)
input: EEI, trend removal=No WELL_2 Litho: max corr. EEI=0.45 @ X=-121 WELL_2 Fluid: max corr. EEI=0.59 @ X=-151 input: EEI, trend removal=Yes WELL_2 Litho: max corr. EEI=0.39 @ X=-107 WELL_2 Fluid: max corr. EEI=0.44 @ X=-136 input: ln(EEI), trend removal=Yes WELL_2 Litho: max corr. EEI=0.39 @ X=70 WELL_2 Fluid: max corr. EEI=0.43 @ X=-143
Finally, the eei_scanlogs
function plots the target logs, i.e. the lithology and fluid logs (VSH
and SWE
respectively) with EEI curves at different angles superimposed.
eei_scanlogs(wells,names,markers,'top_sand',100,scal0,k0,chi_min=-180,chi_max=180,brine=False,removetrend=False)
def eei_optimalchi_interactive(filename,depth_range,scal,k,ln=False,brine=False,litholog='VSH',fluidlog='SWE',removetrend=False):
ww=pd.read_csv(filename, index_col=0)
label_eei = 'EEI'
if ln:
label_eei = 'ln('+label_eei+')'
if removetrend:
label_eei = 'de-trended '+label_eei
if brine:
idvp,idvs,idrho=['VP_FRMB','VS_FRMB','RHO_FRMB']
else:
idvp,idvs,idrho=['VP','VS','RHO']
chi=np.linspace(-180,180)
corr_litho, corr_fluid = (np.zeros(chi.size) for _ in range(2))
fig,ax=plt.subplots(nrows=1,ncols=2,figsize=(12,6))
z1 = depth_range[0]
z2 = depth_range[1]
ww = ww.loc[(ww.index>=z1) & (ww.index<=z2)]
ww = ww.interpolate(limit_direction='both')
litho,fluid=ww[litholog],ww[fluidlog]
for n,X in enumerate(chi):
eei = calc_eei(ww[idvp],ww[idvs],ww[idrho],X,scal=scal,k=k)
if ln:
eei = np.log(eei)
if removetrend:
eei = pd.Series(data=detrend(eei), index=ww.index)
litho = pd.Series(data=detrend(litho), index=ww.index)
fluid = pd.Series(data=detrend(fluid), index=ww.index)
corr_litho[n]=litho.corr(eei)
corr_fluid[n]=fluid.corr(eei)
maxcorr_litho = np.max(np.abs(corr_litho))
maxcorr_litho_id=np.argmax(np.abs(corr_litho))
maxcorr_fluid = np.max(np.abs(corr_fluid))
maxcorr_fluid_id=np.argmax(np.abs(corr_fluid))
print('Litho: max corr. EEI={:.2f} @ X={:.0f}'.format(maxcorr_litho,chi[maxcorr_litho_id]))
print('Fluid: max corr. EEI={:.2f} @ X={:.0f}'.format(maxcorr_fluid,chi[maxcorr_fluid_id]))
ax[0].plot(chi,corr_litho)
ax[1].plot(chi,corr_fluid)
ax[0].set_title('{}-{}'.format(label_eei,litholog))
ax[1].set_title('{}-{}'.format(label_eei,fluidlog))
for aa in ax:
aa.grid(b=True)
aa.set_xlabel('CHI angle'), aa.set_ylabel('cross-correlation')
aa.set_ylim(-1,1)
plt.tight_layout()
interact(eei_optimalchi_interactive,filename='qsiwell2_frm.csv',
depth_range=widgets.FloatRangeSlider(value=[1100, 1300],min=1050,max=1400.0,step=10,description='depth range',
continuous_update=False,orientation='horizontal',readout=True,readout_format='.0f'),
scal=fixed(scal0),k=fixed(k0),
ln=widgets.Checkbox(value=False,description='Logarithm EEI'),
brine=widgets.Checkbox(value=False,description='Use brine-replaced logs'),
removetrend=widgets.Checkbox(value=False,description='Remove trend'),
litholog=fixed('VSH'),fluidlog=fixed('SWE'));
Litho: max corr. EEI=nan @ X=-180 Fluid: max corr. EEI=nan @ X=-180
ggg={'Gamma Ray':'GR','Shale Volumerosity':'PHIE','Water Saturation':'SWE'}
interact(eei_diagnostic,ww=fixed(w2),ztop=(2000,2100,10),zbot=(2300,2400,10),scal=fixed(scal0),k=fixed(k0),
ln=widgets.Checkbox(value=False,description='Logarithm EEI'),
brine=widgets.Checkbox(value=False,description='Use brine-replaced logs'),
removetrend=widgets.Checkbox(value=False,description='Remove trend'),
chi_min=fixed(-180),chi_max=fixed(180),name=fixed('EEI Diagnostic'),
targetlog=widgets.RadioButtons(options=ggg,description='Target log'));
max corr. EEI-GR=0.67 @ X=121
(array([[7428.42630264, 7435.45872116, 7425.73291988, ..., 7365.1181096 , 7404.79823094, 7428.42630264], [7495.42225411, 7482.5136826 , 7452.02939007, ..., 7467.71951901, 7490.45190418, 7495.42225411], [7513.9552278 , 7498.2078052 , 7464.60830614, ..., 7490.83755298, 7511.47895923, 7513.9552278 ], ..., [6023.46950442, 6067.24320665, 6117.96311962, ..., 5958.92050795, 5987.2091003 , 6023.46950442], [5950.60131951, 6006.29358257, 6070.08845732, ..., 5866.15390373, 5903.70086319, 5950.60131951], [5810.51213916, 5870.32862874, 5940.41053779, ..., 5724.21442582, 5761.63167342, 5810.51213916]]), array([ 0.56141335, 0.47160573, 0.29526148, 0.00438255, -0.29279421, -0.47504501, -0.56447414, -0.60788964, -0.62981623, -0.64111897, -0.64678172, -0.64927174, -0.64990583, -0.64942924, -0.64827605, -0.64669174, -0.64479095, -0.64258012, -0.63995585, -0.63667775, -0.63230233, -0.62604143, -0.61645495, -0.6007446 , -0.5730102 , -0.51970996, -0.40879206, -0.18253254, 0.14655723, 0.40148682, 0.52972325, 0.58999803, 0.62066129, 0.63782733, 0.64825489, 0.65501674, 0.6596313 , 0.66289751, 0.66524566, 0.66689371, 0.66791836, 0.66828378, 0.66784433, 0.6663231 , 0.66325689, 0.65787936, 0.64887528, 0.63384245, 0.60804653, 0.56141335]))