diff --git a/simpegDCIP/Dev/DC2D_fwr_obs.py b/simpegDCIP/Dev/DC2D_fwr_obs.py
index a83507bd..029695d3 100644
--- a/simpegDCIP/Dev/DC2D_fwr_obs.py
+++ b/simpegDCIP/Dev/DC2D_fwr_obs.py
@@ -2,31 +2,130 @@ import os
home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegdc\simpegDCIP\Dev'
-inpfile = 'MAG3Cfwr.inp'
-
os.chdir(home_dir)
#%%
-from SimPEG import np, Utils, Mesh
+from SimPEG import np, Utils, Mesh, mkvc, SolverLU
import simpegDCIP as DC
import pylab as plt
+#from ipywidgets import interact, IntSlider
+from matplotlib import animation
+from JSAnimation import HTMLWriter
from readUBC_DC2DMesh import readUBC_DC2DMesh
from readUBC_DC2DModel import readUBC_DC2DModel
from readUBC_DC2DLoc import readUBC_DC2DLoc
+# Number of padding cells to remove from plotting
+padc = 16
+
# Load UBC mesh 2D
-mesh = readUBC_DC2DMesh('mesh2d.txt')
+mesh = readUBC_DC2DMesh('mesh2d_fine.txt')
# Load model
-model = readUBC_DC2DModel('model2d.con')
+model = readUBC_DC2DModel('model2d_fine.con')
# load obs file
-[tx,rx] = readUBC_DC2DLoc('obs2d_East.loc')
+[txLoc,rxLoc,d,wd] = readUBC_DC2DLoc('obs2d_East.loc')
-#%% Plot model
-fig, axs = plt.subplots(1,1, figsize=(10,7))
-h1 = mesh.plotImage(model, ax = axs)
-plt.ylim([-1000,0])
-plt.xlim([0,2000])
+# Create SimPEG objects
+rx = DC.RxDipole(rxLoc[:,0], rxLoc[:,1])
+#tx = DC.SrcDipole([rx],txLoc[200,0],txLoc[200,1])
+# Create sub-mesh for plotting
+hx = mesh.hx
+hy = mesh.hy
+
+hx_sub = hx[padc:-padc]
+hy_sub = hy[padc:]
+mesh_sub = Mesh.TensorMesh([hx_sub,hy_sub],(hx_sub[0], -sum(hy_sub)))
+model_sub = model.reshape(mesh.nCy,mesh.nCx)
+model_sub = mkvc(model_sub[padc:,padc:-padc].T)
+
+xx = mesh_sub.vectorCCx
+yy = mesh_sub.vectorCCy
+
+#%% Solve
+txii = range(500,1500,200)
+jx_CC_sub = np.zeros((len(txii),mesh_sub.nCx,mesh_sub.nCy))
+jy_CC_sub = np.zeros((len(txii),mesh_sub.nCx,mesh_sub.nCy))
+
+
+#==============================================================================
+# line1, = mesh_sub.plotImage([], ax = axs)
+# line2, = axs.streamplot([], [], [], [],color='k')
+#==============================================================================
+
+
+
+#==============================================================================
+# def init():
+# line1.set_data([])
+# line2.set_data([], [], [], [])
+# return line1, line2,
+#
+# def animate(ii):
+#==============================================================================
+
+
+
+for ii in range(len(txii)):
+
+ fig, axs = plt.subplots(1,1, figsize=(14,7))
+
+ tx = DC.SrcDipole([rx],txii[ii],txii[ii])
+
+ survey = DC.SurveyDC([tx])
+ problem = DC.ProblemDC_CC(mesh)
+ problem.pair(survey)
+ problem.Solver = SolverLU
+
+ u1 = problem.fields(model)
+
+ Msig1 = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
+
+ j = Msig1*mesh.cellGrad*u1[tx, 'phi_sol']
+ j_CC = mesh.aveF2CCV*j
+
+ # Compute charge density solving div*grad*phi
+ Q = mesh.faceDiv*mesh.cellGrad*u1[tx, 'phi_sol']
+
+ jx_CC = j_CC[0:mesh.nC].reshape(mesh.nCy,mesh.nCx).T
+ jy_CC = j_CC[mesh.nC:].reshape(mesh.nCy,mesh.nCx).T
+
+ #%% Grab only the core for presentation
+ jx_CC_sub[ii,:,:] = jx_CC[padc:-padc,padc:]
+ jy_CC_sub[ii,:,:] = jy_CC[padc:-padc,padc:]
+ Q_sub = Q.reshape(mesh.nCy,mesh.nCx)
+ Q_sub = Q_sub[padc:,padc:-padc]
+
+ J_rho = np.sqrt(jx_CC_sub[ii,:,:]**2 + jy_CC_sub[ii,:,:]**2)
+ lw = np.log10(J_rho/J_rho.min())
+
+#==============================================================================
+# line1.set_data(model_sub)
+# line2.set_data(mesh_sub.vectorCCx, mesh_sub.vectorCCy, jx_CC_sub[2,:,:].T, jy_CC_sub[2,:,:].T)
+#==============================================================================
+ axs.imshow(Q_sub,alpha=0.75,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',vmin=-1e-2, vmax=1e-2)
+ axs.imshow(np.log10(model_sub.reshape(mesh_sub.nCy,mesh_sub.nCx)),alpha=0.2,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',cmap='gray')
+ axs.streamplot(mesh_sub.vectorCCx, mesh_sub.vectorCCy, jx_CC_sub[ii,:,:].T, jy_CC_sub[ii,:,:].T,color='k',linewidth = lw.T)
+ plt.ylim([-800,0])
+ plt.xlim([25,2000])
+ plt.show()
+
+#%% Create widget
+
+#def viewInv(msh,iteration):
+
+
+
+#, linewidth=lw.T
+#%%
+#interact(viewInv,msh = mesh_sub, iteration = IntSlider(min=0, max=len(txii)-1 ,step=1, value=0))
+# set embed_frames=True to embed base64-encoded frames directly in the HTML
+#==============================================================================
+# anim = animation.FuncAnimation(fig, animate, init_func=init,
+# frames=100, interval=20, blit=True)
+#
+# anim.save('animation.html', writer=HTMLWriter(embed_frames=True))
+#==============================================================================
diff --git a/simpegDCIP/Dev/JSAnimation/IPython_display.py b/simpegDCIP/Dev/JSAnimation/IPython_display.py
new file mode 100644
index 00000000..70dcb869
--- /dev/null
+++ b/simpegDCIP/Dev/JSAnimation/IPython_display.py
@@ -0,0 +1,92 @@
+from .html_writer import HTMLWriter
+from matplotlib.animation import Animation
+import matplotlib.pyplot as plt
+import tempfile
+import random
+import os
+
+
+__all__ = ['anim_to_html', 'display_animation']
+
+
+class _NameOnlyTemporaryFile(object):
+ """A context-managed temporary file which is not opened.
+
+ The file should be accessible by name on any system.
+
+ Parameters
+ ----------
+ suffix : string
+ The suffix of the temporary file (default = '')
+ prefix : string
+ The prefix of the temporary file (default = '_tmp_')
+ hash_length : string
+ The length of the random hash. The size of the hash space will
+ be 16 ** hash_length (default=8)
+ seed : integer
+ the seed for the random number generator. If not specified, the
+ system time will be used as a seed.
+ absolute : boolean
+ If true, return an absolute path to a temporary file in the current
+ working directory.
+
+ Example
+ -------
+
+ >>> with _NameOnlyTemporaryFile(seed=0, absolute=False) as f:
+ ... print(f)
+ ...
+ _tmp_d82c07cd
+ >>> os.path.exists('_tmp_d82c07cd') # file removed after context
+ False
+
+ """
+ def __init__(self, prefix='_tmp_', suffix='', hash_length=8,
+ seed=None, absolute=True):
+ rng = random.Random(seed)
+ self.name = '%s%0*x%s' % (prefix, hash_length,
+ rng.getrandbits(4 * hash_length), suffix)
+ if absolute:
+ self.name = os.path.abspath(self.name)
+
+ def __enter__(self):
+ return self
+
+ def __exit__(self, *exc_info):
+ if os.path.exists(self.name):
+ os.remove(self.name)
+
+
+def anim_to_html(anim, fps=None, embed_frames=True, default_mode='loop'):
+ """Generate HTML representation of the animation"""
+ if fps is None and hasattr(anim, '_interval'):
+ # Convert interval in ms to frames per second
+ fps = 1000. / anim._interval
+
+ plt.close(anim._fig)
+ if hasattr(anim, "_html_representation"):
+ return anim._html_representation
+ else:
+ # tempfile can't be used here: we need a filename, and this
+ # fails on windows. Instead, we use a custom filename generator
+ #with tempfile.NamedTemporaryFile(suffix='.html') as f:
+ with _NameOnlyTemporaryFile(suffix='.html') as f:
+ anim.save(f.name, writer=HTMLWriter(fps=fps,
+ embed_frames=embed_frames,
+ default_mode=default_mode))
+ html = open(f.name).read()
+
+ anim._html_representation = html
+ return html
+
+
+def display_animation(anim, **kwargs):
+ """Display the animation with an IPython HTML object"""
+ from IPython.display import HTML
+ return HTML(anim_to_html(anim, **kwargs))
+
+
+# This is the magic that makes animations display automatically in the
+# IPython notebook. The _repr_html_ method is a special method recognized
+# by IPython.
+Animation._repr_html_ = anim_to_html
diff --git a/simpegDCIP/Dev/JSAnimation/__init__.py b/simpegDCIP/Dev/JSAnimation/__init__.py
new file mode 100644
index 00000000..a684391d
--- /dev/null
+++ b/simpegDCIP/Dev/JSAnimation/__init__.py
@@ -0,0 +1 @@
+from .html_writer import HTMLWriter
diff --git a/simpegDCIP/Dev/JSAnimation/examples.py b/simpegDCIP/Dev/JSAnimation/examples.py
new file mode 100644
index 00000000..066ebdf4
--- /dev/null
+++ b/simpegDCIP/Dev/JSAnimation/examples.py
@@ -0,0 +1,97 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib import animation
+from JSAnimation import IPython_display
+
+def basic_animation(frames=100, interval=30):
+ """Plot a basic sine wave with oscillating amplitude"""
+ fig = plt.figure()
+ ax = plt.axes(xlim=(0, 10), ylim=(-2, 2))
+ line, = ax.plot([], [], lw=2)
+
+ x = np.linspace(0, 10, 1000)
+
+ def init():
+ line.set_data([], [])
+ return line,
+
+ def animate(i):
+ y = np.cos(i * 0.02 * np.pi) * np.sin(x - i * 0.02 * np.pi)
+ line.set_data(x, y)
+ return line,
+
+ return animation.FuncAnimation(fig, animate, init_func=init,
+ frames=frames, interval=interval)
+
+
+def lorenz_animation(N_trajectories=20, rseed=1, frames=200, interval=30):
+ """Plot a 3D visualization of the dynamics of the Lorenz system"""
+ from scipy import integrate
+ from mpl_toolkits.mplot3d import Axes3D
+ from matplotlib.colors import cnames
+
+ def lorentz_deriv(coords, t0, sigma=10., beta=8./3, rho=28.0):
+ """Compute the time-derivative of a Lorentz system."""
+ x, y, z = coords
+ return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
+
+ # Choose random starting points, uniformly distributed from -15 to 15
+ np.random.seed(rseed)
+ x0 = -15 + 30 * np.random.random((N_trajectories, 3))
+
+ # Solve for the trajectories
+ t = np.linspace(0, 2, 500)
+ x_t = np.asarray([integrate.odeint(lorentz_deriv, x0i, t)
+ for x0i in x0])
+
+ # Set up figure & 3D axis for animation
+ fig = plt.figure()
+ ax = fig.add_axes([0, 0, 1, 1], projection='3d')
+ ax.axis('off')
+
+ # choose a different color for each trajectory
+ colors = plt.cm.jet(np.linspace(0, 1, N_trajectories))
+
+ # set up lines and points
+ lines = sum([ax.plot([], [], [], '-', c=c)
+ for c in colors], [])
+ pts = sum([ax.plot([], [], [], 'o', c=c, ms=4)
+ for c in colors], [])
+
+ # prepare the axes limits
+ ax.set_xlim((-25, 25))
+ ax.set_ylim((-35, 35))
+ ax.set_zlim((5, 55))
+
+ # set point-of-view: specified by (altitude degrees, azimuth degrees)
+ ax.view_init(30, 0)
+
+ # initialization function: plot the background of each frame
+ def init():
+ for line, pt in zip(lines, pts):
+ line.set_data([], [])
+ line.set_3d_properties([])
+
+ pt.set_data([], [])
+ pt.set_3d_properties([])
+ return lines + pts
+
+ # animation function: called sequentially
+ def animate(i):
+ # we'll step two time-steps per frame. This leads to nice results.
+ i = (2 * i) % x_t.shape[1]
+
+ for line, pt, xi in zip(lines, pts, x_t):
+ x, y, z = xi[:i + 1].T
+ line.set_data(x, y)
+ line.set_3d_properties(z)
+
+ pt.set_data(x[-1:], y[-1:])
+ pt.set_3d_properties(z[-1:])
+
+ ax.view_init(30, 0.3 * i)
+ fig.canvas.draw()
+ return lines + pts
+
+ return animation.FuncAnimation(fig, animate, init_func=init,
+ frames=frames, interval=interval)
diff --git a/simpegDCIP/Dev/JSAnimation/html_writer.py b/simpegDCIP/Dev/JSAnimation/html_writer.py
new file mode 100644
index 00000000..063d7aa7
--- /dev/null
+++ b/simpegDCIP/Dev/JSAnimation/html_writer.py
@@ -0,0 +1,327 @@
+import os
+import sys
+import random
+import string
+import warnings
+if sys.version_info < (3, 0):
+ from cStringIO import StringIO as InMemory
+else:
+ from io import BytesIO as InMemory
+from matplotlib.animation import writers, FileMovieWriter
+from base64 import b64encode
+
+
+ICON_DIR = os.path.join(os.path.dirname(__file__), 'icons')
+
+
+class _Icons(object):
+ """This class is a container for base64 representations of the icons"""
+ icons = ['first', 'prev', 'reverse', 'pause', 'play', 'next', 'last']
+
+ def __init__(self, icon_dir=ICON_DIR, extension='png'):
+ self.icon_dir = icon_dir
+ self.extension = extension
+ for icon in self.icons:
+ setattr(self, icon,
+ self._load_base64('{0}.{1}'.format(icon, extension)))
+
+ def _load_base64(self, filename):
+ data = open(os.path.join(self.icon_dir, filename), 'rb').read()
+ return 'data:image/{0};base64,{1}'.format(self.extension,
+ b64encode(data).decode('ascii'))
+
+
+JS_INCLUDE = """
+
+"""
+
+
+DISPLAY_TEMPLATE = """
+
+
![]()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
![]()
+
+
+
+
+
+
+
+
+
+
+
+
+
+