diff --git a/.gitignore b/.gitignore
index 15e8813d..f8a55e4a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,17 +6,23 @@ __pycache__
 .DS_Store
 
 # the default build locations
+/build
 /builds
 /products
 
 # the synthetic data from the mogi model
 /models/mogi/examples/**/*.txt
 /models/mogi/examples/**/*.csv
+/models/mogi/examples
 
 # the synthetic data from the cdm model
 /models/cdm/examples/**/*.txt
 /models/cdm/examples/**/*.csv
 
+# reverso examples folder
+/models/reverso/examples
+/models/reverso/examples/*.asc
+
 # auctex lint
 /doc/**/auto
 # tex lint
diff --git a/models/mogi/examples/mogi.pfg b/models/mogi/examples/mogi.pfg
deleted file mode 100644
index 0c941a56..00000000
--- a/models/mogi/examples/mogi.pfg
+++ /dev/null
@@ -1,79 +0,0 @@
-;
-; michael a.g. aïvázis
-; orthologue
-; (c) 1998-2020 all rights reserved
-;
-
-mogi:
-    ; shell
-    ; shell = mpi.shells.mpirun ; for running with mpi
-
-    controller:
-        ; saving the results of the calculation
-        archiver:
-            theta = theta.asc
-            sigma = sigma.asc
-            llk = posterior.asc
-
-    ; job layout
-    job:
-        tasks = 1       ; number of tasks per host
-        gpus = 0        ; number of gpus per task
-        chains = 2**12  ; number of chains per task
-        steps = 20      ; the length of each chain
-
-    ; event handlers
-    monitors:
-       prof = altar.bayesian.profiler
-
-    ; model configuration
-    model:
-        ; the mogi engine: native or fast
-        mode = fast
-
-        ; parameter sets
-        psets:
-            location = contiguous
-            depth = contiguous
-            source = contiguous
-            offsets = contiguous
-
-            location:
-                count = 2
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (-5000, 5000)
-                prior:
-                    support = (-5000, 5000)
-            depth:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (2000, 4000)
-                prior:
-                    support = (2000, 4000)
-            source:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (9, 11)
-                prior:
-                    support = (9, 11)
-            offsets:
-                count = 2
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (-0.1, 0.1)
-                prior:
-                    support = (-0.1, 0.1)
-
-
-; for parallel runs
-mpi.shells.mpirun # altar.plexus.shell:
-    extra = -mca btl self,tcp
-
-; end of file
diff --git a/models/mogi/examples/synthetic/mogi.py b/models/mogi/examples/synthetic/mogi.py
deleted file mode 100755
index 753ee361..00000000
--- a/models/mogi/examples/synthetic/mogi.py
+++ /dev/null
@@ -1,167 +0,0 @@
-#!/usr/bin/env python3
-# -*- python -*-
-# -*- coding: utf-8 -*-
-#
-# michael a.g. aïvázis <michael.aivazis@para-sim.com>
-#
-# (c) 2013-2020 parasim inc
-# (c) 2010-2020 california institute of technology
-# all rights reserved
-#
-
-
-# externals
-from math import sin, cos, pi as π
-# the framework
-import altar
-# my model
-import altar.models.mogi
-
-
-# app
-class Mogi(altar.application, family="altar.applications.mogi"):
-    """
-    A generator of synthetic data for Mogi sources
-    """
-
-    # user configurable state
-    x = altar.properties.float(default=0)
-    x.doc = "the x coördinate of the Mogi source"
-
-    y = altar.properties.float(default=0)
-    y.doc = "the y coördinate of the Mogi source"
-
-    d = altar.properties.float(default=3000)
-    d.doc = "the depth of the Mogi source"
-
-    dV = altar.properties.float(default=1e10)
-    dV.doc = "the strength of the Mogi source"
-
-    nu = altar.properties.float(default=.25)
-    nu.doc = "the Poisson ratio"
-
-
-    # protocol obligation
-    @altar.export
-    def main(self, *args, **kwds):
-        """
-        The main entry point
-        """
-        # compute the displacements
-        data, covariance = self.mogi()
-        # dump the displacements in a CSV file
-        data.write(uri="displacements.csv")
-        # and the covariance in an ascii file
-        covariance.save(filename=altar.primitives.path("cd.txt"))
-        # all done
-        return 0
-
-
-    # meta-methods
-    def __init__(self, **kwds):
-        # chain up
-        super().__init__(**kwds)
-        # create my stations
-        self.stations = self.makeStations()
-        # all done
-        return
-
-
-    # implementation details
-    def mogi(self):
-        """
-        Synthesize displacements for a grid of stations given a specific source location and
-        strength
-        """
-        # get the stations
-        stations = self.stations
-        # dedcue the number of observations
-        observations = len(stations)
-        # make a source
-        source = altar.models.mogi.source(x=self.x, y=self.y, d=self.d, dV=self.dV, nu=self.nu)
-
-        # observe all displacements from the same angle for now
-        theta = π/4 # the azimuthal angle
-        phi = π     # the polar angle
-        # build the common projection vector
-        s = sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)
-
-        # allocate a matrix to hold the projections
-        los = altar.matrix(shape=(observations,3))
-        # go through the observations
-        for obs in range(observations):
-            # store the LOS vector
-            los[obs, 0] = s[0]
-            los[obs, 1] = s[1]
-            los[obs, 2] = s[2]
-
-        # compute the displacements
-        u = source.displacements(locations=stations, los=los)
-
-        # prepare the dataset
-        # rows: one for each location
-        # columns: observation id, u.s, x, y, theta, phi
-        # observation id simulates data that come from different sources and therefore require
-        # a different offset
-        data = altar.models.mogi.data(name="displacements")
-
-        # go through the observation locations
-        for idx, (x,y) in enumerate(stations):
-            # make a new entry in the data sheet
-            observation = data.pyre_new()
-            # western stations
-            if x < 0:
-                # come from a different data set
-                observation.oid = 1
-            # than
-            else:
-                # eastern stations
-                observation.oid = 0
-
-            # record the location of this observation
-            observation.x = x
-            observation.y = y
-
-            # project the displacement
-            observation.d = u[idx]
-            # save the direction of the projection vector
-            observation.theta = theta
-            observation.phi = phi
-
-        # the length of the data sheet is the number of observations
-        observations = len(data)
-        # allocate a matrix for the data correlation
-        correlation = altar.matrix(shape=[observations]*2).zero()
-
-        # go through the observations
-        for idx, observation in enumerate(data):
-            # set the covariance to a fraction of the "observed" displacement
-            correlation[idx,idx] = 1.0 #.01 * observation.d
-
-        # all done
-        return data, correlation
-
-
-    def makeStations(self):
-        """
-        Create a set of station coordinate
-        """
-        # get some help
-        import itertools
-        # build a set of points on a grid
-        stations = itertools.product(range(-5,6), range(-5,6))
-        # and return it
-        return tuple((x*1000, y*1000) for x,y in stations)
-
-
-# bootstrap
-if __name__ == "__main__":
-    # instantiate
-    app = Mogi(name="mogi")
-    # invoke
-    status = app.run()
-    # share
-    raise SystemExit(status)
-
-
-# end of file
diff --git a/models/reverso/examples/localhost b/models/reverso/examples/localhost
deleted file mode 100644
index 08660ef4..00000000
--- a/models/reverso/examples/localhost
+++ /dev/null
@@ -1 +0,0 @@
-localhost slots=4
diff --git a/models/reverso/examples/reverso.pfg b/models/reverso/examples/reverso.pfg
deleted file mode 100644
index c8aae6be..00000000
--- a/models/reverso/examples/reverso.pfg
+++ /dev/null
@@ -1,107 +0,0 @@
-;
-; michael a.g. aïvázis (michael.aivazis@para-sim.com)
-; grace bato           (mary.grace.p.bato@jpl.nasa.gov)
-; eric m. gurrola      (eric.m.gurrola@jpl.nasa.gov)
-;
-; (c) 2013-2020 parasim inc
-; (c) 2010-2020 california institute of technology
-; all rights reserved
-
-
-reverso:
-    ; for running with MPI
-    ; shell = mpi.mpirun
-
-    controller:
-        ; saving the results of the calculation
-        archiver:
-            theta = theta.asc
-            sigma = sigma.asc
-            llk = posterior.asc
-
-    ; job layout
-    job:
-        hosts = 1       ; number of MPI hosts
-        tasks = 1       ; number of tasks per host
-        gpus = 0        ; number of gpus per task
-        chains = 2**12  ; number of chains per task
-        steps = 20      ; the length of each chain
-
-    monitors:
-       prof = altar.bayesian.profiler
-
-    ; model configuration
-    model:
-        ; the cdm engine: native or fast
-        mode = fast
-
-        ; parameter sets
-        psets:
-            Qin = contiguous
-            H_s = contiguous
-            H_d = contiguous
-            a_s = contiguous
-            a_d = contiguous
-            a_c = contiguous
-
-            Qin:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (.4, .8)
-                prior:
-                    support = (.4, .8)
-
-            H_s:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (2500, 3500)
-                prior:
-                    support = (2500, 3500)
-
-            H_d:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (3500, 4500)
-                prior:
-                    support = (3500, 4500)
-
-            a_s:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (1500, 2500)
-                prior:
-                    support = (1500, 2500)
-
-            a_d:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (1500, 2500)
-                prior:
-                    support = (1500, 2500)
-
-            a_c:
-                count = 1
-                prep = uniform
-                prior = uniform
-                prep:
-                    support = (1, 2)
-                prior:
-                    support = (1, 2)
-
-
-; for parallel runs
-mpi.shells.mpirun # cdm.shell:
-    hostfile = localhost
-    extra = -mca btl self,tcp
-
-; end of file
diff --git a/models/reverso/examples/synthetic/reverso.py b/models/reverso/examples/synthetic/reverso.py
deleted file mode 100755
index aa5463a7..00000000
--- a/models/reverso/examples/synthetic/reverso.py
+++ /dev/null
@@ -1,173 +0,0 @@
-#!/usr/bin/env python3
-# -*- python -*-
-# -*- coding: utf-8 -*-
-#
-# michael a.g. aïvázis <michael.aivazis@para-sim.com>
-#
-# (c) 2013-2020 parasim inc
-# (c) 2010-2020 california institute of technology
-# all rights reserved
-
-
-# externals
-import math
-# framework
-import altar
-# my model
-import altar.models.reverso
-
-
-# app
-class Reverso(altar.application, family="altar.applications.reverso"):
-    """
-    A generator of synthetic data for {reverso} sources
-    """
-
-
-    # public data
-    H_s = altar.properties.float(default=3.0e3)
-    H_s.doc = "depth of the shallow reservoir"
-
-    H_d = altar.properties.float(default=4.0e3)
-    H_d.doc = "depth of the deep reservoir"
-
-    a_s = altar.properties.float(default=2.0e3)
-    a_s.doc = "radius of the shallow magma reservoir"
-
-    a_d = altar.properties.float(default=2.2e3)
-    a_d.doc = "radius of the deep magma reservoir"
-
-    a_c = altar.properties.float(default=1.5)
-    a_c.doc = "radius of the hydraulic pipe connecting two magma reservoirs"
-
-    Qin = altar.properties.float(default=0.6)
-    Qin.doc = "basal magma inflow rate"
-
-    # physical parameters
-    G = altar.properties.float(default=20.0E9)
-    G.doc = "shear modulus, [Pa, kg-m/s**2]"
-
-    v = altar.properties.float(default=0.25)
-    v.doc = "Poisson's ratio"
-
-    mu = altar.properties.float(default=2000.0)
-    mu.doc = "viscosity [Pa-s]"
-
-    drho = altar.properties.float(default=300.0)
-    drho.doc = "density difference (ρ_r-ρ_m), [kg/m**3]"
-
-    g = altar.properties.float(default=9.81)
-    g.doc = "gravitational acceleration [m/s**2]"
-
-
-    # obligations
-    @altar.export
-    def main(self, *args, **kwds):
-        """
-        The main entry point
-        """
-        # build the data records
-        data = self.reverso()
-        # dump the displacements in a CSV file
-        data.write(uri="displacements.csv")
-        # all done
-        return 0
-
-
-    # meta methods
-    def __init__(self, **kwds):
-        # chain up
-        super().__init__(**kwds)
-        # generate my the observation locations and times
-        self.ticks = tuple(self.makeTicks())
-        # all done
-        return
-
-
-    # implementation details
-    def reverso(self):
-        """
-        The generator
-        """
-        # make a source
-        source = altar.models.reverso.source(
-            H_s=self.H_s, H_d=self.H_d, a_s=self.a_s, a_d=self.a_d, a_c=self.a_c,
-            Qin=self.Qin,
-            G=self.G, v=self.v, mu=self.mu, drho=self.drho, g=self.g,
-            )
-
-        # prep the dataset; the layout is baked in
-        # rows: one for each location, time
-        # columns: oid, t,x,y,  u.E,u.N,u.U, σ.E,σ.N,σ.U
-        #
-        # the observation id simulates observations from different sensors
-        data = altar.models.reverso.data(name="displacements")
-
-        # get the observation locations and times
-        ticks = self.ticks
-        # prime the displacement calculator
-        displacements = source.displacements(locations=ticks)
-
-        # compute the displacements
-        for i,((t,x,y), (u_r, u_Z)) in enumerate(zip(ticks, displacements)):
-            # make a new entry in the data sheer
-            rec = data.pyre_new()
-
-            # record the observation id
-            rec.oid = 0
-            # record time and location
-            rec.t = t
-            rec.x = x
-            rec.y = y
-
-            # find the polar angle of the vector to the observation location
-            phi = math.atan2(y,x)
-            # compute the E and N components
-            u_E = u_r * math.sin(phi)
-            u_N = u_r * math.cos(phi)
-
-            # record the displacements
-            rec.uE = u_E
-            rec.uN = u_N
-            rec.uZ = u_Z
-
-            # estimate the variance base on a 5% deviation from the mean value
-            rec.σE = max(0.05*u_E, .01)**2
-            rec.σN = max(0.05*u_N, .01)**2
-            rec.σZ = max(0.05*u_Z, .01)**2
-
-        # all done
-        return data
-
-
-    def makeTicks(self):
-        """
-        Generate times and locations for the observations
-        """
-        # get time
-        year = altar.units.time.year.value
-        # max time
-        tMax = 1 * year
-        # build the time value
-        for exponent in range(-6,1):
-            # compute the time mark
-            t = 10**exponent * tMax
-            # build the distance
-            for r in range(1000, 6000, 1000):
-                # assemble the tick mark
-                yield (t, r, 0)
-        # all done
-        return
-
-
-# bootstrap
-if __name__ == "__main__":
-    # instantiate
-    app = Reverso(name="reverso")
-    # invoke
-    status = app.run()
-    # share
-    raise SystemExit(status)
-
-
-# end of file
diff --git a/models/reverso/lib/libreverso/Source.cc b/models/reverso/lib/libreverso/Source.cc
index 905dcfa9..1cc8cb8b 100644
--- a/models/reverso/lib/libreverso/Source.cc
+++ b/models/reverso/lib/libreverso/Source.cc
@@ -80,7 +80,7 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
                 H_s, H_d, a_s, a_d, a_c,
                 Qin,
                 _G, _v, _mu, _drho, _g,
-                predicted);
+                predicted);        
     }
 
     // all done
diff --git a/models/reverso/lib/libreverso/reverso.cc b/models/reverso/lib/libreverso/reverso.cc
index aa26641e..600e3faa 100644
--- a/models/reverso/lib/libreverso/reverso.cc
+++ b/models/reverso/lib/libreverso/reverso.cc
@@ -10,6 +10,7 @@
 #include <array>
 #include <stdexcept>
 #include <gsl/gsl_matrix.h>
+
 // declarations
 #include "reverso.h"
 
@@ -57,7 +58,7 @@ reverso(int sample, const gsl_matrix * locations,
 
     // the analytic solution
     // the characteristic time constant (eq. 10)
-    auto tau = (8.0 * mu * std::pow(H_c, gamma_s) * gamma_d * k * std::pow(a_s,3))
+    auto tau = (8.0 * mu * H_c * gamma_s * gamma_d * k * std::pow(a_s,3))
         / (G * std::pow(a_c,4) * gamma_r);
 
     auto A  = gamma_d*k / gamma_r;
diff --git a/models/reverso/reverso/libreverso.py b/models/reverso/reverso/libreverso.py
index 88f62db9..ea8702ba 100644
--- a/models/reverso/reverso/libreverso.py
+++ b/models/reverso/reverso/libreverso.py
@@ -53,7 +53,7 @@ def REVERSO(locations,
 
     # the analytic solution
     # the characteristic time constant (eq. 10)
-    tau = (8.0 * mu * H_c**gamma_s * gamma_d * k * a_s**3) / (G * a_c**4 * gamma_r)
+    tau = (8.0 * mu * H_c *gamma_s * gamma_d * k * a_s**3) / (G * a_c**4 * gamma_r)
 
     A  = gamma_d*k / gamma_r
     A *= dPd0 - dPs0 + drho*g*H_c - 8*gamma_s*mu*Qin*H_c / (pi * a_c**4 * gamma_r)
diff --git a/models/reverso/reverso/libreverso.pyc b/models/reverso/reverso/libreverso.pyc
new file mode 100644
index 00000000..ff3d18fd
Binary files /dev/null and b/models/reverso/reverso/libreverso.pyc differ
diff --git a/models/reverso/reverso/libreverso_v1.py b/models/reverso/reverso/libreverso_v1.py
new file mode 100644
index 00000000..88645425
--- /dev/null
+++ b/models/reverso/reverso/libreverso_v1.py
@@ -0,0 +1,248 @@
+#!/usr/bin/env python3
+
+# -*- python -*-
+# -*- coding: utf-8 -*-
+#
+# mary grace bato <mary.grace.p.bato@jpl.nasa.gov>
+#
+# (c) 2019 jet propulsion laboratory
+# (c) 2019 california institute of technology
+# all rights reserved
+#
+
+# The Two-Reservoir Model
+# This model assumes an elastic half-space and incompressible magma. The two magma
+# reservoirs comprise a deep reservoir connected to a shallow reservoir by an
+# hydraulic pipe ["A two-magma chamber model as a source of deformation at Grimvsvotn
+# Volcano, Iceland by Reverso etal (2014) Journal of Geophysical Research: Solid Earth
+
+import numpy
+from matplotlib import pyplot
+
+def Urmat(Hs, Hd, gammas, gammad, r, G, a_s, a_d, v):
+    Rs = numpy.sqrt(r**2 + Hs**2)
+    Rd = numpy.sqrt(r**2 + Hd**2)
+
+    if gammas == 1.0:
+        alphas = 1.0
+    else:
+        alphas = 4.*Hs**2/(numpy.pi*Rs**2)
+
+    if gammad == 1.0:
+        alphad = 1.0
+    else:
+        alphad = 4.*Hd**2/(numpy.pi*Rd**2)
+
+    H = ([
+          [r*a_s**3*alphas*(1-v)/(G*(Hs**2+r**2)**1.5),
+           r*a_d**3*alphad*(1-v)/(G*(Hd**2+r**2)**1.5)]
+         ])
+
+    return H
+
+def Uzmat(Hs, Hd, gammas, gammad, r, G, a_s, a_d, v):
+    Rs = numpy.sqrt(r**2 + Hs**2)
+    Rd = numpy.sqrt(r**2 + Hd**2)
+
+    if gammas == 1.0:
+        alphas = 1.0
+    else:
+        alphas = 4.*Hs**2/(numpy.pi*Rs**2)
+
+    if gammad == 1.0:
+        alphad = 1.0
+    else:
+        alphad = 4.*Hd**2/(numpy.pi*Rd**2)
+
+    H = ([
+          [Hs*a_s**3*alphas*(1-v)/(G*(Hs**2+r**2)**1.5),
+           Hd*a_d**3*alphad*(1-v)/(G*(Hd**2+r**2)**1.5)]
+         ])
+
+    return H
+
+def losmat(Hs, Hd, gammas, gammad, x, y, G, a_s, a_d, v, theta, phi):
+    r = numpy.sqrt(x**2 + y**2)
+    Rs = numpy.sqrt(r**2 + Hs**2)
+    Rd = numpy.sqrt(r**2 + Hd**2)
+
+    if gammas == 1.0:
+        alhpas = 1.0
+    else:
+        alphas = (4.0*Hs**2)/(numpy.pi*Rs**2)
+
+    if gammad == 1.0:
+        alhpad = 1.0
+    else:
+        alphad = (4.0*Hd**2)/(numpy.pi*Rd**2)
+
+    # Constants
+    GAMMA = (1.0-v)/G
+    Ds = alphas * (a_s/Rs)**3
+    Dd = alphad * (a_d/Rd)**3
+
+    H = ([
+           [GAMMA*Ds*(numpy.sin(theta)*(numpy.sin(phi)*y -
+                      numpy.sin(theta)*numpy.cos(phi)*x) +
+                      numpy.cos(theta)*Hs)               ,
+            GAMMA*Dd*(numpy.sin(theta)*(numpy.sin(phi)*y -
+                     numpy.sin(theta)*numpy.cos(phi)*x)  +
+                      numpy.cos(theta)*Hd)
+           ]
+         ]
+        )
+
+    return H
+
+def main():
+    ## Physical parameters
+    # shear modulus, [Pa, kg-m/s**2]
+    G = 20.0E9
+    # Poisson's ratio
+    v = 0.25
+    # Viscosity [Pa-s]
+    mu = 2000.0
+    # Density difference (ρ_r-ρ_m), [kg/m**3]
+    drho = 300.0
+    # Gravitational acceleration [m/s**2]
+    g = 9.81
+
+    # Basal conditions
+    Qin = 0.6   # Basal magma inflow rate [m**3/s]
+
+    # Initial conditions
+    # Shallow reservoir overpressure at t=0 [Pa]
+    dPs0 = 0.0
+    # Deep reservoir overpressure at t=0 [Pa]
+    dPd0 = 0.0
+
+    # Geometry
+    # radius of the hydraulic pipe
+    ac = 1.5
+    # radius of the shallow reservoir
+    a_s = 2.0e3
+    # radius of the deep reservoir
+    a_d = 2.2e3
+    # ratio of the reservoir volumes
+    k = (a_d/a_s)**3
+    # depth of the deep reservoir
+    Hd = 4.0e3
+    # depth of the shallow reservoir
+    Hs = 3.0e3
+    # length of the hydraulic connection (no vertical extensions of the reservoirs? Fig 6.)
+    Hc = Hd - Hs
+    gammas = 8.0*(1.0-v)/(3.*numpy.pi)
+    gammad = 8.0*(1.0-v)/(3.*numpy.pi)
+
+    # time-step in seconds (1 day)
+    dt = 86400.0
+    # max time (1 year) in seconds
+    tmax = dt*365.0 *1.0
+
+    ## differential equation
+    # the time array in seconds
+    t = numpy.arange(0, tmax, dt)
+    # the time array as fraction of total duration
+    tfrac = t/tmax
+    ns = len(t)
+
+    ## Initialization
+    dPs = numpy.zeros(ns)
+    dPs[0] = dPs0
+    dPd = numpy.zeros(ns)
+    dPd[0] = dPd0
+
+    # Simplifying the equations
+    # Hydraulic strength. Eq. (10) ξ*γ_d*k/(γ_s+γ_d*k)
+    C1 = (G*ac**4)/(8*mu*Hc*a_s**3*gammas)
+    # A in eq (11) modified to incorporate initial overpressores
+    A1 = drho*g*Hc + dPd0 - dPs0
+    A2 = G*Qin / (gammad*numpy.pi*a_d**3)
+    C2 = gammas / (gammad*k)
+
+    for i in range(1, ns):
+        dPs[i] = dt*C1*(A1 + dPd[i-1] - dPs[i-1]) + dPs[i-1]
+        dPd[i] = A2*dt - C2*(dPs[i] - dPs[i-1])  + dPd[i-1]
+
+    pyplot.plot(tfrac, dPs/1.0e6, label='Differential')
+    pyplot.legend(loc=2, prop={'size':14}, framealpha=0.5)
+    pyplot.show()
+
+    pyplot.plot(tfrac, dPd/1.e6, label='Differential')
+    pyplot.legend(loc=2, prop={'size':14}, framealpha=0.5)
+    pyplot.show()
+
+    # Analytic solution
+    # Calculate the characteristic time constant
+    tau = (8.0*mu*Hc**gammas*gammad*k*a_s**3)/(G*ac**4*(gammas+gammad*k))     # eq. (10)
+
+    A =  gammad*k/(gammas + gammad*k)
+    A *= dPd0 - dPs0 + drho*g*Hc - 8.*gammas*mu*Qin*Hc/(numpy.pi*ac**4*(gammas+gammad*k))
+
+    f0 = A*(1. - numpy.exp(-t/tau))
+    f1 = G*Qin*t/(numpy.pi*a_s**3*(gammas+gammad*k))
+    dPs_anal = f0 + f1 + dPs0
+    dPd_anal = -f0*gammas/(gammad*k) + f1 + dPd0
+
+    # Comparing Analytical solution with the differential equation
+    pyplot.plot(tfrac, dPs/1.0e6, label='Differential')
+    pyplot.plot(tfrac, dPs_anal/1.0e6, ls='--', lw=6, alpha=0.6, label='Analytical')
+    pyplot.legend(loc=2, prop={'size':14}, framealpha=0.5)
+    pyplot.title('Shallow Overpressure (MPa)')
+    pyplot.show()
+
+    pyplot.plot(tfrac, dPd/1.0e6, label='Differential')
+    pyplot.plot(tfrac, dPd_anal/1.0e6, ls='--', lw=6, alpha=0.6, label='Analytical')
+    pyplot.legend(loc=2, prop={'size':14}, framealpha=0.5)
+    pyplot.title('Deep Overpressure (MPa)')
+    pyplot.show()
+
+
+    # Generate r-array of GNSS stations.
+    # r is the distance from the center of the volcano and the GNSS station (or InSAR points)
+    rr = numpy.arange(1000, 6000, 1000)
+
+    # Number of observations
+    nObs = len(rr)
+
+    # H-matrix for the radial displacement
+    H_Ur = numpy.squeeze([Urmat(Hs, Hd, gammas, gammad, r, G, a_s, a_d, v) for i, r in enumerate(rr)])
+    H_Uz = numpy.squeeze([Uzmat(Hs, Hd, gammas, gammad, r, G, a_s, a_d, v) for i, r in enumerate(rr)])
+
+    # Generate the corresponding displacements
+    # H-matrix for the radial displacement
+    Ur = numpy.squeeze([numpy.mat(H_Ur) * numpy.mat([[dPs[i]], [dPd[i]]]) for i in range(ns)])
+    Uz = numpy.squeeze([numpy.mat(H_Uz) * numpy.mat([[dPs[i]], [dPd[i]]]) for i in range(ns)])
+
+    #Plot radial displacement
+    for i, label in enumerate(rr):
+        pyplot.plot(tfrac, Ur[:,i], label='r={}'.format(label/1000.)+' km')
+    pyplot.legend()
+    pyplot.title('Radial Displacement (m)')
+    pyplot.show()
+
+    #Plot vertical displacement
+    for i, label in enumerate(rr):
+        pyplot.plot(tfrac, Uz[:,i], label='r={}'.format(label/1000.)+' km')
+    pyplot.legend()
+    pyplot.title('Vertical Displacement (m)')
+    pyplot.show()
+
+    # Synthetic InSAR dataset
+    # Incidence angle, theta
+    theta = numpy.pi * (41./180.)
+    # Azimuth angle, phi
+    phi = numpy.pi * (-169./180.)
+
+    # create meshgrid
+    x = numpy.arange(-5000, 5100, 1000)
+    y = x
+    X, Y = numpy.meshgrid(x, y)
+    H_los = [losmat(Hs, Hd, gammas, gammad, x, y, G, a_s, a_d, v, theta, phi)]
+
+
+if __name__ == "__main__":
+    status = main()
+    raise SystemExit(status)
+
+# end-of-file