diff --git a/doc/source/library/examples/gf_forward.rst b/doc/source/library/examples/gf_forward.rst
index 60a17139d..2ebd29251 100644
--- a/doc/source/library/examples/gf_forward.rst
+++ b/doc/source/library/examples/gf_forward.rst
@@ -14,70 +14,14 @@ such stores.
 Further API documentation for the utilized objects can be found at :class:`~pyrocko.gf.targets.Target`,
 :class:`~pyrocko.gf.seismosizer.LocalEngine` and :class:`~pyrocko.gf.seismosizer.DCSource`.
 
-:download:`gf_forward_example1.py </static/gf_forward_example1.py>`
+Download :download:`gf_forward_example1.py </static/gf_forward_example1.py>`
 
-::
-
-    from pyrocko.gf import LocalEngine, Target, DCSource
-    from pyrocko import trace
-    from pyrocko.gui_util import PhaseMarker
-
-    # We need a pyrocko.gf.Engine object which provides us with the traces
-    # extracted from the store. In this case we are going to use a local
-    # engine since we are going to query a local store.
-    engine = LocalEngine(store_superdirs=['/media/usb/gf_stores'])
-
-    # The store we are going extract data from:
-    store_id = 'crust2_dd'
-
-    # Define a list of pyrocko.gf.Target objects, representing the recording
-    # devices. In this case one station with a three component sensor will
-    # serve fine for demonstation.
-    channel_codes = 'ENZ'
-    targets = [
-        Target(
-            lat=10.,
-            lon=10.,
-            store_id=store_id,
-            codes=('', 'STA', '', channel_code))
-        for channel_code in channel_codes]
-
-    # Let's use a double couple source representation.
-    source_dc = DCSource(
-        lat=11.,
-        lon=11.,
-        depth=10000.,
-        strike=20.,
-        dip=40.,
-        rake=60.,
-        magnitude=4.)
-
-    # Processing that data will return a pyrocko.gf.Reponse object.
-    response = engine.process(source_dc, targets)
-
-    # This will return a list of the requested traces:
-    synthetic_traces = response.pyrocko_traces()
-
-    # In addition to that it is also possible to extract interpolated travel times
-    # of phases which have been defined in the store's config file.
-    store = engine.get_store(store_id)
-
-    markers = []
-    for t in targets:
-        dist = t.distance_to(source_dc)
-        depth = source_dc.depth
-        arrival_time = store.t('any_P', (depth, dist))
-        m = PhaseMarker(tmin=arrival_time,
-                        tmax=arrival_time,
-                        phasename='P',
-                        nslc_ids=(t.codes,))
-        markers.append(m)
-
-    # Finally, let's scrutinize these traces.
-    trace.snuffle(synthetic_traces, markers=markers)
+.. literalinclude :: /static/gf_forward_example1.py
+    :language: python
 
 .. figure :: /static/gf_synthetic.png
     :align: center
+    :width: 90%
     :alt: Synthetic seismograms calculated through pyrocko.gf
 
     Synthetic seismograms calculated through :class:`pyrocko.gf` displayed in :doc:`/apps/snuffler/index`. The three traces show the east, north and vertical synthetical displacement stimulated by a double-couple source at 155 km distance.
@@ -92,98 +36,15 @@ We will utilize :class:`~pyrocko.gf.seismosizer.LocalEngine`, :class:`~pyrocko.g
 
 .. figure:: /static/gf_static_displacement.png
     :align: center
+    :width: 90%
     :alt: Static displacement from a strike-slip fault calculated through pyrocko
 
     Synthetic surface displacement from a vertical strike-slip fault, with a N104W azimuth, in the Line-of-sight (LOS), east, north and vertical directions. LOS as for Envisat satellite (Look Angle: 23., Heading:-76). Positive motion toward the satellite. 
 
-:download:`gf_forward_example2.py </static/gf_forward_example2.py>`
-
-::
-
-    from pyrocko.gf import LocalEngine, SatelliteTarget, RectangularSource
-    import numpy as num
-
-    # distance in kilometer
-    km = 1e3
-
-    # Ignite the LocalEngine and point it to fomosto stores stored on a
-    # USB stick, for this example we use a static store with id 'static_store'
-    engine = LocalEngine(store_superdirs=['/media/usb/stores'])
-    store_id = 'static_store'
-
-    # We define an extended source, in this case a rectangular geometry
-    # Centroid UTM position is defined relatively to geographical lat, lon position
-    # Purely lef-lateral strike-slip fault with an N104W azimuth.
-    rect_source = RectangularSource(
-        lat=0., lon=0.,
-        north_shift=0., east_shift=0., depth=6.5*km,
-        width=5*km, length=8*km,
-        dip=90., rake=0., strike=104.,
-        slip=1.)
-
-    # We will define 1000 randomly distributed targets.
-    ntargets = 1000
-
-    # We initialize the satellite target and set the line of sight vectors
-    # direction, example of the Envisat satellite
-    look = 23.     # angle between the LOS and the vertical
-    heading = -76  # angle between the azimuth and the east (anti-clock)
-    theta = num.empty(ntargets)  # vertical LOS from horizontal
-    theta.fill(num.deg2rad(90. - look))
-    phi = num.empty(ntargets)  # horizontal LOS from E in anti-clokwise rotation
-    phi.fill(num.deg2rad(-90-heading))
-
-    satellite_target = SatelliteTarget(
-        north_shifts=(num.random.rand(ntargets)-.5) * 30. * km,
-        east_shifts=(num.random.rand(ntargets)-.5) * 30. * km,
-        tsnapshot=60,
-        interpolation='nearest_neighbor',
-        phi=phi,
-        theta=theta,
-        store_id=store_id)
-
-    # The computation is performed by calling process on the engine
-    result = engine.process(rect_source, [satellite_target])
-
-
-    def plot_static_los_result(result, target=0):
-        '''Helper function for plotting the displacement'''
-
-        import matplotlib.pyplot as plt
-
-        N = result.request.targets[target].coords5[:, 2]
-        E = result.request.targets[target].coords5[:, 3]
-        result = result.results_list[0][target].result
-
-        # get the component names
-        components = result.keys()
-        fig, _ = plt.subplots(int(len(components)/2), int(len(components)/2))
-
-        vranges = [(result[k].max(),
-                    result[k].min()) for k in components]
-
-        for dspl, ax, vrange in zip(components, fig.axes, vranges):
-
-            lmax = num.abs([num.min(vrange), num.max(vrange)]).max()
-            levels = num.linspace(-lmax, lmax, 50)
-
-            # plot interpolated points in map view with tricontourf
-            cmap = ax.tricontourf(E, N, result[dspl],
-                                  cmap=plt.get_cmap('seismic'), levels=levels)
-
-            ax.set_title(dspl+' [m]')
-            ax.set_aspect('equal')
-
-            # We plot the modeled fault
-            n, e = rect_source.outline(cs='xy').T
-            ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
-
-            fig.colorbar(cmap, ax=ax, aspect=5)
-
-        plt.show()
-
-    plot_static_los_result(result)
+Download :download:`gf_forward_example2.py </static/gf_forward_example2.py>`
 
+.. literalinclude :: /static/gf_forward_example2.py
+    :language: python
 
 Calculate forward model of thrust event and display wrapped phase
 -----------------------------------------------------------------
@@ -192,146 +53,16 @@ In this example we compare the synthetic unwappred and wrapped LOS displacements
 
 .. figure:: /static/gf_static_wrapper.png
     :align: center
+    :width: 90%
     :alt: Static displacement from a thrust fault calculated through pyrocko
 
     Synthetic LOS displacements from a south-dipping thrust fault. LOS as for Sentinel-1 satellite (Look Angle: 36., Heading:-76). Positive motion toward the satellite. Left: unwrapped phase. Right: Wrapped phase.
 
 
-:download:`gf_forward_example3.py </static/gf_forward_example3.py>`
-
-::
-
-    from pyrocko.gf import LocalEngine, SatelliteTarget, RectangularSource
-    import numpy as num
-
-    # distance in kilometer
-    km = 1e3
-
-    # Ignite the LocalEngine and point it to fomosto stores stored on a
-    # USB stick, for this example we use a static store with id 'static_store'
-    store_id = 'static_store'
-    engine = LocalEngine(store_superdirs=['/media/usb/gf_stores'])
-
-    # We want to reproduce the USGS Solution of the event
-    d, strike, dip, l, W, rake, slip = 10.5, 90., 40., 10., 10., 90., 5.
-
-    # We compute the magnitude of the event
-    potency = l*km*W*km*slip
-    m0 = potency*31.5e9
-    mw = (2./3) * num.log10(m0) - 6.07
-
-    # We define an extended source, in this case a rectangular geometry
-    # horizontal distance
-    # The centorid north position depends on its dip angle and its width.
-    n = num.cos(num.deg2rad(dip))*W/2
-
-    thrust = RectangularSource(
-        north_shift=n*km, east_shift=0.,
-        depth=d*km, width=W*km, length=l*km,
-        dip=dip, rake=rake, strike=strike,
-        slip=slip)
-
-    # We define a grid for the targets.
-    left, right, bottom, top = -15*km, 15*km, -15*km, 15*km
-    ntargets = 10000
-
-    # We initialize the satellite target and set the line of site vectors
-    # Case example of the Sentinel-1 satellite:
-    # Heading: -166 (anti clokwise rotation from east)
-    # Average Look Angle: 36 (from vertical)
-    heading = -76.
-    look = 36.
-    phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clokwise rotation
-    theta = num.empty(ntargets)  # Vertical LOS from horizontal
-    phi.fill(num.deg2rad(-90-heading))
-    theta.fill(num.deg2rad(90.-look))
-
-    satellite_target = SatelliteTarget(
-        north_shifts=num.random.uniform(bottom, top, ntargets),
-        east_shifts=num.random.uniform(left, right, ntargets),
-        tsnapshot=60,
-        interpolation='nearest_neighbor',
-        phi=phi,
-        theta=theta,
-        store_id=store_id)
-
-    # The computation is performed by calling process on the engine
-    result = engine.process(thrust, [satellite_target])
-
-
-    def plot_static_los_result(result, target=0):
-        '''Helper function for plotting the displacement'''
-        import matplotlib.pyplot as plt
-
-        # get forward model from engine
-        N = result.request.targets[target].coords5[:, 2]
-        E = result.request.targets[target].coords5[:, 3]
-        result = result.results_list[0][target].result
-
-        fig, _ = plt.subplots(1, 2, figsize=(8, 4))
-        fig.suptitle(
-            "thrust: depth={:0.2f}, l={}, w={:0.2f},strike={}, "
-            "rake={}, dip={}, slip={}\n"
-            "heading={}, look angle={}, Mw={:0.3f}"
-            .format(d, l, W, strike, rake, dip, slip, heading, look, mw),
-            fontsize=14,
-            fontweight='bold')
-
-        # Plot unwrapped LOS displacements
-        ax = fig.axes[0]
-        # We shift the relative LOS displacements
-        los = result['displacement.los'] - result['displacement.los'].min()
-        losrange = [(los.max(), los.min())]
-        losmax = num.abs([num.min(losrange), num.max(losrange)]).max()
-        levels = num.linspace(0, losmax, 50)
-
-        cmap = ax.tricontourf(
-            E, N, los,
-            cmap=plt.get_cmap('seismic'),
-            levels=levels)
+Download :download:`gf_forward_example3.py </static/gf_forward_example3.py>`
 
-        ax.set_title('los')
-        ax.set_aspect('equal')
-
-        # We plot the fault projection to the surface
-        n, e = thrust.outline(cs='xy').T
-        ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
-        # We underline the tip of the thrust
-        ax.plot(e[:2], n[:2], linewidth=2., color='black', alpha=0.5)
-
-        fig.colorbar(cmap, ax=ax, orientation='vertical', aspect=5, shrink=0.5)
-
-        # We plot wrapped phase
-        ax = fig.axes[1]
-        # We wrap the phase between 0 and 0.028 mm
-        wavelenght = 0.028
-        wrapped_los = num.mod(los, wavelenght)
-        levels = num.linspace(0, wavelenght, 50)
-
-        # ax.tricontour(E, N, wrapped_los,
-        #   map='gist_rainbow', levels=levels, colors='k')
-        cmap = ax.tricontourf(
-            E, N, wrapped_los,
-            cmap=plt.get_cmap('gist_rainbow'),
-            levels=levels,
-            interpolation='bicubic')
-
-        ax.set_xlim(left, right)
-        ax.set_ylim(bottom, top)
-
-        ax.set_title('wrapped los')
-        ax.set_aspect('equal')
-
-        # We plot the fault projection to the surface
-        n, e = thrust.outline(cs='xy').T
-        ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
-        # We underline the tiip of the fault
-        ax.plot(e[:2], n[:2], linewidth=2., color='black', alpha=0.5)
-
-        fig.colorbar(cmap, orientation='vertical', shrink=0.5, aspect=5)
-        plt.show()
-
-    plot_static_los_result(result)
+.. literalinclude :: /static/gf_forward_example3.py
+    :language: python
 
 
 Combining severals sources 
@@ -340,200 +71,13 @@ In this example we combine two rectangular sources and plot the forward model in
 
 .. figure:: /static/gf_static_several.png
     :align: center
+    :width: 90%
 
     Synthetic LOS displacements from a flower-structure made of one strike-slip
     fault and one thrust fault. LOS as for Sentinel-1 satellite (Look Angle:
     36., Heading:-76). Positive motion toward the satellite. 
 
-:download:`gf_forward_example4.py </static/gf_forward_example4.py>`
-
-::
-
-    from pyrocko.gf import LocalEngine, SatelliteTarget, RectangularSource
-    import numpy as num
-    from pyrocko import gf
-    from pyrocko.guts import List
-
-    km = 1e3
-
-
-    class CombiSource(gf.Source):
-        '''Composite source model.'''
-
-        discretized_source_class = gf.DiscretizedMTSource
-
-        subsources = List.T(gf.Source.T())
-
-        def __init__(self, subsources=[], **kwargs):
-            if subsources:
-
-                lats = num.array(
-                    [subsource.lat for subsource in subsources], dtype=num.float)
-                lons = num.array(
-                    [subsource.lon for subsource in subsources], dtype=num.float)
-
-                assert num.all(lats == lats[0]) and num.all(lons == lons[0])
-                lat, lon = lats[0], lons[0]
-
-                # if not same use:
-                # lat, lon = center_latlon(subsources)
-
-                depth = float(num.mean([p.depth for p in subsources]))
-                t = float(num.mean([p.time for p in subsources]))
-                kwargs.update(time=t, lat=float(lat), lon=float(lon), depth=depth)
-
-            gf.Source.__init__(self, subsources=subsources, **kwargs)
-
-        def get_factor(self):
-            return 1.0
-
-        def discretize_basesource(self, store, target=None):
-
-            dsources = []
-            t0 = self.subsources[0].time
-            for sf in self.subsources:
-                assert t0 == sf.time
-                ds = sf.discretize_basesource(store, target)
-                ds.m6s *= sf.get_factor()
-                dsources.append(ds)
-
-            return gf.DiscretizedMTSource.combine(dsources)
-
-    # distance in kilometer
-    km = 1e3
-    # We define a grid for the targets.
-    left, right, bottom, top = -10*km, 10*km, -10*km, 10*km
-    ntargets = 1000
-
-    # Ignite the LocalEngine and point it to fomosto stores stored on a
-    # USB stick, for this example we use a static store with id 'static_store'
-    store_id = 'static_store'
-    engine = LocalEngine(store_superdirs=['/media/usb/gf_stores'])
-
-    # We define two finite sources
-    # The first one is a purely vertical strike-slip fault
-    strikeslip = RectangularSource(
-        north_shift=0., east_shift=0.,
-        depth=6*km, width=4*km, length=10*km,
-        dip=90., rake=0., strike=90.,
-        slip=1.)
-
-    # The second one is a ramp connecting to the root of the strike-slip fault
-    # ramp north shift (n) and width (w) depend on its dip angle and on
-    # the strike slip fault width
-    n, w = 2/num.tan(num.deg2rad(45.)), 2.*(2./(num.sin(num.deg2rad(45.))))
-    thrust = RectangularSource(
-        north_shift=n*km, east_shift=0.,
-        depth=6*km, width=w*km, length=10*km,
-        dip=45., rake=90., strike=90.,
-        slip=0.5)
-
-    # We initialize the satellite target and set the line of site vectors
-    # Case example of the Sentinel-1 satellite:
-    # Heading: -166 (anti clokwise rotation from east)
-    # Average Look Angle: 36 (from vertical)
-    heading = -76
-    look = 36.
-    phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clokwise rotation
-    theta = num.empty(ntargets)  # Vertical LOS from horizontal
-    phi.fill(num.deg2rad(-90. - heading))
-    theta.fill(num.deg2rad(90. - look))
-
-    satellite_target = SatelliteTarget(
-        north_shifts=num.random.uniform(bottom, top, ntargets),
-        east_shifts=num.random.uniform(left, right, ntargets),
-        tsnapshot=60,
-        interpolation='nearest_neighbor',
-        phi=phi,
-        theta=theta,
-        store_id=store_id)
-
-    # We combine the two sources here
-    patches = [strikeslip, thrust]
-    sources = CombiSource(subsources=patches)
-
-    # The computation is performed by calling process on the engine
-    result = engine.process(sources, [satellite_target])
-
-
-    def plot_static_los_profile(result, strike, l, w, x0, y0):
-        import matplotlib.pyplot as plt
-        import matplotlib.cm as cm
-        import matplotlib.colors as mcolors
-        fig, _ = plt.subplots(1, 2, figsize=(8, 4))
-
-        # strike,l,w,x0,y0: strike, length, width, x, and y position
-        # of the profile
-        strike = num.deg2rad(strike)
-        # We define the parallel and perpendicular vectors to the profile
-        s = [num.sin(strike), num.cos(strike)]
-        n = [num.cos(strike), -num.sin(strike)]
-
-        # We define the boundaries of the profile
-        ypmax, ypmin = l/2, -l/2
-        xpmax, xpmin = w/2, -w/2
-
-        # We define the corners of the profile
-        xpro, ypro = num.zeros((7)), num.zeros((7))
-        xpro[:] = x0-w/2*s[0]-l/2*n[0], x0+w/2*s[0]-l/2*n[0], \
-            x0+w/2*s[0]+l/2*n[0], x0-w/2*s[0]+l/2*n[0], x0-w/2*s[0]-l/2*n[0], \
-            x0-l/2*n[0], x0+l/2*n[0]
-
-        ypro[:] = y0-w/2*s[1]-l/2*n[1], y0+w/2*s[1]-l/2*n[1], \
-            y0+w/2*s[1]+l/2*n[1], y0-w/2*s[1]+l/2*n[1], y0-w/2*s[1]-l/2*n[1], \
-            y0-l/2*n[1], y0+l/2*n[1]
-
-        # We get the forward model from the engine
-        N = result.request.targets[0].coords5[:, 2]
-        E = result.request.targets[0].coords5[:, 3]
-        result = result.results_list[0][0].result
-
-        # We first plot the surface displacements in map view
-        ax = fig.axes[0]
-        los = result['displacement.los']
-        # losrange = [(los.max(),los.min())]
-        # losmax = num.abs([num.min(losrange), num.max(losrange)]).max()
-        levels = num.linspace(los.min(), los.max(), 50)
-
-        cmap = ax.tricontourf(E, N, los, cmap=plt.get_cmap('seismic'), levels=levels)
-
-        for sourcess in patches:
-            fn, fe = sourcess.outline(cs='xy').T
-            ax.fill(fe, fn, color=(0.5, 0.5, 0.5), alpha=0.5)
-            ax.plot(fe[:2], fn[:2], linewidth=2., color='black', alpha=0.5)
-
-        # We plot the limits of the profile in map view
-        ax.plot(xpro[:], ypro[:], color='black', lw=1.)
-        # plot colorbar
-        fig.colorbar(cmap, ax=ax, orientation='vertical', aspect=5)
-        ax.set_title('Map view')
-        ax.set_aspect('equal')
-
-        # We plot displacements in profile
-        ax = fig.axes[1]
-        # We compute the perpandicular and parallel components in the profile basis
-        yp = (E-x0)*n[0]+(N-y0)*n[1]
-        xp = (E-x0)*s[0]+(N-y0)*s[1]
-        los = result['displacement.los']
-
-        # We select data encompassing the profile
-        index = num.nonzero(
-            (xp > xpmax) | (xp < xpmin) | (yp > ypmax) | (yp < ypmin))
-
-        ypp, losp = num.delete(yp, index), \
-            num.delete(los, index)
-
-        # We associate the same color scale to the scatter plot
-        norm = mcolors.Normalize(vmin=los.min(), vmax=los.max())
-        m = cm.ScalarMappable(norm=norm, cmap=plt.get_cmap('seismic'))
-        facelos = m.to_rgba(losp)
-        ax.scatter(
-            ypp, losp,
-            s=0.3, marker='o', color=facelos, label='LOS displacements')
-
-        ax.legend(loc='best')
-        ax.set_title('Profile')
-
-        plt.show()
+Download :download:`gf_forward_example4.py </static/gf_forward_example4.py>`
 
-    plot_static_los_profile(result, 110., 18*km, 5*km, 0., 0.)
+.. literalinclude :: /static/gf_forward_example4.py
+    :language: python
diff --git a/doc/source/library/examples/plotting.rst b/doc/source/library/examples/plotting.rst
index c7c67d171..9675a1e4d 100644
--- a/doc/source/library/examples/plotting.rst
+++ b/doc/source/library/examples/plotting.rst
@@ -22,48 +22,10 @@ Beachballs from moment tensors
 
 This example demonstrates how to create beachballs from (random) moment tensors.  
 
-::
-    
-    import random
-    import logging
-    import sys
-    from matplotlib import pyplot as plt
-    from pyrocko import beachball, moment_tensor as pmt
-    from pyrocko import util
+Download :download:`beachball_example01.py </static/beachball_example01.py>`
 
-    logger = logging.getLogger(sys.argv[0])
-
-    util.setup_logging()
-
-    fig = plt.figure(figsize=(10., 4.))
-    fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
-    axes = fig.add_subplot(1, 1, 1)
-
-    for i in xrange(200):
-
-        # create random moment tensor
-        mt = pmt.MomentTensor.random_mt()
-
-        try:
-            # create beachball from moment tensor
-            beachball.plot_beachball_mpl(
-                mt, axes,
-                # type of beachball: deviatoric, full or double couple (dc)
-                beachball_type='full',
-                size=random.random()*120.,
-                position=(random.random()*10., random.random()*10.),
-                alpha=random.random(),
-                linewidth=1.0)
-
-        except beachball.BeachballError, e:
-            logger.error('%s for MT:\n%s' % (e, mt))
-
-    axes.set_xlim(0., 10.)
-    axes.set_ylim(0., 10.)
-    axes.set_axis_off()
-    fig.savefig('beachball-example01.pdf')
-
-    plt.show()
+.. literalinclude :: /static/beachball_example01.py
+    :language: python
 
 .. figure :: /static/beachball-example01.png
     :align: center
@@ -76,31 +38,10 @@ This example demonstrates how to create beachballs from (random) moment tensors.
 This example shows how to plot a full, a deviatoric and a double-couple beachball
 for a moment tensor.
 
-::
+Download :download:`beachball_example03.py </static/beachball_example03.py>`
 
-    from matplotlib import pyplot as plt
-    from pyrocko import beachball, moment_tensor as pmt, plot
-    
-    fig = plt.figure(figsize=(4., 2.))
-    fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
-    axes = fig.add_subplot(1, 1, 1)
-    axes.set_xlim(0., 4.)
-    axes.set_ylim(0., 2.)
-    axes.set_axis_off()
-    
-    for i, beachball_type in enumerate(['full', 'deviatoric', 'dc']):
-        beachball.plot_beachball_mpl(
-                pmt.as_mt((124654616., 370943136., -6965434.0,
-                           553316224., -307467264., 84703760.0)),
-                axes,
-                beachball_type=beachball_type,
-                size=60.,
-                position=(i+1, 1),
-                color_t=plot.mpl_color('scarletred2'),
-                linewidth=1.0)
-    
-    fig.savefig('beachball-example03.pdf')
-    plt.show()
+.. literalinclude :: /static/beachball_example03.py
+    :language: python
 
 .. figure :: /static/beachball-example03.png
     :align: center
@@ -122,133 +63,11 @@ Creating the beachball this ways allows for finer control over their location
 based on their size (in display units) which allows for a round beachball even
 if the axis are not 1:1.
 
-::
-
-    from matplotlib import transforms, pyplot as plt
-    from pyrocko import beachball, gf
-
-    # create source object
-    source1 = gf.DCSource(depth=35e3, strike=0., dip=90., rake=0.)
-
-    # set size of beachball
-    sz = 20.
-    # set beachball offset in points (one point from each axis)
-    szpt = (sz / 2.) / 72. + 1. / 72.
-
-    fig = plt.figure(figsize=(10., 4.))
-    ax = fig.add_subplot(1, 1, 1)
-    ax.set_xlim(0, 10)
-    ax.set_ylim(0, 10)
-
-    # get the bounding point (left-top)
-    x0 = ax.get_xlim()[0]
-    y1 = ax.get_ylim()[1]
-
-    # create a translation matrix, based on the final figure size and
-    # beachball location
-    move_trans = transforms.ScaledTranslation(szpt, -szpt, fig.dpi_scale_trans)
-
-    # get the inverse matrix for the axis where the beachball will be plotted
-    inv_trans = ax.transData.inverted()
-
-    # set the bouding point relative to the plotted axis of the beachball
-    x0, y1 = inv_trans.transform(move_trans.transform(
-        ax.transData.transform((x0, y1))))
-
-    # plot beachball
-    beachball.plot_beachball_mpl(source1.pyrocko_moment_tensor(), ax,
-                                 beachball_type='full', size=sz,
-                                 position=(x0, y1), linewidth=1.)
-
-
-    # create source object
-    source2 = gf.RectangularExplosionSource(depth=35e3, strike=0., dip=90.)
-
-    # set size of beachball
-    sz = 30.
-    # set beachball offset in points (one point from each axis)
-    szpt = (sz / 2.) / 72. + 1. / 72.
-
-    # get the bounding point (right-upper)
-    x1 = ax.get_xlim()[1]
-    y1 = ax.get_ylim()[1]
-
-    # create a translation matrix, based on the final figure size and
-    # beachball location
-    move_trans = transforms.ScaledTranslation(-szpt, -szpt, fig.dpi_scale_trans)
-
-    # get the inverse matrix for the axis where the beachball will be plotted
-    inv_trans = ax.transData.inverted()
-
-    # set the bouding point relative to the plotted axis of the beachball
-    x1, y1 = inv_trans.transform(move_trans.transform(
-        ax.transData.transform((x1, y1))))
-
-    # plot beachball
-    beachball.plot_beachball_mpl(source2.pyrocko_moment_tensor(), ax,
-                                 beachball_type='full', size=sz,
-                                 position=(x1, y1), linewidth=1.)
-
-
-    # create source object
-    source3 = gf.CLVDSource(amplitude=35e6, azimuth=30., dip=30.)
-
-    # set size of beachball
-    sz = 40.
-    # set beachball offset in points (one point from each axis)
-    szpt = (sz / 2.) / 72. + 1. / 72.
+Download :download:`beachball_example02.py </static/beachball_example02.py>`
 
-    # get the bounding point (left-bottom)
-    x0 = ax.get_xlim()[0]
-    y0 = ax.get_ylim()[0]
+.. literalinclude :: /static/beachball_example02.py
+    :language: python
 
-    # create a translation matrix, based on the final figure size and
-    # beachball location
-    move_trans = transforms.ScaledTranslation(szpt, szpt, fig.dpi_scale_trans)
-
-    # get the inverse matrix for the axis where the beachball will be plotted
-    inv_trans = ax.transData.inverted()
-
-    # set the bouding point relative to the plotted axis of the beachball
-    x0, y0 = inv_trans.transform(move_trans.transform(
-        ax.transData.transform((x0, y0))))
-
-    # plot beachball
-    beachball.plot_beachball_mpl(source3.pyrocko_moment_tensor(), ax,
-                                 beachball_type='full', size=sz,
-                                 position=(x0, y0), linewidth=1.)
-
-    # create source object
-    source4 = gf.DoubleDCSource(depth=35e3, strike1=0., dip1=90., rake1=0.,
-                                strike2=45., dip2=74., rake2=0.)
-
-    # set size of beachball
-    sz = 50.
-    # set beachball offset in points (one point from each axis)
-    szpt = (sz / 2.) / 72. + 1. / 72.
-
-    # get the bounding point (right-bottom)
-    x1 = ax.get_xlim()[1]
-    y0 = ax.get_ylim()[0]
-
-    # create a translation matrix, based on the final figure size and
-    # beachball location
-    move_trans = transforms.ScaledTranslation(-szpt, szpt, fig.dpi_scale_trans)
-
-    # get the inverse matrix for the axis where the beachball will be plotted
-    inv_trans = ax.transData.inverted()
-
-    # set the bouding point relative to the plotted axis of the beachball
-    x1, y0 = inv_trans.transform(move_trans.transform(
-        ax.transData.transform((x1, y0))))
-
-    # plot beachball
-    beachball.plot_beachball_mpl(source4.pyrocko_moment_tensor(), ax,
-                                 beachball_type='full', size=sz,
-                                 position=(x1, y0), linewidth=1.)
-
-    fig.savefig('beachball-example02.pdf')
-    plt.show()
 
 .. figure :: /static/beachball-example02.png
     :align: center
@@ -271,71 +90,10 @@ The take-off angles needed can be computed with some help of the
 :py:mod:`pyrocko.cake` module. Azimuth and distance computations are done with
 functions from :py:mod:`pyrocko.orthodrome`.
 
-:download:`beachball_example04.py </static/beachball_example04.py>`
-
-::
-
-    import numpy as num
-    from matplotlib import pyplot as plt
-    from pyrocko import moment_tensor as pmt, beachball, cake, orthodrome
-
-    km = 1000.
-
-    # source position and mechanism
-    slat, slon, sdepth = 0., 0., 10.*km
-    mt = pmt.MomentTensor.random_dc()
-
-    # receiver positions
-    rdepth = 0.0
-    rlatlons = [(50., 10.), (60., -50.), (-30., 60.)]
-
-    # earth model and phase for takeoff angle computations
-    mod = cake.load_model('ak135-f-continental.m')
-    phases = cake.PhaseDef.classic('P')
-
-    # setup figure with aspect=1.0/1.0, ranges=[-1.1, 1.1]
-    fig = plt.figure(figsize=(2., 2.))  # size in inch
-    fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
-    axes = fig.add_subplot(1, 1, 1, aspect=1.0)
-    axes.set_axis_off()
-    axes.set_xlim(-1.1, 1.1)
-    axes.set_ylim(-1.1, 1.1)
-
-    projection = 'lambert'
-
-    beachball.plot_beachball_mpl(
-        mt, axes,
-        position=(0., 0.),
-        size=2.0,
-        color_t=(0.7, 0.4, 0.4),
-        projection=projection,
-        size_units='data')
-
-    for rlat, rlon in rlatlons:
-        distance = orthodrome.distance_accurate50m(slat, slon, rlat, rlon)
-        rays = mod.arrivals(
-            phases=cake.PhaseDef('P'),
-            zstart=sdepth, zstop=rdepth, distances=[distance*cake.m2d])
-
-        if not rays:
-            continue
-
-        takeoff = rays[0].takeoff_angle()
-        azi = orthodrome.azimuth(slat, slon, rlat, rlon)
-
-        # to spherical coordinates, r, theta, phi in radians
-        rtp = num.array([[1., num.deg2rad(takeoff), num.deg2rad(90.-azi)]])
-
-        # to 3D coordinates (x, y, z)
-        points = beachball.numpy_rtp2xyz(rtp)
-
-        # project to 2D with same projection as used in beachball
-        x, y = beachball.project(points, projection=projection).T
-
-        axes.plot(x, y, '+', ms=10., mew=2.0, mec='black', mfc='none')
-
-    fig.savefig('beachball-example04.png')
+Download :download:`beachball_example04.py </static/beachball_example04.py>`
 
+.. literalinclude :: /static/beachball_example04.py
+    :language: python
 
 .. figure :: /static/beachball-example04.png
     :align: center
diff --git a/doc/source/static/gf_forward_example2.py b/doc/source/static/gf_forward_example2.py
index a390072f2..37f090dc7 100644
--- a/doc/source/static/gf_forward_example2.py
+++ b/doc/source/static/gf_forward_example2.py
@@ -1,26 +1,39 @@
-from pyrocko.gf import LocalEngine, SatelliteTarget, RectangularSource
+from pyrocko import gf
 import numpy as num
 
 # distance in kilometer
 km = 1e3
+# many seconds make a day
+day= 24.*3600.
 
-# Ignite the LocalEngine and point it to fomosto stores stored on a
-# USB stick, for this example we use a static store with id 'static_store'
-engine = LocalEngine(store_superdirs=['/media/usb/stores'])
-store_id = 'static_store'
+# Ignite the LocalEngine and point it to your fomosto store, e.g. stored on a
+# USB stick, which for example has the id 'Abruzzo_Ameri_static_nearfield'
+# (download at http://kinherd.org:8080/gfws/static/stores/)
+engine = gf.LocalEngine(store_superdirs=['/media/usb/stores'])
+store_id = 'Abruzzo_Ameri_static_nearfield'
 
 # We define an extended source, in this case a rectangular geometry
 # Centroid UTM position is defined relatively to geographical lat, lon position
 # Purely lef-lateral strike-slip fault with an N104W azimuth.
-rect_source = RectangularSource(
+rect_source = gf.RectangularSource(
     lat=0., lon=0.,
     north_shift=0., east_shift=0., depth=6.5*km,
     width=5*km, length=8*km,
     dip=90., rake=0., strike=104.,
     slip=1.)
 
-# We will define 1000 randomly distributed targets.
-ntargets = 1000
+# We will define a grid of targets
+# number in east and north directions, and total
+ngrid = 80
+# extension from origin in all directions
+obs_size = 20.*km
+ntargets = ngrid**2
+# make regular line vector
+norths = num.linspace(-obs_size, obs_size, ngrid)
+easts  = num.linspace(-obs_size, obs_size, ngrid)
+# make regular grid
+norths2d = num.repeat(norths, len(easts))
+easts2d = num.tile(easts, len(norths))
 
 # We initialize the satellite target and set the line of sight vectors
 # direction, example of the Envisat satellite
@@ -31,10 +44,10 @@
 phi = num.empty(ntargets)  # horizontal LOS from E in anti-clokwise rotation
 phi.fill(num.deg2rad(-90-heading))
 
-satellite_target = SatelliteTarget(
-    north_shifts=(num.random.rand(ntargets)-.5) * 30. * km,
-    east_shifts=(num.random.rand(ntargets)-.5) * 30. * km,
-    tsnapshot=60,
+satellite_target = gf.SatelliteTarget(
+    north_shifts=norths2d,
+    east_shifts=easts2d,
+    tsnapshot=1.*day,
     interpolation='nearest_neighbor',
     phi=phi,
     theta=theta,
@@ -48,30 +61,31 @@ def plot_static_los_result(result, target=0):
     '''Helper function for plotting the displacement'''
 
     import matplotlib.pyplot as plt
-
+   
+    # get target coordinates and displacements from results 
     N = result.request.targets[target].coords5[:, 2]
     E = result.request.targets[target].coords5[:, 3]
-    result = result.results_list[0][target].result
+    synth_disp = result.results_list[0][target].result
 
-    # get the component names
-    components = result.keys()
+    # get the component names of displacements
+    components = synth_disp.keys()
     fig, _ = plt.subplots(int(len(components)/2), int(len(components)/2))
 
-    vranges = [(result[k].max(),
-                result[k].min()) for k in components]
+    vranges = [(synth_disp[k].max(),
+                synth_disp[k].min()) for k in components]
 
-    for dspl, ax, vrange in zip(components, fig.axes, vranges):
+    for comp, ax, vrange in zip(components, fig.axes, vranges):
 
         lmax = num.abs([num.min(vrange), num.max(vrange)]).max()
-        levels = num.linspace(-lmax, lmax, 50)
-
-        # plot interpolated points in map view with tricontourf
-        cmap = ax.tricontourf(E, N, result[dspl],
-                              cmap=plt.get_cmap('seismic'), levels=levels)
-
-        ax.set_title(dspl+' [m]')
+        # plot displacements at targets as colored points 
+        cmap = ax.scatter(E, N, c=synth_disp[comp], s = 10., marker = 's', 
+                          edgecolor='face', cmap='seismic',
+                          vmin= -1.5*lmax, vmax=1.5*lmax)
+        
+        ax.set_title(comp+' [m]')
         ax.set_aspect('equal')
-
+        ax.set_xlim(-obs_size, obs_size)
+        ax.set_ylim(-obs_size, obs_size)
         # We plot the modeled fault
         n, e = rect_source.outline(cs='xy').T
         ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
@@ -80,5 +94,4 @@ def plot_static_los_result(result, target=0):
 
     plt.show()
 
-
 plot_static_los_result(result)
diff --git a/doc/source/static/gf_forward_example3.py b/doc/source/static/gf_forward_example3.py
index 4adaf87e7..235fdc4c2 100644
--- a/doc/source/static/gf_forward_example3.py
+++ b/doc/source/static/gf_forward_example3.py
@@ -1,58 +1,68 @@
-from pyrocko.gf import LocalEngine, SatelliteTarget, RectangularSource
+from pyrocko import gf
 import numpy as num
 
 # distance in kilometer
 km = 1e3
+# many seconds make a day
+day= 24.*3600.
 
-# Ignite the LocalEngine and point it to fomosto stores stored on a
-# USB stick, for this example we use a static store with id 'static_store'
-store_id = 'static_store'
-engine = LocalEngine(store_superdirs=['/media/usb/gf_stores'])
+# Ignite the LocalEngine and point it to your fomosto store, e.g. stored on a
+# USB stick, which for example has the id 'Abruzzo_Ameri_static_nearfield'
+# (download at http://kinherd.org:8080/gfws/static/stores/)
+engine = gf.LocalEngine(store_superdirs=['/media/usb/stores'])
+store_id = 'Abruzzo_Ameri_static_nearfield'
 
-# We want to reproduce the USGS Solution of the event
-d, strike, dip, l, W, rake, slip = 10.5, 90., 40., 10., 10., 90., 5.
+# We want to reproduce the USGS Solution of an event, e.g.
+dep, strike, dip, leng, wid, rake, slip = 10.5, 90., 40., 10., 10., 90., .5
 
 # We compute the magnitude of the event
-potency = l*km*W*km*slip
-m0 = potency*31.5e9
+potency = leng*km*wid*km*slip
+rigidity = 31.5e9
+m0 = potency*rigidity
 mw = (2./3) * num.log10(m0) - 6.07
 
-# We define an extended source, in this case a rectangular geometry
-# horizontal distance
-# The centorid north position depends on its dip angle and its width.
-n = num.cos(num.deg2rad(dip))*W/2
-
-thrust = RectangularSource(
-    north_shift=n*km, east_shift=0.,
-    depth=d*km, width=W*km, length=l*km,
+# We define an extended rectangular source
+thrust = gf.RectangularSource(
+    north_shift=0., east_shift=0.,
+    depth=dep*km, width=wid*km, length=leng*km,
     dip=dip, rake=rake, strike=strike,
     slip=slip)
 
-# We define a grid for the targets.
-left, right, bottom, top = -15*km, 15*km, -15*km, 15*km
-ntargets = 10000
+# We will define a grid of targets
+# number in east and north directions, and total
+ngrid = 90
+# extension from origin in all directions
+obs_size = 20.*km
+ntargets = ngrid**2
+# make regular line vector
+norths = num.linspace(-obs_size, obs_size, ngrid)
+easts  = num.linspace(-obs_size, obs_size, ngrid)
+# make regular grid
+norths2d = num.repeat(norths, len(easts))
+easts2d = num.tile(easts, len(norths))
+
 
 # We initialize the satellite target and set the line of site vectors
 # Case example of the Sentinel-1 satellite:
-# Heading: -166 (anti clokwise rotation from east)
+# Heading: -166 (anti-clockwise rotation from east)
 # Average Look Angle: 36 (from vertical)
 heading = -76.
 look = 36.
-phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clokwise rotation
+phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clockwise rotation
 theta = num.empty(ntargets)  # Vertical LOS from horizontal
 phi.fill(num.deg2rad(-90-heading))
 theta.fill(num.deg2rad(90.-look))
 
-satellite_target = SatelliteTarget(
-    north_shifts=num.random.uniform(bottom, top, ntargets),
-    east_shifts=num.random.uniform(left, right, ntargets),
-    tsnapshot=60,
+satellite_target = gf.SatelliteTarget(
+    north_shifts=norths2d,
+    east_shifts=easts2d,
+    tsnapshot=1.*day,
     interpolation='nearest_neighbor',
     phi=phi,
     theta=theta,
     store_id=store_id)
 
-# The computation is performed by calling process on the engine
+# The computation is performed by calling 'process' on the engine
 result = engine.process(thrust, [satellite_target])
 
 
@@ -60,38 +70,42 @@ def plot_static_los_result(result, target=0):
     '''Helper function for plotting the displacement'''
     import matplotlib.pyplot as plt
 
-    # get forward model from engine
+    # get synthetic displacements and target coordinates from engine's 'result'
     N = result.request.targets[target].coords5[:, 2]
     E = result.request.targets[target].coords5[:, 3]
-    result = result.results_list[0][target].result
-
+    synth_disp = result.results_list[0][target].result
+    
+    # we get the fault projection to the surface for plotting 
+    n, e = thrust.outline(cs='xy').T
+    
     fig, _ = plt.subplots(1, 2, figsize=(8, 4))
     fig.suptitle(
-        "thrust: depth={:0.2f}, l={}, w={:0.2f},strike={}, "
-        "rake={}, dip={}, slip={}\n"
-        "heading={}, look angle={}, Mw={:0.3f}"
-        .format(d, l, W, strike, rake, dip, slip, heading, look, mw),
+        "fault: dep={:0.2f}, l={}, w={:0.2f},str={},"
+        "rake={}, dip={}, slip={}, Mw={:0.3f}\n"
+        "satellite: heading={}, look angle={}"
+        .format(dep, leng, wid, strike, rake, dip, slip, heading, look, mw),
         fontsize=14,
         fontweight='bold')
 
     # Plot unwrapped LOS displacements
     ax = fig.axes[0]
     # We shift the relative LOS displacements
-    los = result['displacement.los'] - result['displacement.los'].min()
+    los = synth_disp['displacement.los']
     losrange = [(los.max(), los.min())]
     losmax = num.abs([num.min(losrange), num.max(losrange)]).max()
-    levels = num.linspace(0, losmax, 50)
 
-    cmap = ax.tricontourf(
-        E, N, los,
-        cmap=plt.get_cmap('seismic'),
-        levels=levels)
+    cmap = ax.scatter(
+        E, N, c=los,
+        s = 10., marker = 's', 
+        edgecolor='face',
+        cmap=plt.get_cmap('seismic'), 
+        vmin= -1.*losmax, vmax=1.*losmax)
 
-    ax.set_title('los')
+    ax.set_title('line-of-sight displacement [m]')
     ax.set_aspect('equal')
-
-    # We plot the fault projection to the surface
-    n, e = thrust.outline(cs='xy').T
+    ax.set_xlim(-obs_size, obs_size)
+    ax.set_ylim(-obs_size, obs_size)
+    # plot fault outline
     ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
     # We underline the tip of the thrust
     ax.plot(e[:2], n[:2], linewidth=2., color='black', alpha=0.5)
@@ -100,29 +114,24 @@ def plot_static_los_result(result, target=0):
 
     # We plot wrapped phase
     ax = fig.axes[1]
-    # We wrap the phase between 0 and 0.028 mm
-    wavelenght = 0.028
-    wrapped_los = num.mod(los, wavelenght)
-    levels = num.linspace(0, wavelenght, 50)
-
-    # ax.tricontour(E, N, wrapped_los,
-    #   map='gist_rainbow', levels=levels, colors='k')
-    cmap = ax.tricontourf(
-        E, N, wrapped_los,
-        cmap=plt.get_cmap('gist_rainbow'),
-        levels=levels,
-        interpolation='bicubic')
-
-    ax.set_xlim(left, right)
-    ax.set_ylim(bottom, top)
-
-    ax.set_title('wrapped los')
+    # We simulate a C-band interferogram for this source
+    c_lambda = 0.056
+    insar_phase = -num.mod(los, c_lambda/2.)/(c_lambda/2.)*2.*num.pi - num.pi
+
+    cmap = ax.scatter(
+        E, N, c= insar_phase,
+        s = 10., marker = 's', 
+        edgecolor='face',
+        cmap=plt.get_cmap('gist_rainbow'))
+
+    ax.set_xlim(-obs_size, obs_size)
+    ax.set_ylim(-obs_size, obs_size)
+    ax.set_title('simulated interferogram')
     ax.set_aspect('equal')
-
-    # We plot the fault projection to the surface
-    n, e = thrust.outline(cs='xy').T
+    
+    # plot fault outline
     ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
-    # We underline the tiip of the fault
+    # We outline the top edge of the fault with a thick line
     ax.plot(e[:2], n[:2], linewidth=2., color='black', alpha=0.5)
 
     fig.colorbar(cmap, orientation='vertical', shrink=0.5, aspect=5)
diff --git a/doc/source/static/gf_forward_example4.py b/doc/source/static/gf_forward_example4.py
index 96703be73..f9915388b 100644
--- a/doc/source/static/gf_forward_example4.py
+++ b/doc/source/static/gf_forward_example4.py
@@ -1,10 +1,10 @@
-from pyrocko.gf import LocalEngine, SatelliteTarget, RectangularSource
 import numpy as num
 from pyrocko import gf
 from pyrocko.guts import List
 
 km = 1e3
 
+day= 24.*3600.
 
 class CombiSource(gf.Source):
     '''Composite source model.'''
@@ -56,12 +56,12 @@ def discretize_basesource(self, store, target=None):
 
 # Ignite the LocalEngine and point it to fomosto stores stored on a
 # USB stick, for this example we use a static store with id 'static_store'
-store_id = 'static_store'
-engine = LocalEngine(store_superdirs=['/media/usb/gf_stores'])
+engine = gf.LocalEngine(store_superdirs=['/media/usb/stores'])
+store_id = 'Abruzzo_Ameri_static_nearfield'
 
 # We define two finite sources
 # The first one is a purely vertical strike-slip fault
-strikeslip = RectangularSource(
+strikeslip = gf.RectangularSource(
     north_shift=0., east_shift=0.,
     depth=6*km, width=4*km, length=10*km,
     dip=90., rake=0., strike=90.,
@@ -71,7 +71,7 @@ def discretize_basesource(self, store, target=None):
 # ramp north shift (n) and width (w) depend on its dip angle and on
 # the strike slip fault width
 n, w = 2/num.tan(num.deg2rad(45.)), 2.*(2./(num.sin(num.deg2rad(45.))))
-thrust = RectangularSource(
+thrust = gf.RectangularSource(
     north_shift=n*km, east_shift=0.,
     depth=6*km, width=w*km, length=10*km,
     dip=45., rake=90., strike=90.,
@@ -79,19 +79,19 @@ def discretize_basesource(self, store, target=None):
 
 # We initialize the satellite target and set the line of site vectors
 # Case example of the Sentinel-1 satellite:
-# Heading: -166 (anti clokwise rotation from east)
+# Heading: -166 (anti clockwise rotation from east)
 # Average Look Angle: 36 (from vertical)
 heading = -76
 look = 36.
-phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clokwise rotation
+phi = num.empty(ntargets)    # Horizontal LOS from E in anti-clockwise rotation
 theta = num.empty(ntargets)  # Vertical LOS from horizontal
 phi.fill(num.deg2rad(-90. - heading))
 theta.fill(num.deg2rad(90. - look))
 
-satellite_target = SatelliteTarget(
+satellite_target = gf.SatelliteTarget(
     north_shifts=num.random.uniform(bottom, top, ntargets),
     east_shifts=num.random.uniform(left, right, ntargets),
-    tsnapshot=60,
+    tsnapshot=1.*day,
     interpolation='nearest_neighbor',
     phi=phi,
     theta=theta,
@@ -140,11 +140,10 @@ def plot_static_los_profile(result, strike, l, w, x0, y0):
     # We first plot the surface displacements in map view
     ax = fig.axes[0]
     los = result['displacement.los']
-    # losrange = [(los.max(),los.min())]
-    # losmax = num.abs([num.min(losrange), num.max(losrange)]).max()
     levels = num.linspace(los.min(), los.max(), 50)
 
-    cmap = ax.tricontourf(E, N, los, cmap=plt.get_cmap('seismic'), levels=levels)
+    cmap = ax.tricontourf(E, N, los, cmap=plt.get_cmap('seismic'), 
+                          levels=levels)
 
     for sourcess in patches:
         fn, fe = sourcess.outline(cs='xy').T
@@ -160,7 +159,7 @@ def plot_static_los_profile(result, strike, l, w, x0, y0):
 
     # We plot displacements in profile
     ax = fig.axes[1]
-    # We compute the perpandicular and parallel components in the profile basis
+    # We compute the perpendicular and parallel components in the profile basis
     yp = (E-x0)*n[0]+(N-y0)*n[1]
     xp = (E-x0)*s[0]+(N-y0)*s[1]
     los = result['displacement.los']
diff --git a/doc/source/static/gf_static_displacement.png b/doc/source/static/gf_static_displacement.png
index e974341cf..a7c995d06 100644
Binary files a/doc/source/static/gf_static_displacement.png and b/doc/source/static/gf_static_displacement.png differ
diff --git a/doc/source/static/gf_static_several.png b/doc/source/static/gf_static_several.png
index 8ab1e1f37..3aab5308e 100644
Binary files a/doc/source/static/gf_static_several.png and b/doc/source/static/gf_static_several.png differ
diff --git a/doc/source/static/gf_static_wrapper.png b/doc/source/static/gf_static_wrapper.png
index 67104c5a0..68de3cc27 100644
Binary files a/doc/source/static/gf_static_wrapper.png and b/doc/source/static/gf_static_wrapper.png differ