In [134]:
import pandas as pd
import matplotlib
matplotlib.rcParams.update({'font.size': 12})
matplotlib.rcParams.update({'grid.color': 'white', 'grid.linewidth':1})
%pylab inline
In [141]:
path = "../data/HVC2014_10Grenon/"
fnameobs = "BC13867-A 2014-10-23.z"
fnamepred = "HVC 2014-10-27simplified.csv"
pathobs = path+fnameobs
pathpred = path+fnamepred
# #windows
# pathforPK = "..\\data\\nt01222a.z"
# pathfordata = "..\\data\\Kimberlite-2015-07-17.xls"
temp = np.loadtxt(pathobs, skiprows=11, delimiter=",")
data = pd.read_csv(pathpred)
data_active = data.loc[np.logical_and(data["Fl"]<100., data["Fm"]<1.)]
# data_active = data.loc[np.logical_and((data['Facies'] == 'XVK')|(data['Facies'] == 'PK')|(data['Facies'] == 'HK')|(data['Facies'] == 'VK'), data.notnull()['Rinf']==True)][["Facies", "0LabID (PCG)", "Peregrine ID", "(Latitude)", "(Longitude)", "Depth (m)","Mag Susc [SI]","Resistivity [Ohm.m]","Geometric Factor [m]","Sat Geometric Dens [g/cc]","Chargeability [ms]","Rinf","Linf","Ro","Rh","Qh","Ph", "Fh","pRh", "pQh","Rm","Qm","Pm", "Fm","pRm", "pQm","Rl","Ql","Pl", "Fl", "pRl", "pQl","Re","Qe","Pe-f","Pe-i"]]
In [142]:
# data_active[data['Facies'] == 'XVK']
In [143]:
data_active
Out[143]:
In [144]:
labid = fnameobs
In [145]:
def readfun(i):
Linf = data_active.iloc[i]['Linf']
Rinf = data_active.iloc[i]['Rinf']
Rh = data_active.iloc[i]['Rh']
Ph = data_active.iloc[i]['Ph']
Fh = data_active.iloc[i]['Fh']
Rm = data_active.iloc[i]['Rm']
Pm = data_active.iloc[i]['Pm']
Fm = data_active.iloc[i]['Fm']
Rl = data_active.iloc[i]['Rl']
Pl = data_active.iloc[i]['Pl']
Fl = data_active.iloc[i]['Fl']
Re = data_active.iloc[i]['Re']
Qe = data_active.iloc[i]['Qe']
Pef = data_active.iloc[i]['Pe-f']
Pei = data_active.iloc[i]['Pe-i']
fname = data_active.iloc[i]['file']
return Linf, Rinf, Rh, Fh, Ph, Rl, Fl, Pl, Rm, Fm, Pm, Re, Qe, Pef, Pei, fname
In [146]:
# Linf = data_active[data['file'] == labid]['Linf'].values[0]
# Rinf = data_active[data['file'] == labid]['Rinf'].values[0]
# Rh = data_active[data['file'] == labid]['Rh'].values[0]
# Ph = data_active[data['file'] == labid]['Ph'].values[0]
# Fh = data_active[data['file'] == labid]['Fh'].values[0]
# Rm = data_active[data['file'] == labid]['Rm'].values[0]
# Pm = data_active[data['file'] == labid]['Pm'].values[0]
# Fm = data_active[data['file'] == labid]['Fm'].values[0]
# Rl = data_active[data['file'] == labid]['Rl'].values[0]
# Pl = data_active[data['file'] == labid]['Pl'].values[0]
# Fl = data_active[data['file'] == labid]['Fl'].values[0]
# Re = data_active[data['file'] == labid]['Re'].values[0]
# Qe = data_active[data['file'] == labid]['Qe'].values[0]
# Pef = data_active[data['file'] == labid]['Pe-f'].values[0]
# Pei = data_active[data['file'] == labid]['Pe-i'].values[0]
In [ ]:
In [147]:
import sys
#windows
# sys.path.append("..\\codes\\")
sys.path.append("../codes/")
from Zarcfit import *
In [193]:
def filehandler(i, viewtype):
Linf, Rinf, Rh, Fh, Ph, Rl, Fl, Pl, Rm, Fm, Pm, Re, Qe, Pef, Pei,predname =readfun(i)
pathobs = path+predname
temp = np.loadtxt(pathobs, skiprows=11, delimiter=",")
# observed impedance
obs = temp[:,4]+1j*temp[:,5]
# frequency
frequency = temp[:,0].copy()
# Initiate Zarfit
zarc = Zarcfit(obs, frequency)
zarc.SetParametersSeries(Linf, Rinf, Rh, Fh, Ph, Rl, Fl, Pl, Rm, Fm, Pm, Re, Qe, Pef, Pei)
predS = zarc.Zseries(frequency)
predP = zarc.Zparallel(frequency)
Ro = zarc.R0
pRinf = 1./(1./Ro+1./zarc.pRl)
m = (Ro-pRinf) / Ro
tau = (Ro*zarc.pQl/m)**(1./zarc.Pl)
if viewtype =="ampphase":
plotAmpPhase(obs, predS, predP, frequency)
elif viewtype =="zarc":
plotZarcdlogZdlogF(obs, predS, predP, frequency)
print (predname)
print (Ro, pRinf, m, tau, Pl)
return
In [194]:
from IPython.html import widgets
In [195]:
ndata = len(data_active.index)
In [196]:
widgets.interact(filehandler, i=widgets.IntSlider(min=0, max=ndata, step=1, value=0),
viewtype=widgets.ToggleButtons(options=['ampphase','zarc']))
In [154]:
In [133]:
plotZarcdlogZdlogF(obs, predS, predP, frequency)
In [186]:
def plotAmpPhase(obs, predS, predP, frequency):
matplotlib.rcParams.update({'font.size': 12})
fig = plt.figure(figsize = (12, 5))
amp = np.sqrt(obs.real**2+obs.imag**2)
amp_predS = np.sqrt(predS.real**2+predS.imag**2)
amp_predP = np.sqrt(predP.real**2+predP.imag**2)
phase = np.angle(obs, deg=True)
phase_predS = np.angle(predS, deg=True)
phase_predP = np.angle(predP, deg=True)
ax = plt.subplot(121)
ax1 = plt.subplot(122)
ax.patch.set_facecolor('black')
ax1.patch.set_facecolor('black')
ax.plot(frequency, amp, 'lime', lw=3, marker=".", ms=8)
ax.plot(frequency, amp_predS, 'deepskyblue', lw=2)
ax.plot(frequency, amp_predP, 'b.')
ax.set_xscale('log')
ax.set_yscale('linear')
ax1.plot(frequency, -phase, 'lime', lw=3, marker=".", ms=8)
ax1.plot(frequency, -phase_predS, 'deepskyblue', lw=2)
ax1.plot(frequency, -phase_predP, 'b.')
ax1.set_xscale('log')
ax1.set_yscale('log')
# ax.set_ylim(amp.min(), amp.max())
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("|Z| (Ohm)")
ax1.set_xlabel("Frequency")
ax1.set_ylabel("Phase (degree)")
ax.grid(True)
ax.yaxis.grid(True, which='minor')
plt.tick_params(axis='y', which='minor')
ax1.grid(True)
ax1.yaxis.grid(True, which='minor')
plt.tick_params(axis='y', which='minor')
leg = ax1.legend(("Observed", "Predicted Series", "Predicted Parallel"), fontsize = 12, bbox_to_anchor = (1.6, 1))
plt.show()
return
In [ ]:
In [187]:
# fig = plt.figure(figsize = (11, 5))
# amp = np.sqrt(obs.real**2+obs.imag**2)
# amp_predS = np.sqrt(predS.real**2+predS.imag**2)
# amp_predP = np.sqrt(predP.real**2+predP.imag**2)
# ax = plt.subplot(121)
# ax1 = plt.subplot(122)
# ax.patch.set_facecolor('black')
# ax1.patch.set_facecolor('black')
# ax.plot(obs.real, -obs.imag, 'lime', lw=3, marker=".", ms=8)
# ax.plot(predS.real, -predS.imag, 'deepskyblue', lw=2)
# ax.plot(predS.real, -predP.imag, 'b.')
# ax1.semilogx(frequency[1:]*0.5+frequency[0:-1]*0.5, diff(np.log10(amp))/diff(np.log10(frequency)), 'lime', lw=3, marker=".", ms=8)
# ax1.semilogx(frequency[1:]*0.5+frequency[0:-1]*0.5, diff(np.log10(amp_predS))/diff(np.log10(frequency)), 'deepskyblue', lw=2)
# ax1.semilogx(frequency[1:]*0.5+frequency[0:-1]*0.5, diff(np.log10(amp_predP))/diff(np.log10(frequency)), 'b.', lw=2)
# ax.set_xlim(obs.real.min(), obs.real.max())
# ax.set_ylim((-obs.imag).min(), (-obs.imag).max())
# ax.set_xlabel("Real (Ohm)")
# ax.set_ylabel("- Imaginary (Ohm)")
# ax1.set_xlabel("Frequency")
# ax1.set_ylabel("dlog(Amp)/dlog(Freq)")
# ax.grid(True)
# ax.yaxis.grid(True, which='minor')
# plt.tick_params(axis='y', which='minor')
# ax1.grid(True)
# ax1.yaxis.grid(True, which='minor')
# plt.tick_params(axis='y', which='minor')
# leg = ax1.legend(("Observed", "Predicted Series", "Predicted Parallel"), fontsize = 12, bbox_to_anchor = (1.6, 1))
In [188]:
def plotZarcdlogZdlogF(obs, predS, predP, frequency):
fig = plt.figure(figsize = (12, 5))
amp = np.sqrt(obs.real**2+obs.imag**2)
amp_predS = np.sqrt(predS.real**2+predS.imag**2)
amp_predP = np.sqrt(predP.real**2+predP.imag**2)
ax = plt.subplot(121)
ax1 = plt.subplot(122)
ax.patch.set_facecolor('black')
ax1.patch.set_facecolor('black')
ax.plot(obs.real, -obs.imag, 'lime', lw=3, marker=".", ms=8)
ax.plot(predS.real, -predS.imag, 'deepskyblue', lw=2)
ax.plot(predS.real, -predP.imag, 'b.')
ax1.semilogx(frequency[1:]*0.5+frequency[0:-1]*0.5, diff(np.log10(amp))/diff(np.log10(frequency)), 'lime', lw=3, marker=".", ms=8)
ax1.semilogx(frequency[1:]*0.5+frequency[0:-1]*0.5, diff(np.log10(amp_predS))/diff(np.log10(frequency)), 'deepskyblue', lw=2)
ax1.semilogx(frequency[1:]*0.5+frequency[0:-1]*0.5, diff(np.log10(amp_predP))/diff(np.log10(frequency)), 'b.', lw=2)
ax.set_xlim(obs.real.min(), obs.real.max())
ax.set_ylim((-obs.imag).min(), (-obs.imag).max())
ax.set_xlabel("Real (Ohm)")
ax.set_ylabel("Frequency (Hz)")
ax1.set_xlabel("Frequency (Hz)")
ax1.set_ylabel("dlog(Amp)/dlog(Freq)")
ax.grid(True)
ax.yaxis.grid(True, which='minor')
plt.tick_params(axis='y', which='minor')
ax1.grid(True)
ax1.yaxis.grid(True, which='minor')
plt.tick_params(axis='y', which='minor')
leg = ax1.legend(("Observed", "Predicted Series", "Predicted Parallel"), fontsize = 12, bbox_to_anchor = (1.6, 1))
plt.show()
return
In [ ]: