<jupyter_start><jupyter_text>Understand Physics of Induced Polarization Effects Airborne Electromagnetic Data**Author**: [<NAME>](https://github.com/sgkang)This notebooks illustrates how we can understand physics behind induced polarization (IP) effects in airborne electromagnetic (AEM) data. Here, we carry out two simulations which includes IP effects or not. Then we show measured voltages at the receiver location, and further explore computed electric fields and db/dt fields. This notebook was used to produce Figures 2, 3, 4 in Kang et al. (2019).<jupyter_code>from SimPEG import *
from SimPEG import EM
from scipy.constants import mu_0
import numpy as np
import scipy.sparse as sp
from simpegEMIP.StretchedExponential import SEInvImpulseProblem, SESurvey
import matplotlib.pyplot as plt
from pymatsolver import PardisoSolver
import matplotlib.patches as mpatches
from matplotlib.collections import PatchCollection
%matplotlib inline
import matplotlib 
figdir = "./figures/"<jupyter_output><empty_output><jupyter_text>Step1: Generate 2D Cylindrical Mesh<jupyter_code>cs, ncx, ncz, npad = 10., 40, 40, 25
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')<jupyter_output><empty_output><jupyter_text>Step2: Generate conductivity and chargeability model<jupyter_code>sig_1 = 0.001
sig_2 = 0.1
sigmaInf = np.ones(mesh.nC) * sig_1
airind = mesh.gridCC[:,2]>0.
z = 0.
radius = 200.
z_top = -50
z_bottom = -150
blkind  = np.logical_and(mesh.gridCC[:,0]<radius, (mesh.gridCC[:,2]<z_top)&(mesh.gridCC[:,2]>z_bottom))
sigmaInf[airind] = 1e-8
sigmaInf[blkind] = sig_2
eta_se, tau_se, c_se = 0.1, 0.001, 0.7
eta = np.zeros(mesh.nC)
eta[blkind] = eta_se
tau = np.ones(mesh.nC) * tau_se
c = np.ones(mesh.nC) * c_se
def rect2D(p1, p2):    
    xy = np.c_[np.r_[p1[0], p2[0], p2[0], p1[0], p1[0]], np.r_[p1[1], p1[1], p2[1], p2[1], p1[1]]]
    return xy
matplotlib.rcParams["font.size"] = 14
fig, ax = plt.subplots(1,1)
mesh.plotImage(eta*np.nan, mirror=True, ax=ax)
z=0
p1 = np.r_[-200, -150-z]
p2 = np.r_[200, -50-z]
xy_box = rect2D(p1, p2)
xy_surface = np.c_[np.r_[-400, 400], np.r_[0, 0]]
ax.plot(xy_box[:,0], xy_box[:,1], 'k-')
ax.plot(xy_surface[:,0], xy_surface[:,1], 'k-')
ax.set_aspect("equal")
ax.text(-300., 50., "Air")
ax.text(-300., -250., "Earth")
# ax.text(-100., -250, ("$\sigma_{half}$ = %.0e S/m")%(sig_1))
ax.text(-100., -250, ("$\sigma_{half}$ = ? $S/m$"))
# ax.text(-150, -120-z, "Chargeable")
ax.text(-300., 300, "SE parameters:")
ax.text(-100., 200+20, ("$\sigma_{\infty}$ = %.0e $S/m$")%(sig_2))
ax.text(-100., 150+20, ("$\eta$ = %.1f")%(eta_se))
ax.text(-100., 100+20, ("$\\tau$ = %.0f $ms$")%(tau_se*1e3))
# ax.text(-100., 200-150-z, ("$c$ = %.1f")%())
ax.set_xlim(-400, 400)
ax.set_ylim(-400, 400)
ax.plot(np.r_[-13., 13.], np.r_[30., 30.], 'g-', lw=4)
ax.text(40., 25., "Loop source", color='g')
ax.set_xlabel('$x \ (m)$')
ax.set_ylabel('$z \ (m)$')
plt.tight_layout()
fig.savefig(figdir+"model_raw", dpi=300)<jupyter_output><empty_output><jupyter_text>Step3: Generate `Problem` and `Survey`, and run simulation<jupyter_code>from simpegEMIP.TDEM import Problem3D_e, Problem3DEM_e
from simpegEMIP.TDEM import Rx, Survey
rxloc = np.array([[0., 0., 30.]])
srcloc = np.array([[0., 0., 30.]])

rx_emip = Rx.Point_dbdt(rxloc, np.logspace(np.log10(2e-5), np.log10(0.5), 31), 'z')
src_emip = EM.TDEM.Src.CircularLoop([rx_emip], waveform=EM.TDEM.Src.StepOffWaveform(), loc=srcloc, radius = 13.)
survey_emip = Survey([src_emip])
prb_emip = Problem3D_e(mesh, sigmaInf=sigmaInf, eta=eta, tau=tau, c=c)
prb_emip.timeSteps = [(1e-06, 5), (2.5e-06, 5), (5e-06, 5), (1e-05, 5), (2e-05, 5), (4e-05, 5), (8e-05, 5), (1.6e-04, 5), (4e-04, 5), (8e-04, 5), (1e-03, 5), (2e-03, 5), (4e-03, 5), (8e-03, 5), (1e-02, 5), (2e-02, 5), (4e-02, 5), (8e-02, 5), (1e-01, 5)]
prb_emip.Solver = PardisoSolver
prb_emip.pair(survey_emip)
m = np.r_[sigmaInf, eta, tau, c]
F = prb_emip.fields(m)
data = survey_emip.dpred(m, f=F)
rx_em = Rx.Point_dbdt(rxloc, np.logspace(np.log10(2e-5), np.log10(0.5), 31), 'z')
src_em = EM.TDEM.Src.CircularLoop([rx_em], waveform=EM.TDEM.Src.StepOffWaveform(), loc=srcloc, radius = 13.)
survey_em = Survey([src_em])
prb_em = Problem3DEM_e(mesh, sigma=sigmaInf)
prb_em.timeSteps = [(1e-06, 5), (2.5e-06, 5), (5e-06, 5), (1e-05, 5), (2e-05, 5), (4e-05, 5), (8e-05, 5), (1.6e-04, 5), (4e-04, 5), (8e-04, 5), (1e-03, 5), (2e-03, 5), (4e-03, 5), (8e-03, 5), (1e-02, 5), (2e-02, 5), (4e-02, 5), (8e-02, 5), (1e-01, 5)]
prb_em.Solver = PardisoSolver
prb_em.pair(survey_em)
m = np.r_[sigmaInf, eta, tau, c]
F_em = prb_em.fields(m)
data_em = survey_em.dpred(m, f=F_em)
data_ip = data - data_em<jupyter_output><empty_output><jupyter_text>Step4: Plot decay curves<jupyter_code>from simpegEMIP.TDEM.Utils import get_dobs_pn
area = np.pi*(13**2)
noise_floor = 1e-4
moment = np.pi*src_emip.radius**2 / 1e12
rectpts = rect2D([0.2, 1e-14/moment],[3., 1e-6/moment])
fig, axs = plt.subplots(1,2, figsize = (12, 4))
ms = 1e3
data_p, data_n = get_dobs_pn(data.flatten()/moment)
data_ip_p, data_ip_n = get_dobs_pn(data_ip.flatten()/moment)
axs[0].loglog(rx_emip.times*ms, data_n, 'k-', lw=2)
axs[0].loglog(rx_emip.times*ms, data_p, 'k--', lw=2)
axs[0].loglog(rx_emip.times*ms, -data_em.flatten()/moment, 'b-')
axs[0].loglog(rx_emip.times*ms, data_ip_n, 'r-')
axs[0].loglog(rx_emip.times*ms, data_ip_p, 'r--')
axs[0].legend(("$d^{obs}$ (+)", "$d^{obs}$ (-)", "$d^{\ F}$ (+)","$d^{IP}$ (+)", "$d^{IP}$ (-)"), fontsize = 12, loc=1)
for ax in axs:
    ax.set_xlabel("Time ($ms$)")    
    ax.set_xlim(rx_emip.times.min()*ms, rx_emip.times.max()*ms)
    ax.grid(True)
#     ax.set_xlim(2e-5*ms, 50*ms)
axs[0].set_ylim(1e-8, 1e3)
axs[1].set_ylim(1e-2, 1e2)
ax.grid(True)
# for t in tc:
tmin = rx_emip.times.min()*ms
tmax = rx_emip.times.max()*ms
axs[0].fill_between([tmin, tmax], np.r_[noise_floor, noise_floor], y2=np.r_[1e-8, 1e-8], color='grey', alpha=0.5)
axs[1].fill_between([tmin, tmax], np.r_[1e0, 1e0], y2=np.r_[1e-2, 1e-2], color='grey', alpha=0.5)
axs[0].text(0.1, 3e-7, "Below noise level")
axs[1].text(0.1, 1e-1, "$R<1$")
axs[0].set_title("(a) $d$, $d^{\ F}$ and $d^{IP}$")
axs[0].set_ylabel("Voltage ($pV/A$-$m^2$)")
axs[1].set_ylabel('$|d^{IP}|/|d^{F}|$ ratio ($R$)')
axs[1].set_title("(b) Ratio")
axs[1].loglog(rx_emip.times*ms, abs(data_ip.flatten()/data_em.flatten()), 'k')
    
plt.tight_layout()
fig.savefig(figdir+"data", dpi=300)<jupyter_output><empty_output><jupyter_text>Step5: Plot electric fields and db/dt fields<jupyter_code>x = np.arange(50)*10. + 10.
z = -np.arange(50)*10.
xz = Utils.ndgrid(x, np.r_[0.], z)
x_mirror = np.r_[-x[::-1], x]
xz_mirror = Utils.ndgrid(x_mirror, np.r_[0.], z)
X, Z = np.meshgrid(x_mirror, z)
Y = np.zeros_like(X)
P = mesh.getInterpolationMat(xz, 'Ey')
Pfx = mesh.getInterpolationMat(xz, 'Fx')
Pfz = mesh.getInterpolationMat(xz, 'Fz')
e_em = P * F_em[:,0,:]
e_emip = P * F[:,0,:]
dbdtx_em = Pfx * (-mesh.edgeCurl*F_em[:,0,:])
dbdtz_em = Pfz * (-mesh.edgeCurl*F_em[:,0,:])
dbdtx_emip = Pfx * (-mesh.edgeCurl*F[:,0,:])
dbdtz_emip = Pfz * (-mesh.edgeCurl*F[:,0,:])
def get_e_dbdt_em(i_time):
    temp_r = dbdtx_em[:,i_time].reshape((50, 50), order="F")
    temp_l = dbdtx_em[:,i_time].reshape((50, 50), order="F")[::-1, :]
    temp_x = np.vstack((-temp_l, temp_r))
    temp_r = dbdtz_em[:,i_time].reshape((50, 50), order="F")
    temp_l = dbdtz_em[:,i_time].reshape((50, 50), order="F")[::-1, :]
    temp_z = np.vstack((temp_l, temp_r))
    temp_abs = np.sqrt(temp_x**2+temp_z**2)
    temp_r = e_em[:,i_time].reshape((50, 50), order="F")
    temp_l = e_em[:,i_time].reshape((50, 50), order="F")[::-1, :]
    temp = np.vstack((-temp_l, temp_r))
    return temp.T, temp_abs.T, temp_x.T, temp_z.T
def get_e_dbdt_emip(i_time):
    temp_r = dbdtx_emip[:,i_time].reshape((50, 50), order="F")
    temp_l = dbdtx_emip[:,i_time].reshape((50, 50), order="F")[::-1, :]
    temp_x = np.vstack((-temp_l, temp_r))
    temp_r = dbdtz_emip[:,i_time].reshape((50, 50), order="F")
    temp_l = dbdtz_emip[:,i_time].reshape((50, 50), order="F")[::-1, :]
    temp_z = np.vstack((temp_l, temp_r))
    temp_abs = np.sqrt(temp_x**2+temp_z**2)
    temp_r = e_emip[:,i_time].reshape((50, 50), order="F")
    temp_l = e_emip[:,i_time].reshape((50, 50), order="F")[::-1, :]
    temp = np.vstack((-temp_l, temp_r))
    return temp.T, temp_abs.T, temp_x.T, temp_z.T

def get_nksip(nkip, M):
    return M[::nksip, ::nskip]
p1 = np.r_[-radius, z_bottom]
p2 = np.r_[radius, z_top]
xy_box = rect2D(p1, p2)
times_for_plot = np.r_[0.01, 0.3, 6, 50] * 1e-3
i_times = [np.argmin(abs(prb_em.times-time)) for time in times_for_plot]
matplotlib.rcParams["font.size"] = 12
matplotlib.rcParams["image.cmap"] = "viridis"

fig, axs = plt.subplots(4,2, figsize = (10, 8))
nskip = 5

for i in range(4):
    temp, temp_abs, temp_x, temp_z = get_e_dbdt_emip(i_times[i])
    im1 = axs[i,0].pcolormesh(x_mirror, z, temp)
    temp_abs_skip = temp_abs[::nskip, ::nskip]    
    im2 = axs[i,1].pcolormesh(x_mirror, z, temp_abs/moment)
    vmax = temp_abs.max()
    x_skip, z_skip, temp_x_skip, temp_z_skip = x_mirror[::nskip], z[::nskip], temp_x[::nskip, ::nskip]/vmax, temp_z[::nskip, ::nskip]/vmax,
    inds_p = temp_z_skip>0.
    temp_x_skip_p = temp_x_skip.copy()
    temp_z_skip_p = temp_z_skip.copy()
    temp_x_skip_n = temp_x_skip.copy()
    temp_z_skip_n = temp_z_skip.copy()
    temp_x_skip_p[~inds_p] = np.nan
    temp_z_skip_p[~inds_p] = np.nan
    temp_x_skip_n[inds_p] = np.nan
    temp_z_skip_n[inds_p] = np.nan    
    
    im2_1 = axs[i,1].quiver(
        x_skip, z_skip, temp_x_skip_p, temp_z_skip_p,
        scale=20, color='r', alpha=0.8,
    )
    im2_2 = axs[i,1].quiver(
        x_skip, z_skip, temp_x_skip_n, temp_z_skip_n,
        scale=20, color='w', alpha=0.8,
    )
    axs[i,0].plot(xy_box[:,0], xy_box[:,1], 'w-')
    axs[i,1].plot(xy_box[:,0], xy_box[:,1], 'w-')
    
    axs[i,0].set_aspect(1)
    axs[i,1].set_aspect(1)
    vmin, vmax = im1.get_clim()
    cb1 = plt.colorbar(im1, ax=axs[i,0], format="%.e", fraction=0.02, ticks=np.linspace(vmin, vmax,3))
    vmin, vmax = im2.get_clim()    
    cb1 = plt.colorbar(im2, ax=axs[i,1], format="%.e", fraction=0.02, ticks=np.linspace(vmin, vmax,3))
    axs[i,1].set_yticks([])
    axs[i,0].text(-400, -450, ("Time at %.2f $ms$")%((prb_em.times[i_times[i]])*1e3), fontsize=12)    
    if i<3:
        axs[i,0].set_xticks([])        
        axs[i,1].set_xticks([]) 
    if i==3:
        axs[i, 0].set_xlabel("$x$ ($m$)")        
        axs[i, 1].set_xlabel("$x$ ($m$)")                
    axs[i, 0].set_ylabel("$z$ ($m$)")
axs[0,0].set_title("(a) $e_y$ ($V/m$)", fontsize=16)
axs[0,1].set_title("(b) $d\\vec{b}/dt$ ($pV/A$-$m^2$)", fontsize=16)

axs[1,0].plot(150, -200, 'o', markersize = 10, markerfacecolor='w', markeredgecolor='k', markeredgewidth=1)
axs[1,0].plot(150, -200, 'kx', markersize = 7, markeredgewidth=2)
axs[1,0].plot(-150, -200, 'o', markersize = 10, markerfacecolor='w', markeredgecolor='k', markeredgewidth=1)
axs[1,0].plot(-150, -200, 'ko', markersize = 5)

axs[2,0].plot(-150, -200, 'o', markersize = 10, markerfacecolor='w', markeredgecolor='k', markeredgewidth=1)
axs[2,0].plot(-150, -200, 'kx', markersize = 7, markeredgewidth=2)
axs[2,0].plot(150, -200, 'o', markersize = 10, markerfacecolor='w', markeredgecolor='k', markeredgewidth=1)
axs[2,0].plot(150, -200, 'ko', markersize = 5)
plt.tight_layout()
fig.savefig(figdir + "fields", dpi=200)<jupyter_output><empty_output><jupyter_text>Step6: Generate movies for EM wave propagation<jupyter_code>## if you want to generate movies set run=True
run = False
if run:
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.animation as animation
    video_names = ['em.mp4', 'emip.mp4']
    for i in range(2):
        fig = plt.figure(figsize=(10, 10))
        temp, temp_abs, temp_x, temp_z = get_e_dbdt_em(10)
        nskip = 3
        ax1 = plt.subplot(211)
        ax2 = plt.subplot(212)
        ax1.plot(xy_box[:,0], xy_box[:,1], 'w-')
        ax2.plot(xy_box[:,0], xy_box[:,1], 'w-')

        ax2.set_aspect(1)
        time_text = ax1.text(-400, -400, '', fontsize=16)
        def init():
            temp, temp_abs, temp_x, temp_z = get_e_dbdt_em(0)
            im1 = ax1.pcolormesh(x_mirror, z, temp)
            temp_abs_skip = temp_abs[::nskip, ::nskip]
            im2 = ax2.pcolormesh(x_mirror, z, temp_abs)
            vmax = temp_abs.max()
            im2_1 = ax2.quiver(
                x_mirror[::nskip], z[::nskip], 
                temp_x[::nskip, ::nskip]/temp_abs_skip, temp_z[::nskip, ::nskip]/temp_abs_skip, 
                scale=30, color='w', alpha=0.8,
            )
            time_text.set_text('')
            return im1, im2, im2_1, time_text

        def animate(i_time):    
            if i == 0:
                temp, temp_abs, temp_x, temp_z = get_e_dbdt_em(i_time)
            else:
                temp, temp_abs, temp_x, temp_z = get_e_dbdt_emip(i_time)
            im1 = ax1.pcolormesh(x_mirror, z, temp)
            temp_abs_skip = temp_abs[::nskip, ::nskip]    
            im2 = ax2.pcolormesh(x_mirror, z, temp_abs)
            vmax = temp_abs.max()
            x_skip, z_skip, temp_x_skip, temp_z_skip = x_mirror[::nskip], z[::nskip], temp_x[::nskip, ::nskip]/vmax, temp_z[::nskip, ::nskip]/vmax,
            inds_p = temp_z_skip>0.
            temp_x_skip_p = temp_x_skip.copy()
            temp_z_skip_p = temp_z_skip.copy()
            temp_x_skip_n = temp_x_skip.copy()
            temp_z_skip_n = temp_z_skip.copy()
            temp_x_skip_p[~inds_p] = np.nan
            temp_z_skip_p[~inds_p] = np.nan
            temp_x_skip_n[inds_p] = np.nan
            temp_z_skip_n[inds_p] = np.nan    

            im2_1 = ax2.quiver(
                x_skip, z_skip, temp_x_skip_p, temp_z_skip_p,
                scale=30, color='r', alpha=0.8,
            )
            im2_2 = ax2.quiver(
                x_skip, z_skip, temp_x_skip_n, temp_z_skip_n,
                scale=30, color='w', alpha=0.8,
            )

            time_text.set_text(("Time at %.1e ms")%((prb_em.times[i_time])*1e3)) 
            ax1.set_aspect(1)
            ax2.set_aspect(1)
            return im1, im2, im2_1, im2_2, time_text

        ani = animation.FuncAnimation(fig, animate, prb_em.times.size-1, interval=200, blit=True,
                                        repeat_delay=1000)
        ani.save(video_names[i], dpi=100)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2049129137464437, 'fraction_numerical': 0.0599542016515162, 'mean_word_length': 3.78963110667996, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '50685', 'license': 'MIT', 'n_tokens_mistral': 7037, 'n_tokens_neox': 6599, 'n_words': 1139, 'path': 'notebooks/2-Understand-Physics.ipynb', 'repo_name': 'simpeg-research/kang-2018-AEM', 'size': '14415'}