<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Where are the diamonds? Using Northern Lights
Step1:  Scipy 2016 Poster
Step2:  Details of the physics at
Step3:  Paraview view
<img src="Figures/NSEMorignalModelPV.png" width=500>
Load stations locations and frequency range
Below with load in the information needed for solving the problem
Step4:  Types of data
The natural source of the MT method is unknown and cannot be seporated from the induced signal due to the Earth's conductivity or explicitly modelled numercially. To address this the ratio's of the meaured fields are used, and assuming that the source is the same, the source cancels out, leaving only the secondary signal in the data. 
Impedance
Is the ratio of the horizontal electrical over magnetic fields. It is a complex tensor, defined as
$$\begin{align}
\textbf{Z}(\omega) =
\begin{bmatrix}
Z_{xx}(\omega) & Z_{xy}(\omega)\
Z_{yx}(\omega) & Z_{yy}(\omega)
\end{bmatrix} =
\begin{bmatrix}
\frac{E_x(\omega)}{H_x(\omega)} & \frac{E_x(\omega)}{H_y(\omega)}\
\frac{E_y(\omega)}{H_x(\omega)} & \frac{E_y(\omega)}{H_y(\omega)}
\end{bmatrix} \label{eq
Step5:  Run the inversions on a cluster
To image the Earth, based on the data we forward modelled above, we use inversion. As with the forward modelling, the inversion is a computationally intensive procedure, and as before we ran yours on a cluster. A single iteration takes at the order of 4 hours to solve using ~ 20Gb of ram. 
We present 2 inversion results, for the impedance and the tipper. 
Links to the files used to run the inversions.
Off-diagonal impedance inversion
Tipper inversion
Load Inversion results
<ASSISTANT_TASK:>
Python Code:
import SimPEG as simpeg
from SimPEG import NSEM
import MT_poster_utils 
from IPython.display import HTML
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Explanation: Where are the diamonds? Using Northern Lights
End of explanation
HTML("Figures/Magnetotelluric_Movie_ThibautAstic.html")
Explanation: Scipy 2016 Poster:
https://drive.google.com/open?id=0B4oB1gnNZnWNcmdncHh3WVZqa1U
What is Magnetotelluric?
Short answer: The magnetotelluric method is a geophysical technique that relays on Natrual Eletromagnetic (EM) waves as a source of energy to image the subsurface. The method is widely used in geophysical exploration, specially imaging deep structures such as in geothermal, kimberlite or techtonic exploration. 
Long answer: Keep on reading...
What is our goal in this notebook?
In this notebook with are sharing the work we did for our SciPy 2016 poster, where we used NSEM code that are built with the SimPEG package (http://simpeg.xyz) to make data from a geological model (forward model) and to image the model again (invert the data). 
This poster is a part of a community effort, using a diamond exploration example as motivation to use the large spectrum of geophysical methods within SimPEG for the same geological setting. We want to demonstrate the advantages of having a shared language between the governing physics of all these problems, the written codes, and other utilitices used in the work. We will show
Explain the background theory's
Show and share the code used 
Visulize the data and models
Conceptual idea
The northern lights are generated by interaction of charged particles and the Earth’s magnetic field, and along with lightning strikes around the world, generate electromagnetic energy that induces electrical currents in the ground. These currents are affected by electrical conductivity contrasts and by measuring electromagnetic fields at the surface, we can image the subsurface structures using geophysical numerical simulation and optimization.
<img src="Figures/MTDiamonds-02.png" width=500a>
Animation of the physical phenomenon
End of explanation
# Load the geological discretized model
mesh, modelDict = simpeg.Mesh.TensorMesh.readVTK('./datafiles/nsmesh_model.vtr')
sigma = modelDict['S/m']
# Print model information
print mesh.nC," cells"
print mesh
# Define the area of interest in UTM (meters)
bw, be = 557100, 557580
bs, bn = 7133340, 7133960
bb, bt = 0,480
# View the model
slice=20
fig,ax=plt.subplots(figsize=(10,8))
modelplot = mesh.plotSlice(np.log10(sigma),normal='Y',ind=slice,ax=ax,grid=True, pcolorOpts={"cmap":"viridis"})
ax.set_xlim([bw,be])
ax.set_ylim([0,bt])
ax.set_aspect('equal')
plt.colorbar(modelplot[0])
ax.set_title("Discretized Model",fontsize=24)
Explanation: Details of the physics at:
http://em.geosci.xyz/content/maxwell3_fdem/MT/MT_N_layered_Earth.html
Physics (quasi-static approximation)
The Maxwell's equations govern the propagation of EM fields in the Earth. In frequency domain (and using the quasi-static approximation and a $e^{i\omega t}$ time relationship), they are given as 
$$
\begin{align}
\nabla \times \textbf{E} + i \omega \textbf{B} &= 0   \
\nabla \times \frac{1}{\mu} \textbf{B} - \sigma \textbf{E} &=   \textbf{s}_E 
\end{align}
$$
The system can be reduced to a single equation 
$$
\begin{align}
(\nabla \times \frac{1}{\mu} \nabla \times + i \omega \sigma) \textbf{E}_s &=   - i \omega \Delta \sigma \textbf{E}_p 
\end{align}
$$
which is then solved numerically using finite volume methods, in SimPEG.
Load the discrezited geological model
Original isosurfaces model
We start with a geological model which is based on a surfaces derived from diamond exploration. 
<img src="Figures/GeoModel.png" width=500>
We need to discritize the Earth and correlate the geological units with physical property (conductivity) need to generate the data. 
Below we use built in functions in SimPEG to read VTK (http://www.vtk.org) files of the physical property model in to mesh and model of the conductivity.
End of explanation
# Load stations locations and frequency range
locs = np.load('./datafiles/MTlocations.npy')
freqList = np.load('./datafiles/MTfrequencies.npy')
# View a scatter plot of the locations at the surface
plt.scatter(locs[:,0],locs[:,1])
# List of the frequencies used for the problem
print freqList
Explanation: Paraview view
<img src="Figures/NSEMorignalModelPV.png" width=500>
Load stations locations and frequency range
Below with load in the information needed for solving the problem
End of explanation
# Load the data - stored as numpy.recArray
mtData = np.load('./datafiles/MTdata.npy')
# Plot data
fig, axes, csList = MT_poster_utils.pseudoSect_OffDiagTip_RealImag(mtData,{'y':7133627.5},colBarMode='each')
# Make the plot
fig, axes, csList = MT_poster_utils.pseudoSect_FullImpTip_RealImag(mtData,{'y':7133627.5},colBarMode='each')
Explanation: Types of data
The natural source of the MT method is unknown and cannot be seporated from the induced signal due to the Earth's conductivity or explicitly modelled numercially. To address this the ratio's of the meaured fields are used, and assuming that the source is the same, the source cancels out, leaving only the secondary signal in the data. 
Impedance
Is the ratio of the horizontal electrical over magnetic fields. It is a complex tensor, defined as
$$\begin{align}
\textbf{Z}(\omega) =
\begin{bmatrix}
Z_{xx}(\omega) & Z_{xy}(\omega)\
Z_{yx}(\omega) & Z_{yy}(\omega)
\end{bmatrix} =
\begin{bmatrix}
\frac{E_x(\omega)}{H_x(\omega)} & \frac{E_x(\omega)}{H_y(\omega)}\
\frac{E_y(\omega)}{H_x(\omega)} & \frac{E_y(\omega)}{H_y(\omega)}
\end{bmatrix} \label{eq:ImpFull}
\end{align}$$
Tipper
Is the ratio of the vertical magnetic field over horizontal. It is a complex vector.
$$\begin{align}
\textbf{T}(\omega) =
\begin{bmatrix}
T_{zx}(\omega) \
T_{zy}(\omega)
\end{bmatrix} =
\begin{bmatrix}
\frac{H_z(\omega)}{H_x(\omega)} \
\frac{H_z(\omega)}{H_y(\omega)}
\end{bmatrix} \label{eq:TipFull}
\end{align}$$
Run Forward modeling script on a cluster
Running the numerical codes to solve the proble are resource heavy. The file below is the script we used to run the forward modelling on a cluster.
Forward Modeling
Data Visualization
The data are visulized as sections.
End of explanation
#Load Model Off-diagonal
mesh,inv= simpeg.Mesh.TensorMesh.readVTK('./datafiles/recoveredMod_off_it18.vtr')
siginvoff=inv['S/m']
#Load Model Tipper
mesh,invtip= simpeg.Mesh.TensorMesh.readVTK('./datafiles/recoveredMod_tip_it36.vtr')
siginvtip=invtip['S/m']
MT_poster_utils.CompareInversion(mesh,sigma,siginvoff,siginvtip,slice_ver=20,slice_hor=40)
Explanation: Run the inversions on a cluster
To image the Earth, based on the data we forward modelled above, we use inversion. As with the forward modelling, the inversion is a computationally intensive procedure, and as before we ran yours on a cluster. A single iteration takes at the order of 4 hours to solve using ~ 20Gb of ram. 
We present 2 inversion results, for the impedance and the tipper. 
Links to the files used to run the inversions.
Off-diagonal impedance inversion
Tipper inversion
Load Inversion results
End of explanation
<END_TASK>