In [1]:
%pylab inline
In [2]:
import pandas as pd
In [3]:
data = pd.read_excel("../data/Kimberlite-2015-07-17.xls")
In [4]:
data_active = data.loc[np.logical_and((data['Facies'] == 'XVK')|(data['Facies'] == 'PK')|(data['Facies'] == 'HK')|(data['Facies'] == 'VK'), data.notnull()['Rinf']==True)][["Facies", "Peregrine ID", "(Latitude)", "(Longitude)", "Depth (m)","Mag Susc [SI]","Resistivity [Ohm.m]","Geometric Factor [m]","Sat Geometric Dens [g/cc]","Chargeability [ms]","Rinf","Ro","Rh","Qh","Ph", "pRh", "pQh","Rm","Qm","Pm", "pRm", "pQm","Rl","Ql","Pl", "pRl", "pQl","Re","Qe","Pe-f","Pe-i"]]
In [5]:
data_active
Out[5]:
In [6]:
def CPEfun(Rx, Fx, px, freq):
out = np.zeros_like(freq, dtype=complex128)
out = 1./(1./Rx + Qx*(np.pi*2*freq*1j)**px)
return out
def CPEfunElec(Rx, Qx, pex, pix, freq):
out = np.zeros_like(freq, dtype=complex128)
out = 1./(1./Rx + (1j)**pix*Qx*(np.pi*2*freq)**pex)
return out
def CPEfunSeries(Rx, Qx, px, freq):
out = np.zeros_like(freq, dtype=complex128)
out = Rx + 1./(Qx*(np.pi*2*freq*1j)**px)
return out
f0peak = lambda R, Q, P: (R*Q)**(-1./P)/np.pi/2.
taupeak = lambda R, Q, P: (R*Q)**(1./P)
rhoinf = lambda rhom, rhol, rho0: 1./(1./rho0+1./rhom+1./rhol)
charg = lambda rhoinf, rho0: (rho0-rhoinf) / rhoinf
In [7]:
def TKCColeColeParallel(frequency, PID, data):
Rh, Qh, Ph = data[data["Peregrine ID"]==PID]['pRh'].values[0], data[data["Peregrine ID"]==PID]['pQh'].values[0], data[data["Peregrine ID"]==PID]['Ph'].values[0]
Rm, Qm, Pm = data[data["Peregrine ID"]==PID]['pRm'].values[0], data[data["Peregrine ID"]==PID]['pQm'].values[0], data[data["Peregrine ID"]==PID]['Pm'].values[0]
Rl, Ql, Pl = data[data["Peregrine ID"]==PID]['pRl'].values[0], data[data["Peregrine ID"]==PID]['pQl'].values[0], data[data["Peregrine ID"]==PID]['Pl'].values[0]
geom = data[data["Peregrine ID"]==PID]['Geometric Factor [m]'].values[0]
fpeakm = f0peak(Rm, Qm, Pm)
fpeakl = f0peak(Rl, Ql, Pl)
# geom = 1.
rhom = CPEfunSeries(Rm, Qm, Pm, frequency)*geom
rhol = CPEfunSeries(Rl, Ql, Pl, frequency)*geom
rho0 = data[data["Peregrine ID"]==PID]['Ro'].values[0]*geom
rho = 1./(1./rho0+1./rhol)
m = (rho.real[0]-rho.real[-1])/rho.real[0]
rhoinf = rho0*(1.-m)
fig, ax = plt.subplots(1, 2, figsize = (15, 3))
ax[0].semilogx(frequency, rho.real, 'k-', lw=2)
ax1 = ax[0].twinx()
ax1.semilogx(frequency, (rho.imag), 'k--', lw=2)
ax1.invert_yaxis()
ax[0].grid(True)
ax[1].plot(rho.real, rho.imag, 'k-')
ax[1].invert_yaxis()
ax[1].grid(True)
print data[data["Peregrine ID"]==PID]['Facies'].values[0], PID
print "R0 = ", rho0
print "Rinf = ", rhoinf
print "Chargeability = ", m
print "Taum = ", 1./fpeakm
print "Taul = ", 1./fpeakl
print Pl, Ql, Rl, geom
return
In [8]:
data_active[data_active['Facies']=='PK']
Out[8]:
In [9]:
frequency = np.logspace(-4, 8, 211)
TKCColeColeParallel(frequency, "K1P-0825", data_active)
In [10]:
frequency = np.logspace(-4, 8, 211)
TKCColeColeParallel(frequency, "K1P-0807", data_active)
In [11]:
data_active[data_active['Facies']=='XVK']
Out[11]:
In [12]:
TKCColeColeParallel(frequency, "K2P-0031", data_active)
In [13]:
data_active[data_active['Facies']=='HK']
Out[13]:
In [14]:
TKCColeColeParallel(frequency, "K1P-0591", data_active)
In [15]:
data_active[data_active['Facies']=='VK']
Out[15]:
In [16]:
TKCColeColeParallel(frequency, "K1P-0589", data_active)
In [17]:
RlPK = data_active[data_active['Facies']=='PK']["pRl"].values[:]
QlPK = data_active[data_active['Facies']=='PK']["pQl"].values[:]
PlPK = data_active[data_active['Facies']=='PK']["Pl"].values[:]
fpeakPK = f0peak(RlPK, QlPK, PlPK)
rhoinfPK = rhoinf(data_active[data_active['Facies']=='PK']["Ro"].values[:], data_active[data_active['Facies']=='PK']["pRm"].values[:],data_active[data_active['Facies']=='PK']["pRl"].values[:])
mPK = charg(rhoinfPK, data_active[data_active['Facies']=='PK']["Ro"].values[:])
In [18]:
RlXVK = data_active[data_active['Facies']=='XVK']["pRl"].values[:]
QlXVK = data_active[data_active['Facies']=='XVK']["pQl"].values[:]
PlXVK = data_active[data_active['Facies']=='XVK']["Pl"].values[:]
fpeakXVK = f0peak(RlXVK, QlXVK, PlXVK)
rhoinfXVK = rhoinf(data_active[data_active['Facies']=='XVK']["Ro"].values[:], data_active[data_active['Facies']=='XVK']["pRm"].values[:],data_active[data_active['Facies']=='XVK']["pRl"].values[:])
mXVK = charg(rhoinfXVK, data_active[data_active['Facies']=='XVK']["Ro"].values[:])
In [19]:
RlVK = data_active[data_active['Facies']=='VK']["pRl"].values[:]
QlVK = data_active[data_active['Facies']=='VK']["pQl"].values[:]
PlVK = data_active[data_active['Facies']=='VK']["Pl"].values[:]
fpeakVK = f0peak(RlVK, QlVK, PlVK)
rhoinfVK = rhoinf(data_active[data_active['Facies']=='VK']["Ro"].values[:], data_active[data_active['Facies']=='VK']["pRm"].values[:],data_active[data_active['Facies']=='VK']["pRl"].values[:])
mVK = charg(rhoinfVK, data_active[data_active['Facies']=='VK']["Ro"].values[:])
In [20]:
RlHK = data_active[data_active['Facies']=='HK']["pRl"].values[:]
QlHK = data_active[data_active['Facies']=='HK']["pQl"].values[:]
PlHK = data_active[data_active['Facies']=='HK']["Pl"].values[:]
fpeakHK = f0peak(RlHK, QlHK, PlHK)
rhoinfHK = rhoinf(data_active[data_active['Facies']=='HK']["Ro"].values[:], data_active[data_active['Facies']=='HK']["pRm"].values[:],data_active[data_active['Facies']=='HK']["pRl"].values[:])
mHK = charg(rhoinfHK, data_active[data_active['Facies']=='HK']["Ro"].values[:])
In [21]:
import matplotlib as mpl
mpl.rcParams["font.size"] = 16
mpl.rcParams["text.usetex"] = True
In [32]:
figsize(5,5)
indactHK = 1./fpeakHK > 0.1
plt.semilogx(1./fpeakXVK*1e6, mXVK, 'kx', ms = 15, lw=3)
plt.semilogx(1./fpeakVK*1e6, mVK, 'bx', ms = 15)
plt.semilogx(1./fpeakHK[~indactHK]*1e6, mHK[~indactHK], 'mx', ms = 15)
plt.semilogx(1./fpeakPK*1e6, mPK, 'rx', ms = 15)
ylim(-0.1, 1.)
xlim(1e-5*1e6, 1e-2*1e6)
plt.grid(True)
plt.ylabel("Chargeability ")
plt.xlabel("Time constant (micro-s)")
plt.legend(("XVK", "VK", "HK", "PK"), bbox_to_anchor=(1.5, 1))
Out[32]:
In [23]:
mus = 1e6
plt.semilogx(1./fpeakXVK*mus, PlXVK,'ko')
plt.semilogx(1./fpeakVK*mus, PlVK, 'bo')
plt.semilogx(1./fpeakHK[~indactHK]*mus, PlHK[~indactHK], 'bo')
plt.semilogx(1./fpeakPK*mus, PlPK, 'mo')
ylim(0., 1.)
xlim(1e-5*mus, 1e-2*mus)
plt.grid(True)
# plt.legend(("XVK", "VK", "HK", "PK"), loc=4)
plt.ylabel("Frequency dependency ")
plt.xlabel("Time constant (micro-s)")
figsize(5,5)
In [24]:
plt.semilogy(data_active[data_active["Facies"]=="XVK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="XVK"]["Resistivity [Ohm.m]"], 'ko')
plt.semilogy(data_active[data_active["Facies"]=="VK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="VK"]["Resistivity [Ohm.m]"], 'bo')
plt.semilogy(data_active[data_active["Facies"]=="HK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="HK"]["Resistivity [Ohm.m]"], 'mo')
plt.semilogy(data_active[data_active["Facies"]=="PK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="PK"]["Resistivity [Ohm.m]"], 'ro')
plt.grid(True)
plt.legend(("XVK", "VK", "HK", "PK"), loc=4)
plt.xlabel("Susceptibility (SI)")
plt.ylabel("Resistivity (ohm-m)")
figsize(5,5)
In [25]:
plt.plot(data_active[data_active["Facies"]=="XVK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="XVK"]["Sat Geometric Dens [g/cc]"], 'ko')
plt.plot(data_active[data_active["Facies"]=="VK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="VK"]["Sat Geometric Dens [g/cc]"], 'bo')
plt.plot(data_active[data_active["Facies"]=="HK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="HK"]["Sat Geometric Dens [g/cc]"], 'ro')
plt.plot(data_active[data_active["Facies"]=="PK"]["Mag Susc [SI]"], data_active[data_active["Facies"]=="PK"]["Sat Geometric Dens [g/cc]"], 'go')
plt.grid(True)
plt.legend(("XVK", "VK", "HK", "PK"), loc=4)
plt.xlabel("Susceptibility (SI)")
plt.ylabel("Density (g/cc)")
figsize(5,5)
In [ ]:
In [ ]:
In [ ]:
In [ ]: