From 08f2020ae98417b0f646af2bfd4a8059cb411921 Mon Sep 17 00:00:00 2001
From: Dave Marchant <dwfmarchant@gmail.com>
Date: Wed, 5 Feb 2014 23:47:21 -0800
Subject: [PATCH] Added analytics directory to Utils.  Finished example of
 TDEM_b.  Seems to be working....

---
 simpegEM/TDEM/TDEM_b.py        | 19 +++++++++++++++----
 simpegEM/Utils/Ana/TEM.py      | 12 ++++++++++++
 simpegEM/Utils/Ana/__init__.py |  1 +
 simpegEM/Utils/__init__.py     |  1 +
 4 files changed, 29 insertions(+), 4 deletions(-)
 create mode 100644 simpegEM/Utils/Ana/TEM.py
 create mode 100644 simpegEM/Utils/Ana/__init__.py

diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py
index 742eccbb..fe560063 100644
--- a/simpegEM/TDEM/TDEM_b.py
+++ b/simpegEM/TDEM/TDEM_b.py
@@ -63,11 +63,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
 if __name__ == '__main__':
     from SimPEG import *
     import simpegEM as EM
+    from simpegEM.Utils.Ana import hzAnalyticDipoleT
+    from scipy.constants import mu_0
+    import matplotlib.pyplot as plt
 
     cs = 5.
     ncx = 20
     ncy = 6
-    npad = 15
+    npad = 20
     hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
     hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
     mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
@@ -81,8 +84,16 @@ if __name__ == '__main__':
     dat = EM.TDEM.DataTDEM1D(txLoc=txLoc, txType=txType, rxLoc=rxLoc, rxType=rxType, timeCh=timeCh)
 
     prb = EM.TDEM.ProblemTDEM_b(mesh, model)
-    prb.setTimes([1e-5, 5e-5, 2.5e-4], [50, 50, 50])
-    prb.sigma = np.ones(mesh.nCz)
+    prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
+    prb.sigma = np.ones(mesh.nCz)*1e-8
+    prb.sigma[mesh.vectorCCz<0] = 0.1
     prb.pair(dat)
 
-    F = prb.field(prb.sigma)
\ No newline at end of file
+    bz_calc = dat.dpred(-999999)
+    bz_ana = mu_0*hzAnalyticDipoleT(dat.rxLoc[0], prb.times, prb.sigma[0])
+
+
+    plt.loglog(prb.times, np.abs(bz_calc.flatten()), label='TDEM_b')
+    plt.loglog(prb.times, np.abs(bz_ana), 'r', label='Analytic')
+    plt.legend()
+    plt.show()
diff --git a/simpegEM/Utils/Ana/TEM.py b/simpegEM/Utils/Ana/TEM.py
new file mode 100644
index 00000000..b20e8191
--- /dev/null
+++ b/simpegEM/Utils/Ana/TEM.py
@@ -0,0 +1,12 @@
+import numpy as np
+from scipy.constants import mu_0, pi
+from scipy.special import erf
+
+def hzAnalyticDipoleT(r, t, sigma):
+    theta = np.sqrt((sigma*mu_0)/(4*t))
+    tr = theta*r
+    etr = erf(tr)
+    t1 = (9/(2*tr**2) - 1)*etr
+    t2 = (1/np.sqrt(pi))*(9/tr + 4*tr)*np.exp(-tr**2)
+    hz = (t1 - t2)/(4*pi*r**3)
+    return hz
diff --git a/simpegEM/Utils/Ana/__init__.py b/simpegEM/Utils/Ana/__init__.py
new file mode 100644
index 00000000..4553f927
--- /dev/null
+++ b/simpegEM/Utils/Ana/__init__.py
@@ -0,0 +1 @@
+from TEM import hzAnalyticDipoleT
\ No newline at end of file
diff --git a/simpegEM/Utils/__init__.py b/simpegEM/Utils/__init__.py
index 94a53c6f..2e690cbb 100644
--- a/simpegEM/Utils/__init__.py
+++ b/simpegEM/Utils/__init__.py
@@ -1 +1,2 @@
 import Sources
+import Ana
\ No newline at end of file