From b0ca7959aba06c38f3a63d72c3804340c07bf4b1 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 16 Feb 2016 09:46:22 -0800 Subject: [PATCH] removed MT1d example from this pr --- .../Examples/MT_1D_analytic_nlayer_Earth.py | 443 ------------------ 1 file changed, 443 deletions(-) delete mode 100644 SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py deleted file mode 100644 index 0cccc5f1..00000000 --- a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py +++ /dev/null @@ -1,443 +0,0 @@ -from scipy.constants import epsilon_0, mu_0 -import matplotlib.pyplot as plt -import numpy as np -from ipywidgets import * -#from SimPEG.EM.Utils import k, omega - -""" -MT1D: n layered earth problem -***************************** - -Author: Thibaut Astic -Contact: thast@eos.ubc.ca -Date: January 2016 - -This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics). - -We start by looking at Maxwell's equations in the electric -field \\\(\\\mathbf{E}\\) and the magnetic flux -\\\(\\\mathbf{H}\\) to write the wave equations -\\(\\ \nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) & -\\(\\ \nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\) - -Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of -\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\) -\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\) - -With U and D the Up and Down components of the E-field. - -The iteration from one layer to another is ensure by: - -\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = - P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\) - -And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0) -and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability) - -The layer 0 is assumed to be the air layer. - -""" - -#Frequency conversion -omega = lambda f: 2.*np.pi*f - -#Evaluate k wavenumber -k = lambda mu,sig,eps,f: np.sqrt(mu*mu_0*eps*epsilon_0*(2.*np.pi*f)**2.-1.j*mu*mu_0*sig*omega(f)) - -#Define a frquency range for a survey -frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.) - -#Functions to create random physical Perties for a n-layered earth -thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]), - np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1) - ,decimals =1)) - -sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]), - np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1) - ,decimals=3)) - -mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]), - np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1) - ,decimals=1)) - -eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]), - np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1) - ,decimals=1)) - -#Evaluate Impedance Z of a layer -ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k - -#Complex Cole-Cole Conductivity - EM utils -PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c))) - - -#Converted thickness array into top of layer array -def top(thick): - topv= np.zeros(len(thick)+1) - - topv[0]=-thick[0] - - for i in range(1,len(topv),1): - topv[i] = topv[i-1] + thick[i-1] - - return topv - -#Propagation Matrix and theirs inverses - -#matrix T for transition of Up and Down components accross a layer -T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_') - -Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_') - -#transition of Up and Down components accross a layer -UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD - - -#matrix P relating Up and Down components with E and H fields -P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_') - -Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2. - - -#Time Variation of E and H -E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D) -H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U) - -#Plot the configuration of the problem -def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): - - topn = top(thick) - widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.) - - ax.set_ylim([z.min(),z.max()]) - ax.set_xlim([-widthg,widthg]) - - ax.set_ylabel("Depth (m)", fontsize=16.) - ax.yaxis.tick_right() - ax.yaxis.set_label_position("right") - - #define filling for the different layers - hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ] - - #Write the physical properties of air - ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)), - xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data', - xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data', - fontsize=14.) - - ax.annotate(("$\epsilon_r$= %1i")%(eps[0]), - xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data', - xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data', - fontsize=14.) - - ax.annotate(("$\mu_r$= %1i")%(mu[0]), - xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data', - xytext=(0, -np.abs(z.max())/3.), textcoords='data', - fontsize=14.) - - #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern - for i in range(1,len(topn)-1,1): - if topn[i] == topn[i+1]: - pass - else: - ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)), - xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data', - xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data', - fontsize=14.) - - ax.annotate(("$\epsilon_r$= %1i")%(eps[i]), - xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data', - xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data', - fontsize=14.) - - ax.annotate(("$\mu_r$= %1.2f")%(mu[i]), - xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data', - xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data', - fontsize=14.) - - ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black') - ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10]) - - #Write the physical properties of the n-th layer and fill it with pattern - ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black') - ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10]) - - ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)), - xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data', - xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data', - fontsize=14.) - - ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]), - xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data', - xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data', - fontsize=14.) - - ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]), - xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data', - xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data', - fontsize=14.) - - #plot Trees! - ax.annotate("", - xy=(widthg/2., -1.*z.max()/5.), xycoords='data', - xytext=(widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data', - xytext=(widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data', - xytext=(widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data', - xytext=(1.2*widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data', - xytext=(1.2*widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data', - xytext=(1.2*widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data', - xytext=(1.5*widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data', - xytext=(1.5*widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.) - ) - - ax.annotate("", - xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data', - xytext=(1.5*widthg/2., 0.), textcoords='data', - arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) - ) - - - ax.invert_yaxis() - - return ax - -#Propagate Up and Down component for a certain frequency & evaluate E and H field - -def Propagate(f,H,sig,chg,taux,c,mu,eps,n): - - sigcm = np.zeros_like(sig,dtype='complex_') - - for j in range(1,len(sig)): - sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f) - - K = k(mu,sigcm,eps,f) - Z = ImpZ(f,mu,K) - - EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_') - UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_') - - UD[1,-1] = 1. - - for i in range(-2,-(n+2),-1): - - UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1] - UD = UD/((np.abs(UD[0,:]+UD[1,:])).max()) - - for j in range(0,n+1): - EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j] - - return UD, EH, Z ,K - - -#Evaluate the apparent resistivity and phase for a frequency range -def appres(F,H,sig,chg,taux,c,mu,eps,n): - - Res = np.zeros_like(F) - Phase = np.zeros_like(F) - App_ImpZ= np.zeros_like(F,dtype='complex_') - - for i in range(0,len(F)): - - UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n) - - App_ImpZ[i] = EH[0,1]/EH[1,1] - - Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i])) - Phase[i] = np.angle(App_ImpZ[i], deg = True) - - return Res,Phase - -#Evaluate Up, Down components, E and H field, for a frequency range, -#a discretized depth range and a time range (use to calculate envelope) -def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample): - - topc = top(H) - - layer = np.zeros(len(zsample),dtype=np.int)-1 - - Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') - Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') - Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') - Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') - UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_') - - for i in range(0,n+1,1): - layer = layer+(zsample>=topc[i])*1 - - for j in range(0,len(F)): - - UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n) - - for p in range(0,len(zsample)): - - UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]]) - - for q in range(0,len(tsample)): - - Exzt[p,q] = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F) - Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F) - Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F) - Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F) - - return Exzt,Hyzt,Uz,Dz,UDaux,layer - - -#Function to Plot Apparent Resistivity and Phase -def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): - - Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n) - - fig,ax = plt.subplots(1,2,figsize=(16,10)) - - ax[0].scatter(Res,F,color='black') - ax[0].set_xscale('Log') - ax[0].set_yscale('Log') - ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)]) - ax[0].set_ylim([F.min(),F.max()]) - ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black") - ax[0].set_ylabel('Frequency (Hz)',fontsize=16.) - ax[0].grid(which='major') - - ax0 = ax[0].twiny() - - ax0.set_xlim([0.,90.]) - ax0.set_ylim([F.min(),F.max()]) - ax0.scatter(Phase,F,color='purple') - ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple") - - zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.) - - ax[0].tick_params(labelsize=16) - ax[1].tick_params(labelsize=16) - ax0.tick_params(labelsize=16) - - if PlotEnvelope: - - widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0) - fenvelope1n=np.ones(100)*fenvelope - ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black') - - tc=np.arange(0.,1./fenvelope,0.01/(fenvelope)) - Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc) - - ax1=ax[1].twiny() - - ax[1].tick_params(labelsize=16) - ax1.tick_params(labelsize=16) - - ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16) - - ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16) - - ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))), - np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))), - color='blue', alpha=0.1) - - ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))), - np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))), - color='red', alpha=0.1) - - ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc) - ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()]) - ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()]) - else: - print 'No envelop (if True, might be slow)' - ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc) - ax[1].get_xaxis().set_ticks([]) - - plt.show() - -#Interactive MT for Notebook -def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope): - - frangn=frange(-5,5,100.) - sig3= np.array([0.,0.001,0.1, 0.001]) - thick3 = np.array([120000.,50.,50.]) - eps3=np.array([1.,1.,1.,1]) - mu3=np.array([1.,1.,1.,1]) - chg3=np.array([0.,0.1,0.,0.2]) - chg3_0=np.array([0.,0.1,0.,0.]) - taux3=np.array([0.,0.1,0.,0.1]) - c3=np.array([1.,1.,1.,1.]) - - sig3[1]=sigl1 - sig3[1]=10.**sig3[1] - sig3[2]=sigl2 - sig3[2]=10.**sig3[2] - sig3[3]=sigl3 - sig3[3]=10.**sig3[3] - mu3[1]=mul1 - mu3[2]=mul2 - mu3[3]=mul3 - eps3[1]=epsl1 - eps3[2]=epsl2 - eps3[3]=epsl3 - thick3[1]=h1 - thick3[2]=h2 - - PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope) - - -def run(n,plotIt=True): - # something to make a plot - - F = frange(-5.,5.,20) - H = thick(50.,100.,n) - sign = sig(-5.,0.,n) - mun = mu(1.,2.,n) - epsn = eps(1.,9.,n) - chg = np.zeros_like(sign) - taux = np.zeros_like(sign) - c = np.zeros_like(sign) - - Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n) - - if plotIt: - - PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True) - - return Res, Phase - -if __name__ == '__main__': - run(3) - - - - -