import numpy as np
from geoscilabs.mag import Mag, Simulator
from SimPEG.potential_fields import magnetics as mag
from SimPEG import utils, data
from discretize import TensorMesh

"""
Explanation: This is the <a href="https://jupyter.org/">Jupyter Notebook</a>, an interactive coding and computation environment. For this lab, you do not have to write any code, you will only be running it. 
To use the notebook:
- "Shift + Enter" runs the code within the cell (so does the forward arrow button near the top of the document)
- You can alter variables and re-run cells
- If you want to start with a clean slate, restart the Kernel either by going to the top, clicking on Kernel: Restart, or by "esc + 00" (if you do this, you will need to re-run the following block of code before running any other cells in the notebook) 
This notebook uses code adapted from 
SimPEG
- Cockett, R., S. Kang, L.J. Heagy, A. Pidlisecky, D.W. Oldenburg (2015, in review), SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers and Geosciences
End of explanation
"""


#Input parameters
fileName = 'https://github.com/geoscixyz/geosci-labs/raw/main/assets/mag/data/DO27_TMI.dat'
xyzd = np.genfromtxt(fileName, skip_header=3)
B = np.r_[60308, 83.8, 25.4]
survey, dobj = Mag.createMagSurvey(xyzd, B)
# View the data and chose a profile
param = Simulator.ViewMagSurvey2D(survey, dobj)
param

# Define the parametric model interactively
model = Simulator.ViewPrism(param.result)
model

"""
Explanation: How do we define direction of an earth magnetic field?
Earth magnetic field is a vector. To define a vector we need to choose a coordinate system. We use right-handed system: 
- X (Easting), 
- Y (Northing), and 
- Z (Up). 
Here we consider an earth magnetic field ($\vec{B_0}$), of which intensity is one. To define this unit vector, we use inclinatino and declination:
- Declination: An angle from geographic North (Ng) (positive clockwise)
- Inclination: Vertical angle from the N-E plane (positive down)
<img src="https://github.com/geoscixyz/geosci-labs/raw/main/images/mag/earthfield.png?raw=true" style="width: 60%; height: 60%"> </img>
What's data: total field anomaly
We consider a typical form of magnetic data. To illustrate this we consider an suceptible object embedded in the earth. 
Based upon the earth magnetic field ($\vec{B}_0$), this object will generate anomalous magnetic field ($\vec{B}_A$). We define an unit vector $\hat{B}_0$ for the earth field as 
$$ \hat{B}_0 = \frac{\vec{B}_0}{|\vec{B}_0|}$$ 
We measure both earth and anomalous magnetic field such that
$$ \vec{B} = \vec{B}_0 + \vec{B}_A$$
Total field anomaly, $\triangle \vec{B}$ can be defined as
$$  |\triangle \vec{B}| = |\vec{B}|-|\vec{B}_E| $$ 
If $|\vec{B}|\ll|\vec{B}_E|$, then that is total field anomaly $\triangle \vec{B}$ is the projection of the anomalous field onto the direction of the earth field:
$$ |\triangle \vec{B}| \simeq \vec{B}_A \cdot \hat{B}_0=|\vec{B}_A|cos\theta$$ 
<img src="https://github.com/geoscixyz/geosci-labs/raw/main/images/mag/totalfieldanomaly.png?raw=true" style="width: 50%; height: 50%">
Define a 3D prism
Our model is a rectangular prism. Parameters to define this prism are given below:

dx: length in Easting (x) direction (meter)
dy: length in Northing (y) direction (meter)
dz: length in Depth (z) direction (meter) below the receiver
depth: top boundary of the prism (meter)
pinc: inclination of the prism (reference is a unit northing vector; degree)
pdec: declination of the prism (reference is a unit northing vector; degree)

You can also change the height of the survey grid above the ground
- rx_h: height of the grid (meter)
Green dots show a plane where we measure data.
End of explanation
"""


Simulator.PFSimulator(model, param)

"""
Explanation: Magnetic applet
Based on the prism that you made above, below Magnetic applet computes magnetic field at receiver locations, and provide both 2D map (left) and profile line (right). 
For the prism, you can alter:
- sus: susceptibility of the prism
Parameters for the earth field are:
- Einc: inclination of the earth field (degree)
- Edec: declination of the earth field (degree)
- Bigrf: intensity of the earth field (nT)
For data, you can view:
- tf: total field anomaly,
- bx :x-component, 
- by :y-component, 
- bz :z-component
You can simulate and view remanent magnetization effect with parameters:
- irt: "induced", "remanent", or "total"
- Q: Koenigsberger ratio ($\frac{M_{rem}}{M_{ind}}$)
- rinc: inclination of the remanent magnetization (degree)
- rdec: declination of the remanent magnetization (degree)
End of explanation
"""


