diff --git a/README.md b/README.md
index 3182aaf..74cb084 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@ RfPy is a software to calculate single event-station receiver functions from the
 Installation, Usage, API documentation and scripts are described at 
 https://paudetseis.github.io/RfPy/.
 
-Authors: [`Pascal Audet`](https://www.uogeophysics.com/authors/admin/) (Developer and Maintainer) and [`Jeremy Gosselin`](https://www.uogeophysics.com/authors/gosselin/) (Contributor)
+Authors: [Pascal Audet](https://www.uogeophysics.com/authors/admin/) and [Wasja Bloch](https://www.eoas.ubc.ca/people/wasjabloch)
 <!-- #### Citing
 
 If you use `SplitPy` in your work, please cite the 
diff --git a/rfpy/binning.py b/rfpy/binning.py
index 4ed72b0..4fe4021 100644
--- a/rfpy/binning.py
+++ b/rfpy/binning.py
@@ -32,7 +32,8 @@
 from scipy.signal import hilbert
 
 
-def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
+def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False,
+        include_empty=False):
     """ 
     Function to stack receiver functions into (baz or slow) bins
     This can be done using a linear stack (i.e., simple
@@ -45,12 +46,17 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
         a single trace.
     stream2 : :class:`~obspy.core.Stream`
         Optionally stack a second stream in the same operation.
-    dbaz : int
-        Number of bazk-azimuth samples in bins
-    dslow : int
-        Number of slowness samples in bins
+    typ: str
+        Attribute to bin
+        'baz': backazimuth (degree)
+        'slow': Horizontal slowness (s/km)
+        'dist': epicentral distance (degree)
+    nbin : int
+        Number of equally sized bins within data range
     pws : bool
         Whether or not to perform phase-weighted stacking
+    include_empty : bool
+        Return empty bins as null traces (default omits them)
 
     Returns
     -------
@@ -58,6 +64,11 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
         Stream containing one or two stacked traces,
         depending on the number of input streams
 
+    Note
+    ----
+    Sets the following attributes of the stack:
+        nbin: Number of bins
+        index_list: Indices of constituent traces in source stream
     """
 
     if not typ in ['baz', 'slow', 'dist']:
@@ -95,6 +106,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
                 nb = 0
                 array = np.zeros(len(stream[0].data))
                 weight = np.zeros(len(stream[0].data), dtype=complex)
+                index_list = []
 
                 # Loop through stat
                 for j, tr in enumerate(stream):
@@ -104,13 +116,16 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
 
                         nb += 1
                         array += tr.data
-                        hilb = hilbert(tr.data)
-                        phase = np.arctan2(hilb.imag, hilb.real)
-                        weight += np.exp(1j*phase)
+                        index_list.append(j)
+
+                        if pws:
+                            hilb = hilbert(tr.data)
+                            phase = np.arctan2(hilb.imag, hilb.real)
+                            weight += np.exp(1j*phase)
                         
                         continue
 
-                if nb > 0:
+                if nb > 0 or include_empty:
 
                     # Average and update stats
                     array /= nb
@@ -118,19 +133,18 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
 
                     trace = Trace(header=stream[0].stats)
                     trace.stats.nbin = nb
+                    trace.stats.index_list = index_list
+
                     if typ == 'baz':
                         trace.stats.baz = bins[i]
                         trace.stats.slow = None
-                        trace.stats.nbin = nb
                     elif typ == 'slow': 
                         trace.stats.slow = bins[i]
                         trace.stats.baz = None
-                        trace.stats.nbin = nb
                     elif typ == 'dist':
                         trace.stats.dist = bins[i]
                         trace.stats.slow = None
                         trace.stats.baz = None
-                        trace.stats.nbin = nb
                     if not pws:
                         weight = np.ones(len(stream[0].data))
                     trace.data = weight*array
@@ -144,7 +158,8 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False):
     return final_stream
 
 
-def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False):
+def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, baz_range=(0, 360),
+        slow_range=(0.04, 0.08), pws=False, include_empty=False):
     """ 
     Function to stack receiver functions into back-azimuth and slowness bins.
     This can be done using a linear stack (i.e., simple
@@ -157,12 +172,18 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False):
         a single trace.
     stream2 : :class:`~obspy.core.Stream`
         Optionally stack a second stream in the same operation.
-    dbaz : int
-        Number of bazk-azimuth samples in bins
-    dslow : int
+    nbaz : int
+        Number of back-azimuth samples in bins
+    nslow : int
         Number of slowness samples in bins
+    baz_range : 2-tuple
+        Minimum and maximum of the back-azimuth bins
+    slow_range : 2-tuple
+        Minimum and maximum of the slowness bins
     pws : bool
         Whether or not to perform phase-weighted stacking
+    include_empty : bool
+        Return empty bins as null traces (default omits them)
 
     Returns
     -------
@@ -170,11 +191,16 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False):
         Stream containing one or two stacked traces,
         depending on the number of input streams
 
+    Note
+    ----
+    Sets the following attributes of the stack:
+        nbin: Number of bins
+        index_list: Indices of constituent traces in source stream
     """
 
     # Define back-azimuth and slowness bins
-    baz_bins = np.linspace(0, 360, nbaz)
-    slow_bins = np.linspace(0.04, 0.08, nslow)
+    baz_bins = np.linspace(*baz_range, nbaz)
+    slow_bins = np.linspace(*slow_range, nslow)
 
     # Extract baz and slowness
     baz = [stream1[i].stats.baz for i in range(len(stream1))]
@@ -185,59 +211,192 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False):
     islow = np.digitize(slow, slow_bins)
 
     final_stream = []
+    if stream2:
+        stream1 = [stream1, stream2]
+    else:
+        stream1 = [stream1]
+
+    for stream in stream1:
+        # try:
+        # Define empty streams
+        binned_stream = Stream()
+
+        # Loop through baz_bins
+        for i in range(nbaz):
+            for j in range(nslow):
+
+                nbin = 0
+                index_list = []
+                array = np.zeros(len(stream[0].data), dtype=type(stream[0].data[0]))
+                weight = np.zeros(len(stream[0].data), dtype=complex)
 
-    for stream in [stream1, stream2]:
-        try:
-            # Define empty streams
-            binned_stream = Stream()
-
-            # Loop through baz_bins
-            for i in range(nbaz):
-                for j in range(nslow):
-
-                    nbin = 0
-                    array = np.zeros(len(stream[0].data))
-                    weight = np.zeros(len(stream[0].data), dtype=complex)
+                # Loop all traces
+                for k, tr in enumerate(stream):
 
-                    # Loop all traces
-                    for k, tr in enumerate(stream):
+                    # If index of baz_bins is equal to ibaz
+                    if i == ibaz[k] and j == islow[k]:
 
-                        # If index of baz_bins is equal to ibaz
-                        if i == ibaz[k] and j == islow[k]:
+                        nbin += 1
+                        array += tr.data
+                        index_list.append(k)
 
-                            nbin += 1
-                            array += tr.data
-                            hilb = hilbert(tr.data)
+                        if pws:
+                            hilb = hilbert(np.real(tr.data))
                             phase = np.arctan2(hilb.imag, hilb.real)
                             weight += np.exp(1j*phase)
-                            
-                            continue
+                        
+                        continue
 
-                    if nbin > 0:
+                if nbin > 0 or include_empty:
 
-                        # Average and update stats
-                        array /= nbin
-                        weight = np.real(abs(weight/nbin))
+                    # Average and update stats
+                    array /= nbin
+                    weight = np.real(abs(weight/nbin))
 
-                        trace = Trace(header=stream[0].stats)
-                        trace.stats.baz = baz_bins[i]
-                        trace.stats.slow = slow_bins[j]
-                        trace.stats.nbin = nbin
+                    trace = Trace(header=stream[0].stats)
+                    trace.stats.baz = baz_bins[i]
+                    trace.stats.slow = slow_bins[j]
+                    trace.stats.nbin = nbin
+                    trace.stats.index_list = index_list
 
-                        if not pws:
-                            weight = np.ones(len(stream[0].data))
-                        trace.data = weight*array
-                        binned_stream.append(trace)
+                    if not pws:
+                        weight = np.ones(len(stream[0].data))
+                    trace.data = weight*array
+                    binned_stream.append(trace)
 
-            final_stream.append(binned_stream)
+        final_stream.append(binned_stream)
 
-        except:
-            continue
+        # except:
+        #     final_stream = [Stream()]
+        #     continue
 
     return final_stream
 
-def bin_all(stream1, stream2=None, pws=False):
+
+def stack(rfdatas, which='rf', pws=False,
+        attributes={'bin': None, 'slow': None, 'index_list': None}):
+    """ 
+    Function to stack receiver functions stored in a RFData list. This can
+    be done using a linear stack (i.e., simple mean), or using phase-weighted
+    stacking.
+
+    Parameters
+    ----------
+    rfdatas : list of :class:`rfpy.rfdata.RFdata`
+        List of RFData objects to be stacked
+    which : str
+        'rf' stack receiver functions
+        'specs' stack spectra
+    pws : bool
+        Whether or not to perform phase-weighted stacking (which=rf only)
+    attributes : dict
+        Attributes passed to the traces of stack
+
+    Returns
+    -------
+    stack : :class:`~obspy.core.Stream`
+        Stream containing stacked traces
+    """
+
+    nbin = len(rfdatas)
+
+    template = rfdatas[0].__dict__[which]
+    typ = float
+    if which == 'specs':
+        typ = complex
+
+    array = np.zeros((len(template), len(template[0].data)), dtype=typ)
+    weight = np.zeros((len(template), len(template[0].data)), dtype=complex)
+
+    for rfdata in rfdatas:
+        stream = rfdata.__dict__[which]
+
+        for n, tr in enumerate(stream):
+            array[n, :] += tr.data
+
+            if pws:
+                hilb = hilbert(np.real(tr.data))
+                phase = np.arctan2(hilb.imag, hilb.real)
+                weight[n, :] += np.exp(1j*phase)
+
+    if pws:
+        weight = np.real(abs(weight/nbin))
+        array *= weight
+    array /= nbin
+
+    stacks = Stream()
+    for n in range(len(template)):
+        stack = Trace(header=template[n].stats)
+        stack.data = array[n, :]
+        stack.stats.nbin = nbin
+        for att in attributes:
+            val = attributes[att]
+            stack.stats.__dict__[att] = val
+        stacks.append(stack)
+
+    return stacks
+
+
+def baz_slow_bin_indices(
+        rfs, nbaz=36+1, nslow=20+1, baz_range=(0, 360), slow_range=(0.04, 0.08), include_empty=False):
     """ 
+    Sort traces of streams into backazimuth (nbaz) and slowness (nslow) bins.
+
+    Parameters
+    ----------
+    rfs : list of :class:`~rfpy.RFData`
+        List of receiver functions to be sorted in bins
+    bazs : int
+        Number of back-azimuth samples in bins
+    slows : int
+        Number of slowness samples in bins
+    baz_range : 2-tuple
+        Minimum and maximum of the back-azimuth bins
+    slow_range : 2-tuple
+        Minimum and maximum of the slowness bins
+    include_empty : bool
+        Include empty bins (default omits them)
+
+    Returns
+    -------
+    index_lists : list of lists
+        Indices that sorts traces of stream into bins defined by nbaz and nslow
+    bazslow_list : list of 2*tuple
+        Backazimuth and slowness values of each element of index list
+    """
+
+    # Define back-azimuth and slowness bins
+    baz_bins = np.linspace(*baz_range, nbaz)
+    slow_bins = np.linspace(*slow_range, nslow)
+
+    # Extract baz and slowness
+    baz = [rf.meta.baz for rf in rfs]
+    slow = [rf.meta.slow for rf in rfs]
+
+    # Digitize baz and slowness
+    ibaz = np.digitize(baz, baz_bins)
+    islow = np.digitize(slow, slow_bins)
+
+    index_lists = []
+    bazslow_list = []
+
+    for i in range(nbaz):
+        for j in range(nslow):
+            index_list = []
+
+            for k in range(len(rfs)):
+                if i == ibaz[k] and j == islow[k]:
+                    index_list.append(k)
+
+            if len(index_list) > 0 or include_empty:
+                index_lists.append(index_list)
+                bazslow_list.append((baz_bins[i], slow_bins[j]))
+
+    return index_lists, bazslow_list
+
+
+def bin_all(stream1, stream2=None, pws=False):
+    """
     Function to bin all streams into a single trace.
     This can be done using a linear stack (i.e., simple
     mean), or using phase-weighted stacking.
@@ -274,9 +433,10 @@ def bin_all(stream1, stream2=None, pws=False):
             # Get phase weights
             for tr in stream:
                 array += tr.data
-                hilb = hilbert(tr.data)
-                phase = np.arctan2(hilb.imag, hilb.real)
-                pweight += np.exp(1j*phase)
+                if pws:
+                    hilb = hilbert(tr.data)
+                    phase = np.arctan2(hilb.imag, hilb.real)
+                    pweight += np.exp(1j*phase)
 
             # Normalize
             array = array/len(stream)
diff --git a/rfpy/plotting.py b/rfpy/plotting.py
index 11ac6c0..14db71a 100644
--- a/rfpy/plotting.py
+++ b/rfpy/plotting.py
@@ -29,12 +29,7 @@
 """
 
 # Import modules and functions
-import os
-import fnmatch
 import numpy as np
-import scipy as sp
-from obspy.core import Stream, Trace, AttribDict
-from scipy.interpolate import griddata
 import matplotlib.pyplot as plt
 
 
@@ -406,8 +401,8 @@ def wiggle_single_event(rfdata, filt=None, pre_filt=None, trange=None):
 
     lqtcopy = rfdata.data.copy()
     rfcopy = rfdata.rf.copy()
-    nn = lqtdata[0].stats.npts
-    sr = lqtdata[0].stats.sampling_rate
+    nn = lqtcopy[0].stats.npts
+    sr = lqtcopy[0].stats.sampling_rate
     taxis = np.arange(-nn/2., nn/2.)/sr
 
     if pre_filt:
@@ -447,6 +442,70 @@ def wiggle_single_event(rfdata, filt=None, pre_filt=None, trange=None):
     plt.show()
 
 
+def panel(stream1, stream2, sort=None, tmin=0., tmax=30., normalize=True,
+           save=False, title=None, form='png'):
+
+    taxis = stream1[0].stats.taxis
+    tselect = (taxis > tmin) & (taxis < tmax)
+    nt = len(taxis[tselect])
+
+    baz_array = np.array([[i, tr.stats.baz] for i, tr in enumerate(stream1)])
+    slow_array = np.array([[i, tr.stats.slow] for i, tr in enumerate(stream1)])
+
+    # build panel
+    panelV = np.zeros((len(stream1), nt))
+    for i, tr in enumerate(stream1):
+        panelV[i,:] = tr.data[tselect]
+    panelH = np.zeros((len(stream2), nt))
+    for i, tr in enumerate(stream2):
+        panelH[i,:] = tr.data[tselect]
+
+    if normalize:
+        vmin = np.mean(panelV) - 4.*np.std(panelV)
+        vmax = np.mean(panelV) + 4.*np.std(panelV)
+    else:
+        vmin = -0.5
+        vmax = 0.5
+
+    x, y = np.meshgrid(np.arange(len(stream1)), taxis[tselect])
+
+    plt.figure(figsize=(8,5))
+
+    plt.subplot(311)
+    plt.pcolor(x, y, panelV.T, shading='auto', cmap='bwr', vmin=vmin, vmax=vmax)
+    plt.ylim(tmin, tmax)
+    plt.gca().invert_yaxis()
+    plt.ylabel('Time (s)')
+    plt.gca().set_xticklabels([])
+
+    plt.subplot(312)
+    plt.pcolor(x, y, panelH.T, shading='auto', cmap='bwr', vmin=vmin, vmax=vmax)
+    plt.ylim(tmin, tmax)
+    plt.gca().invert_yaxis()
+    plt.ylabel('Time (s)')
+    plt.gca().set_xticklabels([])
+
+    plt.subplot(313)
+    plt.plot(baz_array[:,0], baz_array[:,1])
+    plt.xlim((0,len(stream1)))
+    plt.ylim((0,360))
+    plt.ylabel('Back-azimuth (deg)')
+    plt.xlabel('Receiver function #')
+    plt.twinx()
+    plt.plot(slow_array[:,0], slow_array[:,1], '.')
+    plt.ylabel('Slowness (s/km)')
+    plt.grid()
+
+    if title:
+        plt.suptitle(title)
+
+    if save:
+        plt.savefig('RF_PLOTS/' + stream1[0].stats.station +
+                    '.' + title + '.' + form, format=form)
+    else:
+        plt.show()
+
+
 def event_dist(stream, phase='P', save=False, title=None, form='png'):
 
     import cartopy.crs as ccrs
diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py
index 4caf39a..b7a08e7 100644
--- a/rfpy/rfdata.py
+++ b/rfpy/rfdata.py
@@ -22,6 +22,7 @@
 
 # -*- coding: utf-8 -*-
 from math import ceil
+import pickle
 import numpy as np
 from obspy import Trace, Stream, UTCDateTime
 from rfpy import utils
@@ -129,9 +130,9 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'):
             self.phase = phase
             self.accept = True
         else:
-            self.ttime = None
-            self.slow = None
-            self.inc = None
+            self.ttime = np.nan
+            self.slow = np.nan
+            self.inc = np.nan
             self.phase = None
             self.accept = False
 
@@ -142,16 +143,16 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'):
 
         # Attributes that get updated as analysis progresses
         self.rotated = False
-        self.snr = None
-        self.snrh = None
-        self.cc = None
+        self.snr = np.nan
+        self.snrh = np.nan
+        self.cc = np.nan
 
 
 class RFData(object):
     """
     A RFData object contains Class attributes that associate
     station information with a single event (i.e., earthquake) 
-    metadata, corresponding raw and rotated seismograms and 
+    metadata, corresponding raw and rotated seismograms, spectra and 
     receiver functions.
 
     Note
@@ -163,28 +164,37 @@ class RFData(object):
     ----------
     sta : object
         Object containing station information - from :mod:`~stdb` database.
+        When RFData is initialized without sta, functions requiring station info
+        coordinate will not be working.
     meta : :class:`~rfpy.rfdata.Meta`
         Object of metadata information for single event (initially set to None)
     data : :class:`~obspy.core.Stream`
         Stream object containing the three-component seismograms (either
         un-rotated or rotated by the method :func:`~rfpy.rfdata.rotate`)
-
+    rf : :class:`~obspy.core.Stream`
+        Stream object containing the receiver functions
+        (created when executing :func:`~rfpy.rfdata.deconvlove`)
+    specs : :class:`~obspy.core.Stream`
+        Stream object containing the spectra of the data
+        (created when executing :func:`~rfpy.rfdata.calc_spectra`)
     """
 
-    def __init__(self, sta):
+    def __init__(self, sta=None):
 
         # Load example data if initializing empty object
         if sta == 'demo' or sta == 'Demo':
             print("Uploading demo data - station NY.MMPY")
             import os
-            import pickle
             sta = pickle.load(
                 open(os.path.join(
                     os.path.dirname(__file__),
                     "examples/data", "MMPY.pkl"), 'rb'))['NY.MMPY']
 
         # Attributes from parameters
-        self.sta = sta
+        if sta:
+            self.sta = sta
+        else:
+            self.sta = None
 
         # Initialize meta and data objects as None
         self.meta = None
@@ -269,7 +279,6 @@ def add_data(self, stream, returned=False, new_sr=5.):
         # Load demo data
         if stream == 'demo' or stream == 'Demo':
             import os
-            import pickle
             file = open(os.path.join(
                 os.path.dirname(__file__),
                 "examples/data", "ZNE_Data.pkl"), "rb")
@@ -300,7 +309,9 @@ def add_data(self, stream, returned=False, new_sr=5.):
             return self.meta.accept
 
     def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
-                       new_sr=5.,dts=120., returned=False, verbose=False):
+                       new_sr=5.,dts=120., remove_response=False,
+                       local_response_dir='',
+                       returned=False, verbose=False):
         """
         Downloads seismograms based on event origin time and
         P phase arrival.
@@ -317,8 +328,13 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
             Time duration (sec)
         stdata : List
             Station list
+        remove_response : bool
+            Remove instrument response from seismogram and resitute to true ground
+            velocity (m/s) using obspy.core.trace.Trace.remove_response()
         returned : bool
             Whether or not to return the ``accept`` attribute
+        verbose : bool
+            Output diagnostics to screen
 
         Returns
         -------
@@ -332,6 +348,16 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
 
         """
 
+        def _resample(tr, new_sr):
+
+             if tr.stats.sampling_rate > new_sr:
+                 # only filter when downsampling
+                 tr.filter('lowpass', freq=0.5*new_sr, corners=2, zerophase=True)
+
+             tr.resample(new_sr, no_filter=True)
+             return tr
+
+
         if self.meta is None:
             raise(Exception("Requires event data as attribute - aborting"))
 
@@ -351,6 +377,8 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
         err, stream = utils.download_data(
             client=client, sta=self.sta, start=tstart, end=tend,
             stdata=stdata, dtype=dtype, ndval=ndval, new_sr=new_sr,
+            remove_response=remove_response,
+            local_response_dir=local_response_dir,
             verbose=verbose)
 
         # Store as attributes with traces in dictionary
@@ -360,13 +388,19 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
             trZ = stream.select(component='Z')[0]
             self.data = Stream(traces=[trZ, trN, trE])
 
-            # Filter Traces and resample
-            self.data.filter('lowpass', freq=0.5*new_sr,
-                             corners=2, zerophase=True)
+            # Filter Traces when downsampling
+            if self.data[0].stats.sampling_rate > new_sr:
+                self.data.filter('lowpass', freq=0.5*new_sr,
+                                 corners=2, zerophase=True)
+            # Resample
             self.data.resample(new_sr, no_filter=True)
 
         # If there is no ZNE, perhaps there is Z12?
-        except:
+        except Exception as e:
+            if verbose:
+                print('* Met exception: ')
+                print('* ' + str(e))
+                print('* Trying to access Z12 data.')
 
             try:
                 tr1 = stream.select(component='1')[0]
@@ -374,9 +408,11 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
                 trZ = stream.select(component='Z')[0]
                 self.data = Stream(traces=[trZ, tr1, tr2])
 
-                # Filter Traces and resample
-                self.data.filter('lowpass', freq=0.5*new_sr,
-                                 corners=2, zerophase=True)
+                # Filter Traces when downsampling
+                if self.data[0].stats.sampling_rate > new_sr:
+                    self.data.filter('lowpass', freq=0.5*new_sr,
+                                     corners=2, zerophase=True)
+                # Resample
                 self.data.resample(new_sr, no_filter=True)
 
                 # Save Z12 components in case it's necessary for later
@@ -385,7 +421,11 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan,
                 # Rotate from Z12 to ZNE using StDb azcorr attribute
                 self.rotate(align='ZNE')
 
-            except:
+            except Exception as e:
+                if verbose:
+                    print('* Met exception: ')
+                    print('* ' + str(e))
+                    print('* Cannot process')
                 self.meta.accept = False
 
         if returned:
@@ -416,6 +456,11 @@ def rotate(self, vp=None, vs=None, align=None):
         rotated : bool
             Whether or not the object has been rotated
 
+        Note
+        ----
+        'PVH' following Equation 18 in:
+        Kennett (1991). The removal of free surface interactions from three-component
+        seismograms. Geophysical Journal International
         """
 
         if not self.meta.accept:
@@ -552,9 +597,10 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.):
         if not self.meta.accept:
             return
 
-        if self.meta.snr:
-            print("SNR already calculated - continuing")
-            return
+        if self.meta.snr:  # False if None (for backward compatibility)
+            if np.isfinite(self.meta.snr):  # False if nan (None throws error)
+                print("SNR already calculated - continuing")
+                return
 
         t1 = self.meta.time + self.meta.ttime
 
@@ -612,12 +658,13 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.):
         # Calculate signal/noise ratio in dB
         self.meta.snrh = 10*np.log10(srms*srms/nrms/nrms)
 
-
-    def deconvolve(self, phase='P', vp=None, vs=None,
-                   align=None, method='wiener', pre_filt=None,
-                   gfilt=None, wlevel=0.01):
+    def calc_spectra(self, phase='P', vp=None, vs=None,
+                   align=None, method='wiener', wavelet='complete',
+                   envelope_threshold=0.05, time=5., pre_filt=None,
+                   writeto=False, norm=False, length=None):
         """
-        Deconvolves three-component data using one component as the source wavelet.
+        Calculate the numerators and denominator of the spectral division for receiver
+        function calculation using one component as the source wavelet.
         The source component is always taken as the dominant compressional 
         component, which can be either 'Z', 'L', or 'P'. 
 
@@ -631,26 +678,46 @@ def deconvolve(self, phase='P', vp=None, vs=None,
             Alignment of coordinate system for rotation
             ('ZRT', 'LQT', or 'PVH')
         method : str
-            Method for deconvolution. Options are 'wiener', 'water' or 
+            Method for deconvolution. Options are 'wiener', 'water', 'water2' or 
             'multitaper'
+        wavelet : str
+            Type of wavelet for deconvolution. Options are 'complete', 'time' or 
+            'envelope'
+        envelope_threshold : float
+            Threshold [0-1] used in ``wavelet='envelope'``.
+        time : float
+            Window length used in ``wavelet='time'``.
+            Minimum window length for ``wavelet='envelope'``.
         pre_filt : list of 2 floats
             Low and High frequency corners of bandpass filter applied
             before deconvolution
         gfilt : float
             Center frequency of Gaussian filter (Hz). 
         wlevel : float
-            Water level used in ``method='water'``.
+            Water level used in ``method='water'`` and ``method='water2'`` 
+        writeto : str or None
+            Write wavelets for deconvolution to file.
+        length : float or None
+            Time length (s) of the numerator. If None, use entire traces.
 
         Attributes
         ----------
-        rf : :class:`~obspy.core.Stream`
-            Stream containing the receiver function traces
+        specs : :class:`~obspy.core.Stream`
+            Stream containing the cross spectral quantities
 
         """
 
         if not self.meta.accept:
             return
 
+        if not align:
+            try:
+                align = self.meta.align
+            except AttributeError:
+                msg = "No alignment found. "
+                msg += "Please supply via align keyword."
+                raise ValueError(msg)
+
         def _npow2(x):
             return 1 if x == 0 else 2**(x-1).bit_length()
 
@@ -659,23 +726,94 @@ def _pad(array, n):
             tmp[:array.shape[0]] = array
             return tmp
 
-        def _gauss_filt(dt, nft, f0):
-            df = 1./(nft*dt)
-            nft21 = int(0.5*nft + 1)
-            f = df*np.arange(nft21)
-            w = 2.*np.pi*f
-            gauss = np.zeros(nft)
-            gauss[:nft21] = np.exp(-0.25*(w/f0)**2.)/dt
-            gauss[nft21:] = np.flip(gauss[1:nft21-1])
-            return gauss
+        def _Pwavelet(parent, method='complete', overhang=5.,
+                envelope_threshold=0.05, time=5.):
+
+            """
+            Select wavelet from the parent function for deconvolution using method.
+            parent: obspy.Trace
+                wavefrom to extract the wavelet from
+            method: str
+                'complete' use complete parent signal after P arrival  (current
+                    implementation)
+                'envelope' use only the part of the parent signal after the
+                    P arrival where
+                    envelope > envelope_threshold*max(envelope)
+                    fall back to 'complete' if condition not reached
+                'time' use only this many seconds after P arrival
+                    fall back to 'complete' if longer than parent
+                'noise' use the time before the P arrival
+            overhang: float
+                seconds before start and after end of wavelet to be used for
+                tapering
+            envelope_threshold: float
+                fraction of the envelope that defines wavelet (for
+                method='envelope')
+            time: float
+                window (seconds) that defines the wavelet (for method='time')
+                minimum time (seconds) of the wavelet (for method='envelope')
+
+            Return:
+            left, right: (obspy.UTCDateTime) start and end time of wavelet
+            """
+
+            import obspy.signal.filter as osf
+
+            if method not in ['complete', 'envelope', 'time', 'noise']:
+                msg = 'Unknow method for wavelet extraction: ' + method
+                raise NotImplementedError(msg)
+
+            errflag = False
+
+            if method == 'envelope':
+                split = int((self.meta.time + self.meta.ttime +
+                    time - parent.stats.starttime ) *
+                    parent.stats.sampling_rate)
+                env = osf.envelope(parent.data)
+                env /= max(env)  # normalize
+                env = env[split:]  # only look after P + time
+                try:
+                    i = np.nonzero(
+                            np.diff(
+                                np.array(
+                                    env > envelope_threshold, dtype=int))==-1)[0][0]
+                except IndexError:
+                    i = len(parent.data)-1
+                dts = i * parent.stats.delta + time
+                left = self.meta.time + self.meta.ttime - overhang
+                right = self.meta.time + self.meta.ttime + dts + overhang
+                if left < parent.stats.starttime or right > parent.stats.endtime:
+                    print('Envelope wavelet longer than trace.')
+                    print('Falling back to "complete" wavelet')
+                    errflag = True
+
+            if method == 'time':
+                left = self.meta.time + self.meta.ttime - overhang
+                right = self.meta.time + self.meta.ttime + time + overhang
+                if left < parent.stats.starttime or right > parent.stats.endtime:
+                    print('Time wavelet longer than trace.')
+                    print('Falling back to "complete" wavelet')
+                    errflag = True
+
+            if method == 'complete' or errflag:
+                dts = len(parent.data)*parent.stats.delta/2.
+                left = self.meta.time + self.meta.ttime - overhang
+                right = self.meta.time + self.meta.ttime + dts - 2*overhang
+            
+            if method == 'noise':
+                dts = len(parent.data)*parent.stats.delta/2.
+                left = self.meta.time + self.meta.ttime - dts
+                right = self.meta.time + self.meta.ttime - overhang
 
-        def _decon(parent, daughter1, daughter2, noise, nn, method):
+            return left, right
+            
+        def _calc_specs(parent, daughter1, daughter2, noise, nn, method):
 
-            # Get length, zero padding parameters and frequencies
+            # Sampling interval
             dt = parent.stats.delta
 
             # Wiener or Water level deconvolution
-            if method == 'wiener' or method == 'water':
+            if method == 'wiener' or method == 'water' or method == 'water2':
 
                 # npad = _npow2(nn*2)
                 npad = nn
@@ -693,14 +831,6 @@ def _decon(parent, daughter1, daughter2, noise, nn, method):
                 Sd2p = Fd2*np.conjugate(Fp)
                 Snn = np.real(Fn*np.conjugate(Fn))
 
-                # Final processing depends on method
-                if method == 'wiener':
-                    Sdenom = Spp + Snn
-                elif method == 'water':
-                    phi = np.amax(Spp)*wlevel
-                    Sdenom = Spp
-                    Sdenom[Sdenom < phi] = phi
-
             # Multitaper deconvolution
             elif method == 'multitaper':
 
@@ -740,41 +870,19 @@ def _decon(parent, daughter1, daughter2, noise, nn, method):
                 Sd2p = np.sum(Fd2*np.conjugate(Fp), axis=0)
                 Snn = np.sum(np.real(Fn*np.conjugate(Fn)), axis=0)
 
-                # Denominator
-                Sdenom = Spp + Snn
-
             else:
                 print("Method not implemented")
                 pass
 
-            # Apply Gaussian filter?
-            if gfilt:
-                gauss = _gauss_filt(dt, npad, gfilt)
-                gnorm = np.sum(gauss)*(freqs[1]-freqs[0])*dt
-            else:
-                gauss = np.ones(npad)
-                gnorm = 1.
-
-            # Copy traces
-            rfp = parent.copy()
-            rfd1 = daughter1.copy()
-            rfd2 = daughter2.copy()
-
-            # Spectral division and inverse transform
-            rfp.data = np.fft.ifftshift(np.real(np.fft.ifft(
-                gauss*Spp/Sdenom))/gnorm)
-            rfd1.data = np.fft.ifftshift(np.real(np.fft.ifft(
-                gauss*Sd1p/Sdenom))/gnorm)
-            rfd2.data = np.fft.ifftshift(np.real(np.fft.ifft(
-                gauss*Sd2p/Sdenom))/gnorm)
+            return Spp, Sd1p, Sd2p, Snn
 
-            return rfp, rfd1, rfd2
 
         if not self.meta.rotated:
             print("Warning: Data have not been rotated yet - rotating now")
             self.rotate(vp=vp, vs=vs, align=align)
 
-        if not self.meta.snr:
+        #  v--True if None      v--True if nan, error if None
+        if not self.meta.snr or not np.isfinite(self.meta.snr):
             print("Warning: SNR has not been calculated - " +
                   "calculating now using default")
             self.calc_snr()
@@ -795,68 +903,226 @@ def _decon(parent, daughter1, daughter2, noise, nn, method):
         trNl = self.data.select(component=cL)[0].copy()
         trNq = self.data.select(component=cQ)[0].copy()
 
+        trL.stats.channel = 'WV' + self.meta.align[0]
+        trQ.stats.channel = 'WV' + self.meta.align[1]
+        trT.stats.channel = 'WV' + self.meta.align[2]
+        trNl.stats.channel = 'WV' + self.meta.align[0]
+        trNq.stats.channel = 'WV' + self.meta.align[1]
+
         if phase == 'P' or 'PP':
 
             # Get signal length (i.e., seismogram to deconvolve) from trace length
-            dts = len(trL.data)*trL.stats.delta/2.
-            nn = int(round((dts-5.)*trL.stats.sampling_rate)) + 1
+            over = 5
+            dtsqt = len(trL.data) * trL.stats.delta / 2.
+
+            if length and length > dtsqt:
+                msg = f"Warning: Requested length {length} longer than data {dtsqt} - passing"
+                raise ValueError(msg)
+
+            if length is not None:
+                dtsqt = length
+
+            # Traces will be zero-paded to this length (samples)
+            nn = int(round((dtsqt + over) * trL.stats.sampling_rate)) + 1
 
-            # Signal window (-5. to dts-10 sec)
-            sig_left = self.meta.time+self.meta.ttime-5.
-            sig_right = self.meta.time+self.meta.ttime+dts-10.
+            sig_left, sig_right  = _Pwavelet(trL, method=wavelet,
+                envelope_threshold=envelope_threshold, time=time, overhang=over)
+
+            # Trim wavelet
+            trL.trim(sig_left, sig_right, nearest_sample=False, pad=True,
+                fill_value=0.)
+
+            # Signal window (-5. to dtsqt-10 sec)
+            if length is None:
+                sig_left, sig_right  = _Pwavelet(trQ, method='complete', overhang=over)
+            else:
+                sig_left, sig_right  = _Pwavelet(trQ, method='time', time=length, overhang=over)
 
             # Trim signal traces
             [tr.trim(sig_left, sig_right, nearest_sample=False, 
-                pad=nn, fill_value=0.) for tr in [trL, trQ, trT]]
+                pad=True, fill_value=0.) for tr in [trQ, trT]]
 
-            # Noise window (-dts to -5. sec)
-            noise_left = self.meta.time+self.meta.ttime-dts
-            noise_right = self.meta.time+self.meta.ttime-5.
+            # Noise window (-dtsqt to -5. sec)
+            noise_left, noise_right  = _Pwavelet(trQ, method='noise', overhang=over)
 
             # Trim noise traces
             [tr.trim(noise_left, noise_right, nearest_sample=False, 
-                pad=nn, fill_value=0.) for tr in [trNl, trNq]]
+                pad=True, fill_value=0.) for tr in [trNl, trNq]]
 
-        elif phase == 'S' or 'SKS':
+        # Can't deal with S phases right now ####
+        # elif phase == 'S' or 'SKS':
 
-            # Get signal length (i.e., seismogram to deconvolve) from trace length
-            dts = len(trL.data)*trL.stats.delta/2.
-
-            # Trim signal traces (-5. to dts-10 sec)
-            trL.trim(self.meta.time+self.meta.ttime+25.-dts/2.,
-                     self.meta.time+self.meta.ttime+25.)
-            trQ.trim(self.meta.time+self.meta.ttime+25.-dts/2.,
-                     self.meta.time+self.meta.ttime+25.)
-            trT.trim(self.meta.time+self.meta.ttime+25.-dts/2.,
-                     self.meta.time+self.meta.ttime+25.)
+        #     # Get signal length (i.e., seismogram to deconvolve) from trace length
+        #     dts = len(trL.data)*trL.stats.delta/2.
+
+        #     # Trim signal traces (-5. to dts-10 sec)
+        #     trL.trim(self.meta.time+self.meta.ttime+25.-dts/2.,
+        #              self.meta.time+self.meta.ttime+25.)
+        #     trQ.trim(self.meta.time+self.meta.ttime+25.-dts/2.,
+        #              self.meta.time+self.meta.ttime+25.)
+        #     trT.trim(self.meta.time+self.meta.ttime+25.-dts/2.,
+        #              self.meta.time+self.meta.ttime+25.)
             
-            # Trim noise traces (-dts to -5 sec)
-            trNl.trim(self.meta.time+self.meta.ttime-dts,
-                      self.meta.time+self.meta.ttime-dts/2.)
-            trNq.trim(self.meta.time+self.meta.ttime-dts,
-                      self.meta.time+self.meta.ttime-dts/2.)
+        #     # Trim noise traces (-dts to -5 sec)
+        #     trNl.trim(self.meta.time+self.meta.ttime-dts,
+        #               self.meta.time+self.meta.ttime-dts/2.)
+        #     trNq.trim(self.meta.time+self.meta.ttime-dts,
+        #               self.meta.time+self.meta.ttime-dts/2.)
 
         # Taper traces - only necessary processing after trimming
+        # TODO: What does this to the multitaper method
         [tr.taper(max_percentage=0.05, max_length=2.) 
          for tr in [trL, trQ, trT, trNl, trNq]]
 
-        # Pre-filter waveforms before deconvolution
+        # Pre-filter waveforms before calculating spectra
         if pre_filt:
             [tr.filter('bandpass', freqmin=pre_filt[0], freqmax=pre_filt[1],
                        corners=2, zerophase=True) 
              for tr in [trL, trQ, trT, trNl, trNq]]
 
-        # Deconvolve
+        if norm:
+            # Standardize all trace data - useful in simdec
+            meanN = np.mean(trNl.data)
+            stdN = np.std(trNl.data)
+            [(tr.data - meanN)/stdN for tr in [trL, trQ, trT, trNl, trNq]]
+
+        if writeto:
+            with open(writeto, 'wb') as f:
+                pickle.dump(Stream(traces=[trL, trQ, trT]), f)
+
+        # Calculate spectra
         if phase == 'P' or 'PP':
-            rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method)
+            spp, sd1p, sd2p, snn = _calc_specs(trL, trQ, trT, trNl, nn, method)
+
+        # Can't deal with S phases right now ####
+        # elif phase == 'S' or 'SKS':
+        #     rfQ, rfL, rfT = _calc_specs(trQ, trL, trT, trNq, nn, method)
+
+        # Copy traces
+        sLL = trL.copy()
+        sQL = trQ.copy()
+        sTL = trT.copy()
+        sNN = trNl.copy()
+        sLL.data = spp
+        sQL.data = sd1p
+        sTL.data = sd2p
+        sNN.data = snn
+
+        # Update stats of streams
+        sLL.stats.channel = 'S' + self.meta.align[0] + self.meta.align[0]
+        sQL.stats.channel = 'S' + self.meta.align[1] + self.meta.align[0]
+        sTL.stats.channel = 'S' + self.meta.align[2] + self.meta.align[0]
+        sNN.stats.channel = 'SNN'
 
-        elif phase == 'S' or 'SKS':
-            rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method)
+        self.specs = Stream(traces=[sLL, sQL, sTL, sNN])
+
+
+    def deconvolve(self, align=None, method='wiener',
+                   gfilt=None, wlevel=0.01):
+        """
+        Deconvolves three-component data using one component as the source wavelet.
+        The source component is always taken as the dominant compressional 
+        component, which can be either 'Z', 'L', or 'P'. 
+
+        Parameters
+        ----------
+        align : str
+            Alignment of coordinate system for rotation
+            ('ZRT', 'LQT', or 'PVH')
+        method : str
+            Method for deconvolution. Options are 'wiener', 'water' (aka 'wlevel'),
+            'water2' (aka 'wlevel2'), or 'multitaper'
+        gfilt : float
+            Center frequency of Gaussian filter (Hz). 
+        wlevel : float
+            Water level used in ``method='water'``.
+
+        Attributes
+        ----------
+        rf : :class:`~obspy.core.Stream`
+            Stream containing the receiver function traces
+
+        """
+
+        def _gauss_filt(dt, nft, f0):
+            df = 1./(nft*dt)
+            nft21 = int(0.5*nft + 1)
+            f = df*np.arange(nft21)
+            w = 2.*np.pi*f
+            gauss = np.zeros(nft)
+            gauss[:nft21] = np.exp(-0.25*(w/f0)**2.)/dt
+            gauss[nft21:] = np.flip(gauss[1:nft21-1])
+            return gauss
+
+
+        if not hasattr(self, 'specs'):
+            msg = "Spectra have not been calculated."
+            msg += "Call RFData.calc_spectra() first."
+            raise Exception(msg)
+
+        if not align:
+            try:
+                align = self.meta.align
+            except AttributeError:
+                msg = "No alignment found. "
+                msg += "Please supply via align keyword."
+                raise ValueError(msg)
+
+        # Make everything explicit
+        SLL = self.specs[0].copy()
+        SQL = self.specs[1].copy()
+        STL = self.specs[2].copy()
+        SNN = self.specs[3].copy()
+
+        dt = SLL.stats.delta
+        npad = SLL.stats.npts
+        freqs = np.fft.fftfreq(npad, d=dt)
+
+        # Wiener or Water level deconvolution
+        if method in ['wiener', 'multitaper']:
+            # Denominator (Spp + Snn)
+            Sdenom = SLL.data + SNN.data
+
+        elif method == 'water' or method == 'wlevel':
+            phi = np.amax(SLL.data)*wlevel
+            Sdenom = SLL.data
+            Sdenom[Sdenom < phi] = phi
+
+        elif method == 'water2' or method == 'wlevel2':
+            # David Gubbins
+            # Time Series Analysis and Inverse Therory for Geophysicists
+            # "Wiener Filter", Eq. 10.21
+            beta = np.amax(SLL.data)*wlevel
+            Sdenom = SLL.data + beta
+
+        else:
+            raise ValueError('Unknown method: ' + method)
+
+        # Apply Gaussian filter?
+        if gfilt:
+            gauss = _gauss_filt(dt, npad, gfilt)
+            gnorm = np.sum(gauss)*(freqs[1]-freqs[0])*dt
+        else:
+            gauss = np.ones(npad)
+            gnorm = 1.
+
+        # Copy traces
+        rfL = SLL.copy()
+        rfQ = SQL.copy()
+        rfT = STL.copy()
+
+        # Spectral division and inverse transform
+        rfL.data = np.fft.ifftshift(np.real(np.fft.ifft(
+            gauss*SLL.data/Sdenom))/gnorm)
+        rfQ.data = np.fft.ifftshift(np.real(np.fft.ifft(
+            gauss*SQL.data/Sdenom))/gnorm)
+        rfT.data = np.fft.ifftshift(np.real(np.fft.ifft(
+            gauss*STL.data/Sdenom))/gnorm)
 
         # Update stats of streams
-        rfL.stats.channel = 'RF' + self.meta.align[0]
-        rfQ.stats.channel = 'RF' + self.meta.align[1]
-        rfT.stats.channel = 'RF' + self.meta.align[2]
+        rfL.stats.channel = 'RF' + align[0]
+        rfQ.stats.channel = 'RF' + align[1]
+        rfT.stats.channel = 'RF' + align[2]
 
         self.rf = Stream(traces=[rfL, rfQ, rfT])
 
@@ -898,7 +1164,7 @@ def calc_cc(self):
         self.meta.cc = np.corrcoef(obs_Q.data, pred_Q.data)[0][1]
 
 
-    def to_stream(self):
+    def to_stream(self, store=None):
         """
         Method to switch from RFData object to Stream object.
         This allows easier manipulation of the receiver functions
@@ -909,7 +1175,7 @@ def to_stream(self):
         if not self.meta.accept:
             return
 
-        def _add_rfstats(trace):
+        def _add_stats(trace, store):
             trace.stats.snr = self.meta.snr
             trace.stats.snrh = self.meta.snrh
             trace.stats.cc = self.meta.cc
@@ -923,18 +1189,42 @@ def _add_rfstats(trace):
             trace.stats.vp = self.meta.vp
             trace.stats.vs = self.meta.vs
             trace.stats.phase = self.meta.phase
-            trace.stats.is_rf = True
-            nn = self.rf[0].stats.npts
-            sr = self.rf[0].stats.sampling_rate
+            trace.stats.align = self.meta.align
+            if store == 'rf':
+                trace.stats.is_rf = True
+                nn = self.rf[0].stats.npts
+                sr = self.rf[0].stats.sampling_rate
+            elif store == 'specs':
+                trace.stats.is_specs = True
+                nn = self.specs[0].stats.npts
+                sr = self.specs[0].stats.sampling_rate
+            else:
+                nn = self.data[0].stats.npts
+                sr = self.data[0].stats.sampling_rate
             trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn)
+
             return trace
 
-        if not hasattr(self, 'rf'):
-            raise(Exception("Warning: Receiver functions are not available"))
+        if store == None or store == 'data':
+            stream = self.data
+            for tr in stream:
+                tr = _add_stats(tr, store)
+
+        if store == 'rf':
+            if not hasattr(self, 'rf'):
+                raise(Exception("Warning: Receiver functions are not available"))
+            else:
+                stream = self.rf
+                for tr in stream:
+                    tr = _add_stats(tr, store)
 
-        stream = self.rf
-        for tr in stream:
-            tr = _add_rfstats(tr)
+        elif store == 'specs':
+            if not hasattr(self, 'specs'):
+                raise(Exception("Warning: Spectra are not available"))
+            else:
+                stream = self.specs
+                for tr in stream:
+                    tr = _add_stats(tr, store)
 
         return stream
 
@@ -949,7 +1239,7 @@ def save(self, file):
 
         """
 
-        import pickle
         output = open(file, 'wb')
         pickle.dump(self, output)
         output.close()
+
diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py
index 362588f..20aedc4 100644
--- a/rfpy/scripts/rfpy_calc.py
+++ b/rfpy/scripts/rfpy_calc.py
@@ -24,18 +24,13 @@
 
 
 # Import modules and functions
+from argparse import ArgumentParser
+from os.path import exists as exist
 import numpy as np
 import pickle
 import stdb
-from obspy.clients.fdsn import Client
-from obspy import Catalog, UTCDateTime
-from http.client import IncompleteRead
-from rfpy import utils, RFData
+from rfpy import RFData
 from pathlib import Path
-from argparse import ArgumentParser
-from os.path import exists as exist
-from obspy import UTCDateTime
-from numpy import nan
 
 
 def get_calc_arguments(argv=None):
@@ -48,13 +43,12 @@ def get_calc_arguments(argv=None):
 
     parser = ArgumentParser(
         usage="%(prog)s [arguments] <station database>",
-        description="Script used to download and pre-process " +
-        "three-component ('Z', 'N', and 'E'), seismograms for individual " +
-        "events and calculate teleseismic P-wave receiver functions" +
-        "This version requests data on the fly for a given date " +
-        "range. Data are requested from the internet using the " +
-        "client services framework. The stations are processed one " +
-        "by one and the data are stored to disk.")
+        description="Script used to re-calculate receiver functions " +
+        "that already exist on disk, but using different " +
+        "processing options. The stations are processed one " +
+        "by one and the data are stored to disk. " +
+        " \n" +
+        "Note: The sampling rate cannot be changed to a new rate")
 
     # General Settings
     parser.add_argument(
@@ -79,13 +73,6 @@ def get_calc_arguments(argv=None):
         dest="verb",
         default=False,
         help="Specify to increase verbosity.")
-    parser.add_argument(
-        "-O", "--overwrite",
-        action="store_true",
-        dest="ovr",
-        default=False,
-        help="Force the overwriting of pre-existing data. " +
-        "[Default False]")
     parser.add_argument(
         "-L", "--long-name",
         action="store_true",
@@ -96,183 +83,34 @@ def get_calc_arguments(argv=None):
         "names, regardless of the key type of the database."
     )
 
-    # Server Settings
-    ServerGroup = parser.add_argument_group(
-        title="Server Settings",
-        description="Settings associated with which "
-        "datacenter to log into.")
-    ServerGroup.add_argument(
-        "-S", "--Server",
-        action="store",
-        type=str,
-        dest="Server",
-        default="IRIS",
-        help="Specify the server to connect to. Options include: " +
-        "BGR, ETH, GEONET, GFZ, INGV, IPGP, IRIS, KOERI, " +
-        "LMU, NCEDC, NEIP, NERIES, ODC, ORFEUS, RESIF, SCEDC, USGS, USP. " +
-        "[Default IRIS]")
-    ServerGroup.add_argument(
-        "-U", "--User-Auth",
-        action="store",
-        type=str,
-        dest="UserAuth",
-        default="",
-        help="Enter your IRIS Authentification Username and Password " +
-        "(--User-Auth='username:authpassword') to " +
-        "access and download restricted data. " +
-        "[Default no user and password]")
-
-    # Database Settings
-    DataGroup = parser.add_argument_group(
-        title="Local Data Settings",
-        description="Settings associated with defining " +
-        "and using a local data base of pre-downloaded " +
-        "day-long SAC or MSEED files.")
-    DataGroup.add_argument(
-        "--local-data",
-        action="store",
-        type=str,
-        dest="localdata",
-        default=None,
-        help="Specify a comma separated list of paths containing " +
-        "day-long sac or mseed files of data already downloaded. " +
-        "If data exists for a seismogram is already present on disk, " +
-        "it is selected preferentially over downloading " +
-        "the data using the Client interface")
-    DataGroup.add_argument(
-    	"--dtype",
-    	action="store",
-    	type=str,
-    	dest="dtype",
-    	default='SAC',
-    	help="Specify the data archive file type, either SAC " +
-    	" or MSEED. Note the default behaviour is to search for " +
-    	"SAC files. Local archive files must have extensions of '.SAC' "+
-    	" or '.MSEED. These are case dependent, so specify the correct case"+
-        "here.")
-    DataGroup.add_argument(
-        "--no-data-zero",
-        action="store_true",
-        dest="ndval",
-        default=False,
-        help="Specify to force missing data to be set as zero, rather " +
-        "than default behaviour which sets to nan. Note this is applied " +
-        "only to the SAC data")
-    DataGroup.add_argument(
-        "--no-local-net",
-        action="store_false",
-        dest="useNet",
-        default=True,
-        help="Specify to prevent using the Network code in the " +
-        "search for local data (sometimes for CN stations " +
-        "the dictionary name for a station may disagree with that " +
-        "in the filename. [Default Network used]")
-    DataGroup.add_argument(
-        "--save-Z12",
-        action="store_true",
-        dest="saveZ12",
-        default=False,
-        help="Specify to save Z12 (un-rotated) components. [Default " +
-        "False]")
-
-    # Event Selection Criteria
-    EventGroup = parser.add_argument_group(
-        title="Event Settings",
-        description="Settings associated with refining " +
-        "the events to include in matching event-station pairs")
-    EventGroup.add_argument(
-        "--start",
-        action="store",
-        type=str,
-        dest="startT",
-        default="",
-        help="Specify a UTCDateTime compatible string representing " +
-        "the start time for the event search. This will override any " +
-        "station start times. [Default start date of station]")
-    EventGroup.add_argument(
-        "--end",
-        action="store",
-        type=str,
-        dest="endT",
-        default="",
-        help="Specify a UTCDateTime compatible string representing " +
-        "the end time for the event search. This will override any " +
-        "station end times [Default end date of station]")
-    EventGroup.add_argument(
-        "--reverse", "-R",
+    # Constants Settings
+    ConstGroup = parser.add_argument_group(
+        title='Parameter Settings',
+        description="Miscellaneous default values and settings")
+    ConstGroup.add_argument(
+        "--Z12",
         action="store_true",
-        dest="reverse",
+        dest="Z12",
         default=False,
-        help="Reverse order of events. Default behaviour starts at " +
-        "oldest event and works towards most recent. Specify reverse " +
-        "order and instead the program will start with the most recent " +
-        "events and work towards older")
-    EventGroup.add_argument(
-        "--minmag",
-        action="store",
-        type=float,
-        dest="minmag",
-        default=6.0,
-        help="Specify the minimum magnitude of event for which to search. " +
-        "[Default 6.0]")
-    EventGroup.add_argument(
-        "--maxmag",
-        action="store",
-        type=float,
-        dest="maxmag",
-        default=9.0,
-        help="Specify the maximum magnitude of event for which to search. " +
-        "[Default None, i.e. no limit]")
-
-    # Geometry Settings
-    PhaseGroup = parser.add_argument_group(
-        title="Geometry Settings",
-        description="Settings associatd with the "
-        "event-station geometries for the specified phase")
-    PhaseGroup.add_argument(
+        help="Use Z12 data if available. [Default uses ZNE data]")
+    ConstGroup.add_argument(
         "--phase",
         action="store",
         type=str,
         dest="phase",
-        default='P',
+        default='allP',
         help="Specify the phase name to use. Be careful with the distance. " +
-        "setting. Options are 'P' or 'PP'. [Default 'P']")
-    PhaseGroup.add_argument(
-        "--mindist",
-        action="store",
-        type=float,
-        dest="mindist",
-        default=None,
-        help="Specify the minimum great circle distance (degrees) between " +
-        "the station and event. [Default depends on phase]")
-    PhaseGroup.add_argument(
-        "--maxdist",
-        action="store",
-        type=float,
-        dest="maxdist",
-        default=None,
-        help="Specify the maximum great circle distance (degrees) between " +
-        "the station and event. [Default depends on phase]")
-
-    # Constants Settings
-    ConstGroup = parser.add_argument_group(
-        title='Parameter Settings',
-        description="Miscellaneous default values and settings")
-    ConstGroup.add_argument(
-        "--sampling-rate",
-        action="store",
-        type=float,
-        dest="new_sampling_rate",
-        default=10.,
-        help="Specify new sampling rate in Hz. [Default 10.]")
+        "setting. Options are 'P', 'PP', 'allP', 'S', 'SKS' or 'allS'. " +
+        "[Default 'allP']")
     ConstGroup.add_argument(
-        "--dts",
+        "--resample",
         action="store",
         type=float,
-        dest="dts",
-        default=150.,
-        help="Specify the window length in sec (symmetric about arrival " +
-        "time). [Default 150.]")
+        dest="resample",
+        default=None,
+        help="Specify the new sampling-rate for the receiver functions. " +
+        "Note the sampling rate of the original data (ZNE or Z12) stored " +
+        "on disk is unchanged. [Default None]")
     ConstGroup.add_argument(
         "--align",
         action="store",
@@ -319,7 +157,7 @@ def get_calc_arguments(argv=None):
         type=float,
         dest="fmin",
         default=0.05,
-        help="Specify the minimum frequency corner for SNR and CC " +
+        help="Specify the minimum frequency corner for SNR " +
         "filter (Hz). [Default 0.05]")
     ConstGroup.add_argument(
         "--fmax",
@@ -327,10 +165,10 @@ def get_calc_arguments(argv=None):
         type=float,
         dest="fmax",
         default=1.0,
-        help="Specify the maximum frequency corner for SNR and CC " +
+        help="Specify the maximum frequency corner for SNR " +
         "filter (Hz). [Default 1.0]")
 
-    # Constants Settings
+    # Deconvolution Settings
     DeconGroup = parser.add_argument_group(
         title='Deconvolution Settings',
         description="Parameters for deconvolution")
@@ -369,109 +207,16 @@ def get_calc_arguments(argv=None):
     if len(args.stkeys) > 0:
         args.stkeys = args.stkeys.split(',')
 
-    # construct start time
-    if len(args.startT) > 0:
-        try:
-            args.startT = UTCDateTime(args.startT)
-        except:
-            parser.error(
-                "Cannot construct UTCDateTime from start time: " +
-                args.startT)
-    else:
-        args.startT = None
-
-    # construct end time
-    if len(args.endT) > 0:
-        try:
-            args.endT = UTCDateTime(args.endT)
-        except:
-            parser.error(
-                "Cannot construct UTCDateTime from end time: " +
-                args.endT)
-    else:
-        args.endT = None
-
-    # Parse User Authentification
-    if not len(args.UserAuth) == 0:
-        tt = args.UserAuth.split(':')
-        if not len(tt) == 2:
-            parser.error(
-                "Error: Incorrect Username and Password " +
-                "Strings for User Authentification")
-        else:
-            args.UserAuth = tt
-    else:
-        args.UserAuth = []
-
-    # Parse Local Data directories
-    if args.localdata is not None:
-        args.localdata = args.localdata.split(',')
-    else:
-        args.localdata = []
-
-    # Check Datatype specification
-    if (not args.dtype.upper() == 'MSEED') and  (not args.dtype.upper() == 'SAC'):
-    	parser.error(
-    		"Error: Local Data Archive must be of types 'SAC'" +
-    		"or MSEED. These must match the file extensions for " +
-    		" the archived data.")
-
-    # Check NoData Value
-    if args.ndval:
-        args.ndval = 0.0
-    else:
-        args.ndval = nan
-
-    # Check distances for selected phase
-    if args.phase not in ['P', 'PP', 'S', 'SKS']:
+    if args.phase not in ['P', 'PP', 'allP', 'S', 'SKS', 'allS']:
         parser.error(
-            "Error: choose between 'P', 'PP', 'S' and 'SKS'.")
-    if args.phase == 'P':
-        if not args.mindist:
-            args.mindist = 30.
-        if not args.maxdist:
-            args.maxdist = 100.
-        if args.mindist < 30. or args.maxdist > 100.:
-            parser.error(
-                "Distances should be between 30 and 100 deg. for " +
-                "teleseismic 'P' waves.")
-    elif args.phase == 'PP':
-        if not args.mindist:
-            args.mindist = 100.
-        if not args.maxdist:
-            args.maxdist = 180.
-        if args.mindist < 100. or args.maxdist > 180.:
-            parser.error(
-                "Distances should be between 100 and 180 deg. for " +
-                "teleseismic 'PP' waves.")
-    elif args.phase == 'S':
-        if not args.mindist:
-            args.mindist = 55.
-        if not args.maxdist:
-            args.maxdist = 85.
-        if args.mindist < 55. or args.maxdist > 85.:
-            parser.error(
-                "Distances should be between 55 and 85 deg. for " +
-                "teleseismic 'S' waves.")
-    elif args.phase == 'SKS':
-        if not args.mindist:
-            args.mindist = 85.
-        if not args.maxdist:
-            args.maxdist = 115.
-        if args.mindist < 85. or args.maxdist > 115.:
-            parser.error(
-                "Distances should be between 85 and 115 deg. for " +
-                "teleseismic 'SKS' waves.")
-
-    if args.pre_filt is not None:
-        args.pre_filt = [float(val) for val in args.pre_filt.split(',')]
-        args.pre_filt = sorted(args.pre_filt)
-        if (len(args.pre_filt)) != 2:
-            parser.error(
-                "Error: --pre-filt should contain 2 " +
-                "comma-separated floats")
+            "Error: choose between 'P', 'PP', 'allP', 'S', 'SKS' and 'allS'.")
+    if args.phase == 'allP':
+        args.listphase = ['P', 'PP']
+    elif args.phase == 'allS':
+        args.listphase = ['S', 'SKS']
+    else:
+        args.listphase = [args.phase]
 
-    # Check alignment arguments
     if args.align is None:
         args.align = 'ZRT'
     elif args.align not in ['ZRT', 'LQT', 'PVH']:
@@ -479,16 +224,19 @@ def get_calc_arguments(argv=None):
             "Error: Incorrect alignment specifier. Should be " +
             "either 'ZRT', 'LQT', or 'PVH'.")
 
-    if args.dt_snr > args.dts:
-        args.dt_snr = args.dts - 10.
-        print("SNR window > data window. Defaulting to data " +
-              "window minus 10 sec.")
-
     if args.method not in ['wiener', 'water', 'multitaper']:
         parser.error(
             "Error: 'method' should be either 'wiener', 'water' or " +
             "'multitaper'")
 
+    if args.pre_filt is not None:
+        args.pre_filt = [float(val) for val in args.pre_filt.split(',')]
+        args.pre_filt = sorted(args.pre_filt)
+        if (len(args.pre_filt)) != 2:
+            parser.error(
+                "Error: --pre-filt should contain 2 " +
+                "comma-separated floats")
+
     return args
 
 
@@ -496,6 +244,7 @@ def main():
 
     print()
     print("##############################################")
+    print("#                                            #")
     print("#        __                          _       #")
     print("#  _ __ / _|_ __  _   _     ___ __ _| | ___  #")
     print("# | '__| |_| '_ \| | | |   / __/ _` | |/ __| #")
@@ -512,6 +261,9 @@ def main():
     # Load Database
     db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys)
 
+    # Track processed folders
+    procfold = []
+
     # Loop over station keys
     for stkey in list(stkeys):
 
@@ -524,36 +276,12 @@ def main():
             stfld = stkey.split('.')[0]+"."+stkey.split('.')[1]
 
         # Define path to see if it exists
-        if args.phase in ['P', 'PP']:
+        if args.phase in ['P', 'PP', 'allP']:
             datapath = Path('P_DATA') / stfld
-        elif args.phase in ['S', 'SKS']:
+        elif args.phase in ['S', 'SKS', 'allS']:
             datapath = Path('S_DATA') / stfld
-        if not datapath.exists():
-            print('Path to '+str(datapath)+' doesn`t exist - creating it')
-            datapath.mkdir(parents=True)
-
-        # Establish client for data
-        if len(args.UserAuth) == 0:
-            data_client = Client(args.Server)
-        else:
-            data_client = Client(args.Server, user=args.UserAuth[0],
-                                 password=args.UserAuth[1])
-
-        # Establish client for events
-        event_client = Client()
-
-        # Get catalogue search start time
-        if args.startT is None:
-            tstart = sta.startdate
-        else:
-            tstart = args.startT
-
-        # Get catalogue search end time
-        if args.endT is None:
-            tend = sta.enddate
-        else:
-            tend = args.endT
-        if tstart > sta.enddate or tend < sta.startdate:
+        if not datapath.is_dir():
+            print('Path to ' + str(datapath) + ' doesn`t exist - continuing')
             continue
 
         # Temporary print locations
@@ -569,238 +297,117 @@ def main():
         print(" ")
         print(" ")
         print("|===============================================|")
-        print("|===============================================|")
         print("|                   {0:>8s}                    |".format(
             sta.station))
         print("|===============================================|")
-        print("|===============================================|")
         print("|  Station: {0:>2s}.{1:5s}                            |".format(
             sta.network, sta.station))
         print("|      Channel: {0:2s}; Locations: {1:15s}  |".format(
             sta.channel, ",".join(tlocs)))
         print("|      Lon: {0:7.2f}; Lat: {1:6.2f}                |".format(
             sta.longitude, sta.latitude))
-        print("|      Start time: {0:19s}          |".format(
-            sta.startdate.strftime("%Y-%m-%d %H:%M:%S")))
-        print("|      End time:   {0:19s}          |".format(
-            sta.enddate.strftime("%Y-%m-%d %H:%M:%S")))
         print("|-----------------------------------------------|")
-        print("| Searching Possible events:                    |")
-        print("|   Start: {0:19s}                  |".format(
-            tstart.strftime("%Y-%m-%d %H:%M:%S")))
-        print("|   End:   {0:19s}                  |".format(
-            tend.strftime("%Y-%m-%d %H:%M:%S")))
-        if args.maxmag is None:
-            print("|   Mag:   >{0:3.1f}", format(args.minmag) +
-                  "                                 |")
-        else:
-            msg = "|   Mag:   {0:3.1f}".format(args.minmag) + \
-                " - {0:3.1f}".format(args.maxmag) + \
-                "                            |"
-            print(msg)
-
-        print("| ...                                           |")
-
-        # Get catalogue using deployment start and end
-        try:
-            cat = event_client.get_events(
-                starttime=tstart, endtime=tend,
-                minmagnitude=args.minmag, maxmagnitude=args.maxmag)
-        except IncompleteRead:
-            # See http.client.IncompleteRead
-            # https://github.com/obspy/obspy/issues/3028#issue-1179808237
-            
-            # Read yearly chunks
-            print('| Warning: Unable to get entire catalog at once |')
-            print('| Trying to get one year at a time              |')
-            print('|                                               |')
-            
-            chunk = 365 * 86400  # Query length in seconds
-            cat = Catalog()
-            tend = min(tend, UTCDateTime.now())
-            while tstart < tend:
-                print("| Start:   {0:19s}                  |".format(
-                    tstart.strftime("%Y-%m-%d %H:%M:%S")))
-                cat += event_client.get_events(
-                    starttime=tstart, endtime=tstart + chunk,
-                    minmagnitude=args.minmag, maxmagnitude=args.maxmag)
-
-                # Make sure that we go all the way to tend
-                if tstart + chunk > tend:
-                    chunk = tend - tstart
-
-                    # But do not get caught up in an infinite loop due to
-                    # rounding errors
-                    if chunk <= 1:
-                        break
-
-                tstart += chunk
-
-        # Total number of events in Catalogue
-        nevK = 0
-        nevtT = len(cat)
-        print(
-            "|  Found {0:5d}".format(nevtT) +
-            " possible events                  |")
-        ievs = range(0, nevtT)
-
-        # Get Local Data Availabilty
-        if len(args.localdata) > 0:
-            print("|-----------------------------------------------|")
-            print("| Cataloging Local Data...                      |")
-            if args.useNet:
-                stalcllist = utils.list_local_data_stn(
-                    lcldrs=args.localdata, sta=sta.station,
-                    net=sta.network, dtype=args.dtype, altnet=sta.altnet)
-                print("|   {0:>2s}.{1:5s}: {2:6d}".format(
-                    sta.network, sta.station, len(stalcllist)) +
-                    " files                      |")
-                #print(stalcllist[0:10])
-            else:
-                stalcllist = utils.list_local_data_stn(
-                    lcldrs=args.localdata, sta=sta.station, dtype=args.dtype)
-                print("|   {0:5s}: {1:6d} files                " +
-                      "        |".format(
-                          sta.station, len(stalcllist)))
-        else:
-            stalcllist = []
-        print("|===============================================|")
 
-        # Select order of processing
-        if args.reverse:
-            ievs = range(0, nevtT)
-        else:
-            ievs = range(nevtT-1, -1, -1)
+        # Check for folder already processed
+        if stfld in procfold:
+            print('  {0} already processed...skipping   '.format(stfld))
+            continue
 
-        # Read through catalogue
-        for iev in ievs:
+        datafiles = [x for x in datapath.iterdir() if x.is_dir()]
+        for folder in datafiles:
 
-            # Extract event
-            ev = cat[iev]
+            # Skip hidden folders
+            if folder.name.startswith('.'):
+                continue
 
-            # Initialize RF object with station info
+            # Re-initialize RFData object
             rfdata = RFData(sta)
 
-            # Add event to rfdata object
-            accept = rfdata.add_event(
-                ev, gacmin=args.mindist, gacmax=args.maxdist,
-                phase=args.phase, returned=True)
-
-            # Define time stamp
-            yr = str(rfdata.meta.time.year).zfill(4)
-            jd = str(rfdata.meta.time.julday).zfill(3)
-            hr = str(rfdata.meta.time.hour).zfill(2)
-
-            # If event is accepted (data exists)
-            if accept:
-
-                # Display Event Info
-                nevK = nevK + 1
-                if args.reverse:
-                    inum = iev + 1
-                else:
-                    inum = nevtT - iev + 1
-                print(" ")
-                print("*"*50)
-                print("* #{0:d} ({1:d}/{2:d}):  {3:13s} {4}".format(
-                    nevK, inum, nevtT, rfdata.meta.time.strftime(
-                        "%Y%m%d_%H%M%S"), stkey))
-                if args.verb:
-                    print("*   Phase: {}".format(args.phase))
-                    print("*   Origin Time: " +
-                          rfdata.meta.time.strftime("%Y-%m-%d %H:%M:%S"))
-                    print(
-                        "*   Lat: {0:6.2f};        Lon: {1:7.2f}".format(
-                            rfdata.meta.lat, rfdata.meta.lon))
-                    print(
-                        "*   Dep: {0:6.2f} km;     Mag: {1:3.1f}".format(
-                            rfdata.meta.dep, rfdata.meta.mag))
-                    print(
-                        "*   Dist: {0:7.2f} km;".format(rfdata.meta.epi_dist) +
-                        "   Epi dist: {0:6.2f} deg\n".format(rfdata.meta.gac) +
-                        "*   Baz:  {0:6.2f} deg;".format(rfdata.meta.baz) +
-                        "   Az: {0:6.2f} deg".format(rfdata.meta.az))
-
-                # Event Folder
-                timekey = rfdata.meta.time.strftime("%Y%m%d_%H%M%S")
-                evtdir = datapath / timekey
-                RFfile = evtdir / 'RF_Data.pkl'
-                ZNEfile = evtdir / 'ZNE_Data.pkl'
-                metafile = evtdir / 'Meta_Data.pkl'
-                stafile = evtdir / 'Station_Data.pkl'
-
-                # Check if RF data already exist and overwrite has been set
-                if evtdir.exists():
-                    if RFfile.exists():
-                        if not args.ovr:
-                            continue
-
-                # Get data
-                has_data = rfdata.download_data(
-                    client=data_client, dts=args.dts, stdata=stalcllist,
-                    ndval=args.ndval, dtype=args.dtype, new_sr=args.new_sampling_rate,
-                    returned=True, verbose=args.verb)
-
-                if not has_data:
-                    continue
-
-                # Create Folder if it doesn't exist
-                if not evtdir.exists():
-                    evtdir.mkdir(parents=True)
-
-                # Save ZNE Traces
-                pickle.dump(rfdata.data, open(ZNEfile, "wb"))
-
-                # Save Z12 if components exist
-                if hasattr(rfdata, "dataZ12"):
-                    Z12file = evtdir / 'Z12_Data.pkl'
-                    pickle.dump(rfdata.dataZ12, open(Z12file, "wb"))
-
-                # Rotate from ZNE to 'align' ('ZRT', 'LQT', or 'PVH')
-                rfdata.rotate(vp=args.vp, vs=args.vs, align=args.align)
-
-                # Calculate snr over dt_snr seconds
-                rfdata.calc_snr(
-                    dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax)
-                if args.verb:
-                    print("* SNR: {}".format(rfdata.meta.snr))
-
-                # Make sure no processing happens for NaNs
-                if np.isnan(rfdata.meta.snr):
-                    if args.verb:
-                        print("* SNR NaN...Skipping")
-                    print("*"*50)
-                    continue
-
-                # Deconvolve data
-                rfdata.deconvolve(
+            # Load meta data
+            metafile = folder / "Meta_Data.pkl"
+            if not metafile.is_file():
+                continue
+            rfdata.meta = pickle.load(open(metafile, 'rb'))
+
+            # Skip data not in list of phases
+            if rfdata.meta.phase not in args.listphase:
+                continue
+
+            if args.verb:
+                print("* Station: {0}; folder: {1}".format(stkey, folder))
+
+            if args.Z12:
+                try:
+                    # Try loading Z12_Data
+                    Z12file = folder / "Z12_Data.pkl"
+                    rfdata.data = pickle.load(open(Z12file, 'rb'))
+                    # Remove rotated flag and snr flag
+                    rfdata.meta.rotated = False
+                    rfdata.rotate(align='ZNE')
+                except:
+
+                    print("Z12_Data.pkl not available - using ZNE_Data.pkl")
+                    # Load ZNE data
+                    ZNEfile = folder / "ZNE_Data.pkl"
+                    rfdata.data = pickle.load(open(ZNEfile, 'rb'))
+            else:
+                # Load ZNE data
+                ZNEfile = folder / "ZNE_Data.pkl"
+                rfdata.data = pickle.load(open(ZNEfile, 'rb'))
+
+            # Resample if requested
+            if args.resample:
+                rfdata.data.resample(
+                    args.resample, no_filter=False)
+
+            # Remove rotated flag and snr flag
+            rfdata.meta.rotated = False
+            rfdata.meta.snr = None
+
+            # Rotate from ZNE to 'align' ('ZRT', 'LQT', or 'PVH')
+            rfdata.rotate(vp=args.vp, vs=args.vs, align=args.align)
+
+            # Calculate SNR
+            rfdata.calc_snr(dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax)
+
+            if args.verb:
+                print("* SNR: {}".format(rfdata.meta.snr))
+
+            rfdata.calc_spectra(
                     vp=args.vp, vs=args.vs,
                     align=args.align, method=args.method,
-                    gfilt=args.gfilt, wlevel=args.wlevel,
-                    pre_filt=args.pre_filt)
+                    wavelet='complete', envelope_threshold=0.05, 
+                    time=5)
+
+            # Deconvolve data
+            rfdata.deconvolve(
+                    align=args.align, method=args.method,
+                    gfilt=args.gfilt, wlevel=args.wlevel)
+
+            # Get cross-correlation QC
+            rfdata.calc_cc()
 
-                # Get cross-correlation QC
-                rfdata.calc_cc()
-                if args.verb:
-                    print("* CC: {}".format(rfdata.meta.cc))
+            if args.verb:
+                print("* CC: {}".format(rfdata.meta.cc))
 
-                # Convert to Stream
-                rfstream = rfdata.to_stream()
+            # Convert to Stream
+            rfstream = rfdata.to_stream(store='rf')
 
-                # Save event meta data
-                pickle.dump(rfdata.meta, open(metafile, "wb"))
+            # Save RF Traces
+            RFfile = folder / "RF_Data.pkl"
+            pickle.dump(rfstream, open(RFfile, "wb"))
 
-                # Save Station Data
-                pickle.dump(rfdata.sta, open(stafile, "wb"))
+            # Save Meta data
+            metafile = folder / "Meta_Data.pkl"
+            pickle.dump(rfdata.meta, open(metafile, 'wb'))
 
-                # Save RF Traces
-                pickle.dump(rfstream, open(RFfile, "wb"))
+            # Update
+            if args.verb:
+                print("* Output files written")
+                print("**************************************************")
 
-                # Update
-                if args.verb:
-                    print("* Wrote Output Files to: ")
-                    print("*     "+str(evtdir))
-                print("*"*50)
+            # Update processed folders
+            procfold.append(stfld)
 
 
 if __name__ == "__main__":
diff --git a/rfpy/scripts/rfpy_recalc.py b/rfpy/scripts/rfpy_calc_simdec.py
similarity index 58%
rename from rfpy/scripts/rfpy_recalc.py
rename to rfpy/scripts/rfpy_calc_simdec.py
index 74cb57a..0d7a494 100644
--- a/rfpy/scripts/rfpy_recalc.py
+++ b/rfpy/scripts/rfpy_calc_simdec.py
@@ -29,11 +29,12 @@
 import numpy as np
 import pickle
 import stdb
-from rfpy import RFData
+from obspy import Stream
+from rfpy import RFData, binning, plotting, utils
 from pathlib import Path
 
 
-def get_recalc_arguments(argv=None):
+def get_simdec_arguments(argv=None):
     """
     Get Options from :class:`~optparse.OptionParser` objects.
 
@@ -43,12 +44,9 @@ def get_recalc_arguments(argv=None):
 
     parser = ArgumentParser(
         usage="%(prog)s [arguments] <station database>",
-        description="Script used to re-calculate receiver functions " +
-        "that already exist on disk, but using different " +
-        "processing options. The stations are processed one " +
-        "by one and the data are stored to disk. " +
-        " \n" +
-        "Note: The sampling rate cannot be changed to a new rate")
+        description="Script used to calculate receiver functions " +
+        "using simultaneous deconvolution. The stations are processed one " +
+        "by one and the data are stored to disk. ")
 
     # General Settings
     parser.add_argument(
@@ -196,6 +194,100 @@ def get_recalc_arguments(argv=None):
         default=0.01,
         help="Specify the water level, used in the 'water' method. " +
         "[Default 0.01]")
+    DeconGroup.add_argument(
+        "--try-auto-wlevel",
+        action="store_true",
+        dest="auto_wlevel",
+        help="Try to determine waterlevel using Generalized Cross Validation.")
+
+    PreGroup = parser.add_argument_group(
+        title='Pre-processing Settings',
+        description="Options for pre-processing " +
+        "data before stacking")
+    PreGroup.add_argument(
+        "--snr",
+        action="store",
+        type=float,
+        dest="snr",
+        default=-9999.,
+        help="Specify the vertical component SNR threshold for " +
+        "extracting receiver functions. [Default 5.]")
+    PreGroup.add_argument(
+        "--snrh",
+        action="store",
+        type=float,
+        dest="snrh",
+        default=-9999.,
+        help="Specify the horizontal component SNR threshold for " +
+        "extracting receiver functions. [Default None]")
+    PreGroup.add_argument(
+        "--nbaz",
+        action="store",
+        dest="nbaz",
+        type=int,
+        default=36,
+        help="Specify integer number of back-azimuth bins to consider " +
+        "(typically 36 or 72). If not None, the plot will show receiver " +
+        "functions sorted by back-azimuth values. [Default 36]")
+    PreGroup.add_argument(
+        "--nslow",
+        action="store",
+        dest="nslow",
+        type=int,
+        default=20,
+        help="Specify integer number of slowness bins to consider " +
+        "(range of slowness values is 0.04 to 0.08). " + 
+        "If not None, the plot will show receiver " +
+        "functions sorted by slowness values. [Default 20]")
+    PreGroup.add_argument(
+        "--bp",
+        action="store",
+        type=str,
+        dest="bp",
+        default=None,
+        help="Specify the corner frequencies for the bandpass filter. " +
+        "[Default no filtering]")
+
+    PlotGroup = parser.add_argument_group(
+        title='Plot Settings',
+        description="Options for plot format")
+    PlotGroup.add_argument(
+        "--normalize",
+        action="store_true",
+        dest="norm",
+        default=False,
+        help="Set this option to produce receiver functions normalized " +
+        "by the max amplitude of stacked RFs. [Default False]")
+    PlotGroup.add_argument(
+        "--trange",
+        action="store",
+        default=None,
+        type=str,
+        dest="trange",
+        help="Specify the time range for the x-axis (sec). Negative times " +
+        "are allowed [Default 0., 30.]")
+    PlotGroup.add_argument(
+        "--save",
+        action="store_true",
+        dest="saveplot",
+        default=False,
+        help="Set this option if you wish to save the figure. [Default " +
+        "does not save figure]")
+    PlotGroup.add_argument(
+        "--title",
+        action="store",
+        dest="titleplot",
+        type=str,
+        default='',
+        help="Specify title of figure. [Default None]")
+    PlotGroup.add_argument(
+        "--format",
+        action="store",
+        type=str,
+        dest="form",
+        default="png",
+        help="Specify format of figure. Can be any one of the valid" +
+        "matplotlib formats: 'png', 'jpg', 'eps', 'pdf'. [Default 'png']")
 
     args = parser.parse_args(argv)
 
@@ -229,6 +321,24 @@ def get_recalc_arguments(argv=None):
             "Error: 'method' should be either 'wiener', 'water' or " +
             "'multitaper'")
 
+    if args.bp is not None:
+        args.bp = [float(val) for val in args.bp.split(',')]
+        args.bp = sorted(args.bp)
+        if (len(args.bp)) != 2:
+            parser.error(
+                "Error: --bp should contain 2 " +
+                "comma-separated floats")
+
+    if args.trange is None:
+        args.trange = [0., 30.]
+    else:
+        args.trange = [float(val) for val in args.trange.split(',')]
+        args.trange = sorted(args.trange)
+        if (len(args.trange)) != 2:
+            parser.error(
+                "Error: --trange should contain 2 " +
+                "comma-separated floats")
+
     if args.pre_filt is not None:
         args.pre_filt = [float(val) for val in args.pre_filt.split(',')]
         args.pre_filt = sorted(args.pre_filt)
@@ -243,20 +353,20 @@ def get_recalc_arguments(argv=None):
 def main():
 
     print()
-    print("########################################################")
-    print("#                                                      #")
-    print("#        __                                    _       #")
-    print("#  _ __ / _|_ __  _   _     _ __ ___  ___ __ _| | ___  #")
-    print("# | '__| |_| '_ \| | | |   | '__/ _ \/ __/ _` | |/ __| #")
-    print("# | |  |  _| |_) | |_| |   | | |  __/ (_| (_| | | (__  #")
-    print("# |_|  |_| | .__/ \__, |___|_|  \___|\___\__,_|_|\___| #")
-    print("#          |_|    |___/_____|                          #")
-    print("#                                                      #")
-    print("########################################################")
+    print("##################################################################################")
+    print("#                                                                                #")
+    print("#        __                          _              _               _            #")
+    print("#  _ __ / _|_ __  _   _     ___ __ _| | ___     ___(_)_ __ ___   __| | ___  ___  #")
+    print("# | '__| |_| '_ \| | | |   / __/ _` | |/ __|   / __| | '_ ` _ \ / _` |/ _ \/ __| #")
+    print("# | |  |  _| |_) | |_| |  | (_| (_| | | (__    \__ \ | | | | | | (_| |  __/ (__  #")
+    print("# |_|  |_| | .__/ \__, |___\___\__,_|_|\___|___|___/_|_| |_| |_|\__,_|\___|\___| #")
+    print("#          |_|    |___/_____|             |_____|                                #")
+    print("#                                                                                #")
+    print("##################################################################################")
     print()
 
     # Run Input Parser
-    args = get_recalc_arguments()
+    args = get_simdec_arguments()
 
     # Load Database
     db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys)
@@ -308,10 +418,16 @@ def main():
             sta.longitude, sta.latitude))
         print("|-----------------------------------------------|")
 
-        # Check for folder already processed
-        if stfld in procfold:
-            print('  {0} already processed...skipping   '.format(stfld))
-            continue
+        # # Check for folder already processed
+        # if stfld in procfold:
+        #     print('  {0} already processed...skipping   '.format(stfld))
+        #     continue
+
+        #LLstreams = Stream()
+        #QLstreams = Stream()
+        #TLstreams = Stream()
+        #NNstreams = Stream()
+        RFs = np.array([])
 
         datafiles = [x for x in datapath.iterdir() if x.is_dir()]
         for folder in datafiles:
@@ -370,40 +486,93 @@ def main():
             # Calculate SNR
             rfdata.calc_snr(dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax)
 
-            if args.verb:
-                print("* SNR: {}".format(rfdata.meta.snr))
-
-            # Deconvolve data
-            rfdata.deconvolve(
-                vp=args.vp, vs=args.vs,
-                align=args.align, method=args.method,
-                gfilt=args.gfilt, wlevel=args.wlevel,
-                pre_filt=args.pre_filt)
-
-            # Get cross-correlation QC
-            rfdata.calc_cc()
-
-            if args.verb:
-                print("* CC: {}".format(rfdata.meta.cc))
+            # QC Thresholding
+            if rfdata.meta.snrh < args.snrh:
+                continue
+            if rfdata.meta.snr < args.snr:
+                continue
 
-            # Convert to Stream
-            rfstream = rfdata.to_stream()
+            # Calculate spectra
+            rfdata.calc_spectra(
+                    vp=args.vp, vs=args.vs,
+                    align=args.align, method=args.method,
+                    wavelet='complete', envelope_threshold=0.05, 
+                    time=5, norm=True)
+
+            ## Convert to Stream
+            #rfstream = rfdata.to_stream(store='specs')
+
+            ## Add to Stream objects
+            #LLstreams.append(rfstream[0])
+            #QLstreams.append(rfstream[1])
+            #TLstreams.append(rfstream[2])
+            #NNstreams.append(rfstream[3])
+            RFs = np.append(RFs, RF)
+
+        # Bin into baz and slowness bins
+        #LL_tmp, = binning.bin_baz_slow(LLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1,
+            #pws=False, include_empty=False)
+        #QL_tmp, = binning.bin_baz_slow(QLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1,
+            #pws=False, include_empty=False)
+        #TL_tmp, = binning.bin_baz_slow(TLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1,
+            #pws=False, include_empty=False)
+        #NN_tmp, = binning.bin_baz_slow(NNstreams, nbaz=args.nbaz+1, nslow=args.nslow+1,
+            #pws=False, include_empty=False)
+        iibin, baz_slows = binning.baz_slow_bin_indices(RFs, nbaz, nslow)
+
+        # Initialize new empty streams
+        RFL = Stream()
+        RFQ = Stream()
+        RFT = Stream()
+
+        # Cycle through binned spectra
+        #for i in range(len(LL_tmp)):
+        for n, ibin in enumerate(iibin):
+
+            # Receiver functions in bin
+            RFset = RFs[ibin]
+
+            # Initialize new RFData object with empty Meta header
+            rfdata = RFData(sta)
 
-            # Save RF Traces
-            RFfile = folder / "RF_Data.pkl"
-            pickle.dump(rfstream, open(RFfile, "wb"))
+            # Extract the specs from the binned streams
+            #rfdata.specs = Stream(traces=[LL_tmp[i], QL_tmp[i], TL_tmp[i], NN_tmp[i]])
+            rfdata.specs = binning.stack(RFset, which='specs')
 
-            # Save Meta data
-            metafile = folder / "Meta_Data.pkl"
-            pickle.dump(rfdata.meta, open(metafile, 'wb'))
-
-            # Update
-            if args.verb:
-                print("* Output files written")
-                print("**************************************************")
+            # Set waterlevel
+            wlevel = args.wlevel
+            if args.auto_wlevel:
+                try:
+                    wlevel = utils.optimal_wlevel(RFset)
+                except ValueError:
+                    print('Water level GCV function had no minimum.')
+                    print('Falling back to user supplied value.')
+                    wlevel = args.wlevel
 
-            # Update processed folders
-            procfold.append(stfld)
+            # Deconvolve
+            rfdata.deconvolve(
+                align=args.align, method=args.method,
+                gfilt=args.gfilt, wlevel=wlevel)
+
+            # Append RFs to streams
+            RFL.append(rfdata.rf[0])
+            RFQ.append(rfdata.rf[1])
+            RFT.append(rfdata.rf[2])
+
+        # Filter
+        if args.bp:
+            RFQ.filter('bandpass', freqmin=args.bp[0],
+                             freqmax=args.bp[1], corners=2,
+                             zerophase=True)
+            RFT.filter('bandpass', freqmin=args.bp[0],
+                             freqmax=args.bp[1], corners=2,
+                             zerophase=True)
+
+        plotting.panel(RFQ, RFT, tmin=args.trange[0], tmax=args.trange[1],
+            normalize=True, save=args.saveplot,
+                title=args.titleplot, form=args.form)
+        # # Update processed folders
+        # procfold.append(stfld)
 
 
 if __name__ == "__main__":
diff --git a/rfpy/scripts/rfpy_get_data.py b/rfpy/scripts/rfpy_get_data.py
new file mode 100644
index 0000000..6e7a8f8
--- /dev/null
+++ b/rfpy/scripts/rfpy_get_data.py
@@ -0,0 +1,631 @@
+#!/usr/bin/env python
+
+# Copyright 2019 Pascal Audet
+#
+# This file is part of RfPy.
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+
+# Import modules and functions
+import numpy as np
+import pickle
+import stdb
+from obspy.clients.fdsn import Client
+from rfpy import utils, RFData
+from pathlib import Path
+from argparse import ArgumentParser
+from os.path import exists as exist
+from obspy import UTCDateTime
+from numpy import nan
+
+
+def get_data_arguments(argv=None):
+    """
+    Get Options from :class:`~optparse.OptionParser` objects.
+
+    This function is used for data processing on-the-fly (requires web connection)
+
+    """
+
+    parser = ArgumentParser(
+        usage="%(prog)s [arguments] <station database>",
+        description="Script used to download and pre-process " +
+        "three-component ('Z', 'N', and 'E'), seismograms for individual " +
+        "events and calculate teleseismic P-wave receiver functions" +
+        "This version requests data on the fly for a given date " +
+        "range. Data are requested from the internet using the " +
+        "client services framework. The stations are processed one " +
+        "by one and the data are stored to disk.")
+
+    # General Settings
+    parser.add_argument(
+        "indb",
+        help="Station Database to process from.",
+        type=str)
+    parser.add_argument(
+        "--keys",
+        action="store",
+        type=str,
+        dest="stkeys",
+        default="",
+        help="Specify a comma separated list of station keys for " +
+        "which to perform the analysis. These must be " +
+        "contained within the station database. Partial keys will " +
+        "be used to match against those in the dictionary. For " +
+        "instance, providing IU will match with all stations in " +
+        "the IU network [Default processes all stations in the database]")
+    parser.add_argument(
+        "-v", "-V", "--verbose",
+        action="store_true",
+        dest="verb",
+        default=False,
+        help="Specify to increase verbosity.")
+    parser.add_argument(
+        "-O", "--overwrite",
+        action="store_true",
+        dest="ovr",
+        default=False,
+        help="Force the overwriting of pre-existing data. " +
+        "[Default False]")
+    parser.add_argument(
+        "-L", "--long-name",
+        action="store_true",
+        dest="lkey",
+        default=False,
+        help="Force folder names to use long-key form (NET.STN.CHN). " +
+        "Default behaviour uses short key form (NET.STN) for the folder " +
+        "names, regardless of the key type of the database."
+    )
+
+    # Server Settings
+    ServerGroup = parser.add_argument_group(
+        title="Server Settings",
+        description="Settings associated with which "
+        "datacenter to log into.")
+    ServerGroup.add_argument(
+        "-S", "--Server",
+        action="store",
+        type=str,
+        dest="Server",
+        default="IRIS",
+        help="Specify the server to connect to. Options include: " +
+        "BGR, ETH, GEONET, GFZ, INGV, IPGP, IRIS, KOERI, " +
+        "LMU, NCEDC, NEIP, NERIES, ODC, ORFEUS, RESIF, SCEDC, USGS, USP. " +
+        "[Default IRIS]")
+    ServerGroup.add_argument(
+        "-U", "--User-Auth",
+        action="store",
+        type=str,
+        dest="UserAuth",
+        default="",
+        help="Enter your IRIS Authentification Username and Password " +
+        "(--User-Auth='username:authpassword') to " +
+        "access and download restricted data. " +
+        "[Default no user and password]")
+
+    # Database Settings
+    DataGroup = parser.add_argument_group(
+        title="Local Data Settings",
+        description="Settings associated with defining " +
+        "and using a local data base of pre-downloaded " +
+        "day-long SAC or MSEED files.")
+    DataGroup.add_argument(
+        "--local-data",
+        action="store",
+        type=str,
+        dest="localdata",
+        default=None,
+        help="Specify a comma separated list of paths containing " +
+        "day-long sac or mseed files of data already downloaded. " +
+        "If data exists for a seismogram is already present on disk, " +
+        "it is selected preferentially over downloading " +
+        "the data using the Client interface")
+    DataGroup.add_argument(
+    	"--dtype",
+    	action="store",
+    	type=str,
+    	dest="dtype",
+    	default='SAC',
+    	help="Specify the data archive file type, either SAC " +
+    	" or MSEED. Note the default behaviour is to search for " +
+    	"SAC files. Local archive files must have extensions of '.SAC' "+
+    	" or '.MSEED. These are case dependent, so specify the correct case"+
+        "here.")
+    DataGroup.add_argument(
+        "--no-data-zero",
+        action="store_true",
+        dest="ndval",
+        default=False,
+        help="Specify to force missing data to be set as zero, rather " +
+        "than default behaviour which sets to nan. Note this is applied " +
+        "only to the SAC data")
+    DataGroup.add_argument(
+        "--no-local-net",
+        action="store_false",
+        dest="useNet",
+        default=True,
+        help="Specify to prevent using the Network code in the " +
+        "search for local data (sometimes for CN stations " +
+        "the dictionary name for a station may disagree with that " +
+        "in the filename. [Default Network used]")
+    DataGroup.add_argument(
+        "--save-Z12",
+        action="store_true",
+        dest="saveZ12",
+        default=False,
+        help="Specify to save Z12 (un-rotated) components. [Default " +
+        "False]")
+
+    # Event Selection Criteria
+    EventGroup = parser.add_argument_group(
+        title="Event Settings",
+        description="Settings associated with refining " +
+        "the events to include in matching event-station pairs")
+    EventGroup.add_argument(
+        "--start",
+        action="store",
+        type=str,
+        dest="startT",
+        default="",
+        help="Specify a UTCDateTime compatible string representing " +
+        "the start time for the event search. This will override any " +
+        "station start times. [Default start date of station]")
+    EventGroup.add_argument(
+        "--end",
+        action="store",
+        type=str,
+        dest="endT",
+        default="",
+        help="Specify a UTCDateTime compatible string representing " +
+        "the end time for the event search. This will override any " +
+        "station end times [Default end date of station]")
+    EventGroup.add_argument(
+        "--reverse", "-R",
+        action="store_true",
+        dest="reverse",
+        default=False,
+        help="Reverse order of events. Default behaviour starts at " +
+        "oldest event and works towards most recent. Specify reverse " +
+        "order and instead the program will start with the most recent " +
+        "events and work towards older")
+    EventGroup.add_argument(
+        "--minmag",
+        action="store",
+        type=float,
+        dest="minmag",
+        default=6.0,
+        help="Specify the minimum magnitude of event for which to search. " +
+        "[Default 6.0]")
+    EventGroup.add_argument(
+        "--maxmag",
+        action="store",
+        type=float,
+        dest="maxmag",
+        default=9.0,
+        help="Specify the maximum magnitude of event for which to search. " +
+        "[Default None, i.e. no limit]")
+
+    # Geometry Settings
+    PhaseGroup = parser.add_argument_group(
+        title="Geometry Settings",
+        description="Settings associatd with the "
+        "event-station geometries for the specified phase")
+    PhaseGroup.add_argument(
+        "--phase",
+        action="store",
+        type=str,
+        dest="phase",
+        default='P',
+        help="Specify the phase name to use. Be careful with the distance. " +
+        "setting. Options are 'P' or 'PP'. [Default 'P']")
+    PhaseGroup.add_argument(
+        "--mindist",
+        action="store",
+        type=float,
+        dest="mindist",
+        default=None,
+        help="Specify the minimum great circle distance (degrees) between " +
+        "the station and event. [Default depends on phase]")
+    PhaseGroup.add_argument(
+        "--maxdist",
+        action="store",
+        type=float,
+        dest="maxdist",
+        default=None,
+        help="Specify the maximum great circle distance (degrees) between " +
+        "the station and event. [Default depends on phase]")
+
+    # Constants Settings
+    ConstGroup = parser.add_argument_group(
+        title='Parameter Settings',
+        description="Miscellaneous default values and settings")
+    ConstGroup.add_argument(
+        "--sampling-rate",
+        action="store",
+        type=float,
+        dest="new_sampling_rate",
+        default=10.,
+        help="Specify new sampling rate in Hz. [Default 10.]")
+    ConstGroup.add_argument(
+        "--dts",
+        action="store",
+        type=float,
+        dest="dts",
+        default=150.,
+        help="Specify the window length in sec (symmetric about arrival " +
+        "time). [Default 150.]")
+
+    args = parser.parse_args(argv)
+
+    # Check inputs
+    if not exist(args.indb):
+        parser.error("Input file " + args.indb + " does not exist")
+
+    # create station key list
+    if len(args.stkeys) > 0:
+        args.stkeys = args.stkeys.split(',')
+
+    # construct start time
+    if len(args.startT) > 0:
+        try:
+            args.startT = UTCDateTime(args.startT)
+        except:
+            parser.error(
+                "Cannot construct UTCDateTime from start time: " +
+                args.startT)
+    else:
+        args.startT = None
+
+    # construct end time
+    if len(args.endT) > 0:
+        try:
+            args.endT = UTCDateTime(args.endT)
+        except:
+            parser.error(
+                "Cannot construct UTCDateTime from end time: " +
+                args.endT)
+    else:
+        args.endT = None
+
+    # Parse User Authentification
+    if not len(args.UserAuth) == 0:
+        tt = args.UserAuth.split(':')
+        if not len(tt) == 2:
+            parser.error(
+                "Error: Incorrect Username and Password " +
+                "Strings for User Authentification")
+        else:
+            args.UserAuth = tt
+    else:
+        args.UserAuth = []
+
+    # Parse Local Data directories
+    if args.localdata is not None:
+        args.localdata = args.localdata.split(',')
+    else:
+        args.localdata = []
+
+    # Check Datatype specification
+    if (not args.dtype.upper() == 'MSEED') and  (not args.dtype.upper() == 'SAC'):
+    	parser.error(
+    		"Error: Local Data Archive must be of types 'SAC'" +
+    		"or MSEED. These must match the file extensions for " +
+    		" the archived data.")
+
+    # Check NoData Value
+    if args.ndval:
+        args.ndval = 0.0
+    else:
+        args.ndval = nan
+
+    # Check distances for selected phase
+    if args.phase not in ['P', 'PP', 'S', 'SKS']:
+        parser.error(
+            "Error: choose between 'P', 'PP', 'S' and 'SKS'.")
+    if args.phase == 'P':
+        if not args.mindist:
+            args.mindist = 30.
+        if not args.maxdist:
+            args.maxdist = 100.
+        if args.mindist < 30. or args.maxdist > 100.:
+            parser.error(
+                "Distances should be between 30 and 100 deg. for " +
+                "teleseismic 'P' waves.")
+    elif args.phase == 'PP':
+        if not args.mindist:
+            args.mindist = 100.
+        if not args.maxdist:
+            args.maxdist = 180.
+        if args.mindist < 100. or args.maxdist > 180.:
+            parser.error(
+                "Distances should be between 100 and 180 deg. for " +
+                "teleseismic 'PP' waves.")
+    elif args.phase == 'S':
+        if not args.mindist:
+            args.mindist = 55.
+        if not args.maxdist:
+            args.maxdist = 85.
+        if args.mindist < 55. or args.maxdist > 85.:
+            parser.error(
+                "Distances should be between 55 and 85 deg. for " +
+                "teleseismic 'S' waves.")
+    elif args.phase == 'SKS':
+        if not args.mindist:
+            args.mindist = 85.
+        if not args.maxdist:
+            args.maxdist = 115.
+        if args.mindist < 85. or args.maxdist > 115.:
+            parser.error(
+                "Distances should be between 85 and 115 deg. for " +
+                "teleseismic 'SKS' waves.")
+
+    return args
+
+
+def main():
+
+    print()
+    print("#####################################################################")
+    print("#        __                            _          _       _         #")
+    print("#  _ __ / _|_ __  _   _      __ _  ___| |_     __| | __ _| |_ __ _  #")
+    print("# | '__| |_| '_ \| | | |    / _` |/ _ \ __|   / _` |/ _` | __/ _` | #")
+    print("# | |  |  _| |_) | |_| |   | (_| |  __/ |_   | (_| | (_| | || (_| | #")
+    print("# |_|  |_| | .__/ \__, |____\__, |\___|\__|___\__,_|\__,_|\__\__,_| #")
+    print("#          |_|    |___/_____|___/        |_____|                    #")
+    print("#                                                                   #")
+    print("#####################################################################")
+    print()
+
+    # Run Input Parser
+    args = get_data_arguments()
+
+    # Load Database
+    db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys)
+
+    # Loop over station keys
+    for stkey in list(stkeys):
+
+        # Extract station information from dictionary
+        sta = db[stkey]
+
+        # Construct Folder Name
+        stfld = stkey
+        if not args.lkey:
+            stfld = stkey.split('.')[0]+"."+stkey.split('.')[1]
+
+        # Define path to see if it exists
+        if args.phase in ['P', 'PP']:
+            datapath = Path('P_DATA') / stfld
+        elif args.phase in ['S', 'SKS']:
+            datapath = Path('S_DATA') / stfld
+        if not datapath.exists():
+            print('Path to '+str(datapath)+' doesn`t exist - creating it')
+            datapath.mkdir(parents=True)
+
+        # Establish client for data
+        if len(args.UserAuth) == 0:
+            data_client = Client(args.Server)
+        else:
+            data_client = Client(args.Server, user=args.UserAuth[0],
+                                 password=args.UserAuth[1])
+
+        # Establish client for events
+        event_client = Client()
+
+        # Get catalogue search start time
+        if args.startT is None:
+            tstart = sta.startdate
+        else:
+            tstart = args.startT
+
+        # Get catalogue search end time
+        if args.endT is None:
+            tend = sta.enddate
+        else:
+            tend = args.endT
+        if tstart > sta.enddate or tend < sta.startdate:
+            continue
+
+        # Temporary print locations
+        tlocs = sta.location
+        if len(tlocs) == 0:
+            tlocs = ['']
+        for il in range(0, len(tlocs)):
+            if len(tlocs[il]) == 0:
+                tlocs[il] = "--"
+        sta.location = tlocs
+
+        # Update Display
+        print(" ")
+        print(" ")
+        print("|===============================================|")
+        print("|===============================================|")
+        print("|                   {0:>8s}                    |".format(
+            sta.station))
+        print("|===============================================|")
+        print("|===============================================|")
+        print("|  Station: {0:>2s}.{1:5s}                            |".format(
+            sta.network, sta.station))
+        print("|      Channel: {0:2s}; Locations: {1:15s}  |".format(
+            sta.channel, ",".join(tlocs)))
+        print("|      Lon: {0:7.2f}; Lat: {1:6.2f}                |".format(
+            sta.longitude, sta.latitude))
+        print("|      Start time: {0:19s}          |".format(
+            sta.startdate.strftime("%Y-%m-%d %H:%M:%S")))
+        print("|      End time:   {0:19s}          |".format(
+            sta.enddate.strftime("%Y-%m-%d %H:%M:%S")))
+        print("|-----------------------------------------------|")
+        print("| Searching Possible events:                    |")
+        print("|   Start: {0:19s}                  |".format(
+            tstart.strftime("%Y-%m-%d %H:%M:%S")))
+        print("|   End:   {0:19s}                  |".format(
+            tend.strftime("%Y-%m-%d %H:%M:%S")))
+        if args.maxmag is None:
+            print("|   Mag:   >{0:3.1f}", format(args.minmag) +
+                  "                                 |")
+        else:
+            msg = "|   Mag:   {0:3.1f}".format(args.minmag) + \
+                " - {0:3.1f}".format(args.maxmag) + \
+                "                            |"
+            print(msg)
+
+        print("| ...                                           |")
+
+        # Get catalogue using deployment start and end
+        cat = event_client.get_events(
+            starttime=tstart, endtime=tend,
+            minmagnitude=args.minmag, maxmagnitude=args.maxmag)
+
+        # Total number of events in Catalogue
+        nevK = 0
+        nevtT = len(cat)
+        print(
+            "|  Found {0:5d}".format(nevtT) +
+            " possible events                  |")
+        ievs = range(0, nevtT)
+
+        # Get Local Data Availabilty
+        if len(args.localdata) > 0:
+            print("|-----------------------------------------------|")
+            print("| Cataloging Local Data...                      |")
+            if args.useNet:
+                stalcllist = utils.list_local_data_stn(
+                    lcldrs=args.localdata, sta=sta.station,
+                    net=sta.network, dtype=args.dtype, altnet=sta.altnet)
+                print("|   {0:>2s}.{1:5s}: {2:6d}".format(
+                    sta.network, sta.station, len(stalcllist)) +
+                    " files                      |")
+                #print(stalcllist[0:10])
+            else:
+                stalcllist = utils.list_local_data_stn(
+                    lcldrs=args.localdata, sta=sta.station, dtype=args.dtype)
+                print("|   {0:5s}: {1:6d} files                " +
+                      "        |".format(
+                          sta.station, len(stalcllist)))
+        else:
+            stalcllist = []
+        print("|===============================================|")
+
+        # Select order of processing
+        if args.reverse:
+            ievs = range(0, nevtT)
+        else:
+            ievs = range(nevtT-1, -1, -1)
+
+        # Read through catalogue
+        for iev in ievs:
+
+            # Extract event
+            ev = cat[iev]
+
+            # Initialize RF object with station info
+            rfdata = RFData(sta)
+
+            # Add event to rfdata object
+            accept = rfdata.add_event(
+                ev, gacmin=args.mindist, gacmax=args.maxdist,
+                phase=args.phase, returned=True)
+
+            # Define time stamp
+            yr = str(rfdata.meta.time.year).zfill(4)
+            jd = str(rfdata.meta.time.julday).zfill(3)
+            hr = str(rfdata.meta.time.hour).zfill(2)
+
+            # If event is accepted (data exists)
+            if accept:
+
+                # Display Event Info
+                nevK = nevK + 1
+                if args.reverse:
+                    inum = iev + 1
+                else:
+                    inum = nevtT - iev + 1
+                print(" ")
+                print("*"*50)
+                print("* #{0:d} ({1:d}/{2:d}):  {3:13s} {4}".format(
+                    nevK, inum, nevtT, rfdata.meta.time.strftime(
+                        "%Y%m%d_%H%M%S"), stkey))
+                if args.verb:
+                    print("*   Phase: {}".format(args.phase))
+                    print("*   Origin Time: " +
+                          rfdata.meta.time.strftime("%Y-%m-%d %H:%M:%S"))
+                    print(
+                        "*   Lat: {0:6.2f};        Lon: {1:7.2f}".format(
+                            rfdata.meta.lat, rfdata.meta.lon))
+                    print(
+                        "*   Dep: {0:6.2f} km;     Mag: {1:3.1f}".format(
+                            rfdata.meta.dep, rfdata.meta.mag))
+                    print(
+                        "*   Dist: {0:7.2f} km;".format(rfdata.meta.epi_dist) +
+                        "   Epi dist: {0:6.2f} deg\n".format(rfdata.meta.gac) +
+                        "*   Baz:  {0:6.2f} deg;".format(rfdata.meta.baz) +
+                        "   Az: {0:6.2f} deg".format(rfdata.meta.az))
+
+                # Event Folder
+                timekey = rfdata.meta.time.strftime("%Y%m%d_%H%M%S")
+                evtdir = datapath / timekey
+                ZNEfile = evtdir / 'ZNE_Data.pkl'
+                metafile = evtdir / 'Meta_Data.pkl'
+                stafile = evtdir / 'Station_Data.pkl'
+
+                # Check if RF data already exist and overwrite has been set
+                if evtdir.exists():
+                    if RFfile.exists():
+                        if not args.ovr:
+                            continue
+
+                # Get data
+                has_data = rfdata.download_data(
+                    client=data_client, dts=args.dts, stdata=stalcllist,
+                    ndval=args.ndval, dtype=args.dtype, new_sr=args.new_sampling_rate,
+                    returned=True, verbose=args.verb)
+
+                if not has_data:
+                    continue
+
+                # Create Folder if it doesn't exist
+                if not evtdir.exists():
+                    evtdir.mkdir(parents=True)
+
+                # Save ZNE Traces
+                pickle.dump(rfdata.data, open(ZNEfile, "wb"))
+
+                # Save Z12 if components exist
+                if hasattr(rfdata, "dataZ12"):
+                    Z12file = evtdir / 'Z12_Data.pkl'
+                    pickle.dump(rfdata.dataZ12, open(Z12file, "wb"))
+
+                # Save event meta data
+                pickle.dump(rfdata.meta, open(metafile, "wb"))
+
+                # Save Station Data
+                pickle.dump(rfdata.sta, open(stafile, "wb"))
+
+                # Update
+                if args.verb:
+                    print("* Wrote Output Files to: ")
+                    print("*     "+str(evtdir))
+                print("*"*50)
+
+
+if __name__ == "__main__":
+
+    # Run main program
+    main()
diff --git a/rfpy/scripts/rfpy_plot.py b/rfpy/scripts/rfpy_plot.py
index 4715c64..99755b3 100644
--- a/rfpy/scripts/rfpy_plot.py
+++ b/rfpy/scripts/rfpy_plot.py
@@ -149,10 +149,10 @@ def get_plot_arguments(argv=None):
         action="store",
         dest="nbaz",
         type=int,
-        default=None,
+        default=36,
         help="Specify integer number of back-azimuth bins to consider " +
         "(typically 36 or 72). If not None, the plot will show receiver " +
-        "functions sorted by back-azimuth values. [Default None]")
+        "functions sorted by back-azimuth values. [Default 36]")
     PreGroup.add_argument(
         "--nslow",
         action="store",
@@ -301,9 +301,6 @@ def get_plot_arguments(argv=None):
         args.nbaz = 72
         print("'nbaz' or 'nslow' not specified - plotting using " +
               "'nbaz=72'")
-    elif args.nbaz is not None and args.nslow is not None:
-        parser.error(
-            "Error: Cannot specify both 'nbaz' and 'nslow'")
 
     if args.trange is None:
         args.trange = [0., 30.]
@@ -321,8 +318,6 @@ def get_plot_arguments(argv=None):
 
     if args.nbaz is None and args.nslow is None:
         parser.error("Specify at least one of --nbaz or --nslow")
-    elif args.nbaz is not None and args.nslow is not None:
-        parser.error("Specify only one of --nbaz or --nslow")
 
     return args
 
@@ -507,18 +502,24 @@ def main():
         print("Number of transverse RF data: " + str(len(rfTstream)))
         print('')
 
-        if args.nbaz:
+        if args.nbaz is not None and args.nslow is None:
             # Bin according to BAZ
             rf_tmp = binning.bin(rfRstream, rfTstream,
                                  typ='baz', nbin=args.nbaz+1,
                                  pws=args.pws)
 
-        elif args.nslow:
+        elif args.nslow is not None and args.nbaz is None:
             # Bin according to slowness
             rf_tmp = binning.bin(rfRstream, rfTstream,
                                  typ='slow', nbin=args.nslow+1,
                                  pws=args.pws)
 
+        elif args.nbaz is not None and args.nslow is not None:
+            # Bin into baz and slowness bins
+            rf_tmp = binning.bin_baz_slow(rfRstream, rfTstream, 
+                nbaz=args.nbaz+1, nslow=args.nslow+1,
+                pws=False)
+
         # Check bin counts:
         for tr in rf_tmp[0]:
             if (tr.stats.nbin < args.binlim):
@@ -559,19 +560,24 @@ def main():
             tr2 = None
 
         # Now plot
-        if args.nbaz:
+        if args.nbaz and args.nslow is None:
             plotting.wiggle_bins(rf_tmp[0], rf_tmp[1], tr1=tr1, tr2=tr2,
                                  btyp='baz', scale=args.scale,
                                  tmin=args.trange[0], tmax=args.trange[1],
                                  norm=norm, save=args.saveplot,
                                  title=args.titleplot, form=args.form)
-        elif args.nslow:
+        elif args.nslow and args.nbaz is None:
             plotting.wiggle_bins(rf_tmp[0], rf_tmp[1], tr1=tr1, tr2=tr2,
                                  btyp='slow', scale=args.scale,
                                  tmin=args.trange[0], tmax=args.trange[1],
                                  norm=norm, save=args.saveplot,
                                  title=args.titleplot, form=args.form)
 
+        elif args.nbaz is not None and args.nslow is not None:
+            plotting.panel(rf_tmp[0], rf_tmp[1], tmin=args.trange[0], tmax=args.trange[1],
+                normalize=True, save=args.saveplot,
+                title=args.titleplot, form=args.form)
+
         # Update processed folders
         procfold.append(stfld)
 
diff --git a/rfpy/utils.py b/rfpy/utils.py
index 72a98b8..246d6d8 100644
--- a/rfpy/utils.py
+++ b/rfpy/utils.py
@@ -2,7 +2,8 @@
 from obspy import UTCDateTime
 from numpy import nan, isnan, abs
 import numpy as np
-from obspy.core import Stream, read
+from obspy import Stream, Inventory
+from obspy import read, read_inventory
 
 
 def floor_decimal(n, decimals=0):
@@ -79,12 +80,16 @@ def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='SAC', altnet=[])
     fpathmatch = []
     # Loop over all local data directories
     for lcldr in lcldrs:
+        if not lcldr:
+            continue
         # Recursiely walk through directory
         for root, dirnames, filenames in walk(lcldr):
             # Keep paths only for those matching the station
             for sstring in sstrings:
                 for filename in filter(filenames, sstring):
-                    fpathmatch.append(join(root, filename))
+                    # No hidden directories
+                    if not '/.' in root:
+                        fpathmatch.append(join(root, filename))
 
     fpathmatch.sort()
 
@@ -135,49 +140,49 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
         ("*          {0:2s}{1:1s} - Checking Disk".format(sta.channel.upper(),
                                                           comp.upper())))
 
+    # Possible naming conventions
+    f1 = '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'
+    f2 = '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'
+    f3 = '*/{0:4s}.{1:3s}.*.{2:s}.{3:s}.*.{4:2s}{5:1s}*.{6:s}'
+    f4 = '*/{0:4s}.{1:3s}.*.{2:s}.{3:s}.*.*{4:1s}.D.{5:s}'
+    f5 = '*/{0:4s}.{1:3s}.*.??.{2:s}.*.*{3:1s}.D.{4:s}'
+    f6 = '*/{0:}/{1:}/{2:}/{3:}{4:}.D/{1:}.{2:}.*.*.D.{0:}.{5:}'
+
     # Time Window Spans Single Day
     if stjd == edjd:
-        # Format 1
-        lclfiles = list(filter(
-            stdata,
-            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'.format(
-                styr, stjd, sta.network.upper(
-                ), sta.station.upper(), sta.channel.upper()[0:2],
-                comp.upper(), dtype)))
-        # Format 2
-        if len(lclfiles) == 0:
-            lclfiles = list(filter(
-                stdata,
-                '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format(
-                    styr, stjd, sta.network.upper(), sta.station.upper(),
-                    comp.upper(), dtype)))
 
-        # Alternate Nets (for CN/PO issues) Format 1
-        if len(lclfiles) == 0:
-            lclfiles = []
-            for anet in sta.altnet:
-                lclfiles.extend(
-                    list(
-                        filter(
-                            stdata,
-                            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.' +
-                            '{4:2s}{5:1s}.{6:s}'.format(
-                                styr, stjd, anet.upper(), sta.station.upper(),
-                                sta.channel.upper()[0:2], comp.upper(), dtype))))
-
-        # Alternate Nets (for CN/PO issues) Format 2
-        if len(lclfiles) == 0:
-            # Check Alternate Networks
-            lclfiles = []
-            for anet in sta.altnet:
-                lclfiles.extend(
-                    list(
-                        filter(
-                            stdata,
-                            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*' +
-                            '{4:1s}.{5:s}'.format(
-                                styr, stjd, sta.network.upper(),
-                                sta.station.upper(), comp.upper(), dtype))))
+        lclfiles = []
+        nets = [sta.network] + sta.altnet
+        for net in nets:
+            # Start day
+            s1 = f1.format(styr, stjd, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), dtype)
+            s2 = f2.format(styr, stjd, net.upper(), sta.station.upper(),
+                           comp.upper(), dtype)
+            s3 = f3.format(styr, stjd, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), dtype)
+            s4 = f4.format(styr, stjd, net.upper(), sta.station.upper(),
+                           comp.upper(), dtype)
+            s5 = f5.format(styr, stjd, sta.station.upper(),
+                           comp.upper(), dtype)
+            s6 = f6.format(styr,net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), stjd)
+
+            print("*          Trying formats:")
+            print("*          " + s1)
+            print("*          " + s2)
+            print("*          " + s3)
+            print("*          " + s4)
+            print("*          " + s5)
+            print("*          " + s6)
+            print("*          ")
+
+            lclfiles.extend(list(filter(stdata, s1)))
+            lclfiles.extend(list(filter(stdata, s2)))
+            lclfiles.extend(list(filter(stdata, s3)))
+            lclfiles.extend(list(filter(stdata, s4)))
+            lclfiles.extend(list(filter(stdata, s5)))
+            lclfiles.extend(list(filter(stdata, s6)))
 
         # If still no Local files stop
         if len(lclfiles) == 0:
@@ -187,8 +192,12 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
         # Process the local Files
         for sacfile in lclfiles:
             # Read File
-            st = read(sacfile)
-            # st = read(sacfile, format="SAC")
+            try:
+                st = read(sacfile)
+            except OSError:
+                print("*              - Met OSError.")
+                print(f"*              - Possibly corrupt file: {sacfile}.")
+                return False, Stream()
 
             if dtype.upper() == 'MSEED':
                 if len(st) > 1:
@@ -209,7 +218,10 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
                     eddt = False
                     # Check for NoData and convert to NaN if a SAC file
                     if dtype.upper() == 'SAC':
-                        stnd = st[0].stats.sac['user9']
+                        try:
+                            stnd = st[0].stats.sac['user9']
+                        except (KeyError, AttributeError):
+                            stnd = 0.0
                         if (not stnd == 0.0) and (not stnd == -12345.0):
                             st[0].data[st[0].data == stnd] = ndval
                             eddt = True
@@ -244,83 +256,51 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
 
     # Time Window spans Multiple days
     else:
-        # Day 1 Format 1
-        lclfiles1 = list(
-            filter(stdata,
-                   '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'.format(
-                       styr, stjd, sta.network.upper(), sta.station.upper(),
-                       sta.channel.upper()[0:2], comp.upper(), dtype)))
-        # Day 1 Format 2
-        if len(lclfiles1) == 0:
-            lclfiles1 = list(
-                filter(stdata,
-                       '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format(
-                           styr, stjd, sta.network.upper(),
-                           sta.station.upper(), comp.upper(), dtype)))
-        # Day 1 Alternate Nets (for CN/PO issues) Format 1
-        if len(lclfiles1) == 0:
-            lclfiles1 = []
-            for anet in sta.altnet:
-                lclfiles1.extend(
-                    list(
-                        filter(
-                            stdata,
-                            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.' +
-                            '{4:2s}{5:1s}.{6:s}'.format(
-                                styr, stjd, anet.upper(), sta.station.upper(
-                                ), sta.channel.upper()[0:2],
-                                comp.upper(), dtype))))
-        # Day 1 Alternate Nets (for CN/PO issues) Format 2
-        if len(lclfiles1) == 0:
-            lclfiles1 = []
-            for anet in sta.altnet:
-                lclfiles1.extend(
-                    list(
-                        filter(
-                            stdata,
-                            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format(
-                                styr, stjd, anet.upper(),
-                                sta.station.upper(), comp.upper(), dtype))))
-
-        # Day 2 Format 1
-        lclfiles2 = list(
-            filter(stdata,
-                   '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'.format(
-                       edyr, edjd, sta.network.upper(
-                       ), sta.station.upper(), sta.channel.upper()[0:2],
-                       comp.upper(), dtype)))
-        # Day 2 Format 2
-        if len(lclfiles2) == 0:
-            lclfiles2 = list(
-                filter(stdata,
-                       '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*' +
-                       '{4:1s}.{5:s}'.format(
-                           edyr, edjd, sta.network.upper(),
-                           sta.station.upper(),
-                           comp.upper(), dtype)))
-        # Day 2 Alternate Nets (for CN/PO issues) Format 1
-        if len(lclfiles2) == 0:
-            lclfiles2 = []
-            for anet in sta.altnet:
-                lclfiles2.extend(
-                    list(
-                        filter(
-                            stdata,
-                            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.' +
-                            '{4:2s}{5:1s}.{6:s}'.format(
-                                edyr, edjd, anet.upper(), sta.station.upper(),
-                                sta.channel.upper()[0:2], comp.upper(), dtype))))
-        # Day 2 Alternate Nets (for CN/PO issues) Format 2
-        if len(lclfiles2) == 0:
-            lclfiles2 = []
-            for anet in sta.altnet:
-                lclfiles2.extend(
-                    list(
-                        filter(
-                            stdata,
-                            '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format(
-                                edyr, edjd, anet.upper(), sta.station.upper(),
-                                comp.upper(), dtype))))
+        lclfiles1 = []
+        lclfiles2 = []
+        nets = [sta.network] + sta.altnet
+        for net in nets:
+            # Start day
+            s1 = f1.format(styr, stjd, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), dtype)
+            s2 = f2.format(styr, stjd, net.upper(), sta.station.upper(),
+                           comp.upper(), dtype)
+            s3 = f3.format(styr, stjd, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), dtype)
+            s4 = f4.format(styr, stjd, net.upper(), sta.station.upper(),
+                           comp.upper(), dtype)
+            s5 = f5.format(styr, stjd, sta.station.upper(),
+                           comp.upper(), dtype)
+            s6 = f6.format(styr,net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), stjd)
+
+            lclfiles1.extend(list(filter(stdata, s1)))
+            lclfiles1.extend(list(filter(stdata, s2)))
+            lclfiles1.extend(list(filter(stdata, s3)))
+            lclfiles1.extend(list(filter(stdata, s4)))
+            lclfiles1.extend(list(filter(stdata, s5)))
+            lclfiles1.extend(list(filter(stdata, s6)))
+
+            # End day
+            s1 = f1.format(edyr, edjd, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), dtype)
+            s2 = f2.format(edyr, edjd, net.upper(), sta.station.upper(),
+                           comp.upper(), dtype)
+            s3 = f3.format(edyr, edjd, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), dtype)
+            s4 = f4.format(edyr, edjd, net.upper(), sta.station.upper(),
+                           comp.upper(), dtype)
+            s5 = f5.format(edyr, edjd, sta.station.upper(),
+                           comp.upper(), dtype)
+            s6 = f6.format(edyr, net.upper(), sta.station.upper(),
+                           sta.channel.upper()[0:2], comp.upper(), edjd)
+
+            lclfiles2.extend(list(filter(stdata, s1)))
+            lclfiles2.extend(list(filter(stdata, s2)))
+            lclfiles2.extend(list(filter(stdata, s3)))
+            lclfiles2.extend(list(filter(stdata, s4)))
+            lclfiles2.extend(list(filter(stdata, s5)))
+            lclfiles2.extend(list(filter(stdata, s6)))
 
         # If still no Local files stop
         if len(lclfiles1) == 0 and len(lclfiles2) == 0:
@@ -334,8 +314,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
                 st1 = read(sacf1)
                 if dtype.upper() == 'MSEED':
                     if len(st1) > 1:
-                        st1.merge(method=1, interpolation_samples=-
-                                  1, fill_value=-123456789)
+                        st1.merge(method=1, interpolation_samples=-1,
+                                  fill_value=-123456789)
 
                 # Loop over second day file options
                 for sacf2 in lclfiles2:
@@ -343,7 +323,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
                     if dtype.upper() == 'MSEED':
                         if len(st2) > 1:
                             st2.merge(
-                                method=1, interpolation_samples=-1, fill_value=-123456789)
+                                method=1, interpolation_samples=-1,
+                                fill_value=-123456789)
 
                     # Check time overlap of the two files.
                     if st1[0].stats.endtime >= \
@@ -383,7 +364,10 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
                                     eddt = False
                                     # Check for NoData and convert to NaN if a SAC file
                                     if dtype.upper() == 'SAC':
-                                        stnd = st[0].stats.sac['user9']
+                                        try:
+                                            stnd = st[0].stats.sac['user9']
+                                        except KeyError:
+                                            stnd = 0.0
                                         if (not stnd == 0.0) and (not stnd == -12345.0):
                                             st[0].data[st[0].data == stnd] = ndval
                                             eddt = True
@@ -434,8 +418,36 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None,
     return erd, None
 
 
+def attach_local_response_data(stream, local_response_dir):
+    """
+    Function to restrieve response data from locally stored dataless seed
+
+    Parameters
+    ----------
+    stream : obspy.Stream
+        Event seismogram
+    local_response_dir: str
+        Directory holding response information. All files containing the station
+        name are read. If that fails, all files beginning with the network code are
+        read.
+    """
+    stations = set([t.stats.station for t in stream.traces])
+    networks = set([t.stats.network for t in stream.traces])
+    inventory = Inventory()
+    try:
+        for station in stations:
+            inventory += read_inventory(
+                '{:}/*{:}*'.format(local_response_dir, station))
+    except Exception:
+        for net in networks:
+            inventory += read_inventory('{:}/{:}*'.format(local_response_dir, net))
+
+    stream.attach_response(inventories=inventory)
+
+
 def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime,
-                  stdata=[], dtype='SAC', ndval=nan, new_sr=0., verbose=False):
+                  stdata=[], dtype='SAC', ndval=nan, new_sr=0., verbose=False,
+                  remove_response=False, local_response_dir=''):
     """
     Function to build a stream object for a seismogram in a given time window either
     by downloading data from the client object or alternatively first checking if the
@@ -456,9 +468,16 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime,
     end : :class:`~obspy.core.utcdatetime.UTCDateTime`
         End time for request
     stdata : List
-        Station list
+        List of directories holding local waveform data
     ndval : float or nan
         Default value for missing data
+    remove_response : bool
+        Remove instrument response from seismogram and resitute to true ground
+        velocity (m/s) using obspy.core.trace.Trace.remove_response()
+    local_response_dir: str
+        Directory holding response files to be read by obspy.read_inventory().
+        Required when ``remove_response`` and using locally stored waveform data
+        via ``stdata``.
 
     Returns
     -------
@@ -506,6 +525,15 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime,
         if not erd:
             # Combine Data
             st = stZ + stN + stE
+            if remove_response:
+                if local_response_dir:
+                    attach_local_response_data(st, local_response_dir)
+                else:
+                    print('*')
+                    print('* Warning: No local_response_dir given.')
+                    print('* Warning: Continuing without.')
+                    print('*')
+
 
     # No local data? Request using client
     if erd:
@@ -516,115 +544,223 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime,
             # Construct location name
             if len(tloc) == 0:
                 tloc = "--"
-            # Construct Channel List
-            channelsZNE = sta.channel.upper() + 'Z,' + sta.channel.upper() + \
-                'N,' + sta.channel.upper() + 'E'
-            print(("*          {1:2s}[ZNE].{2:2s} - Checking Network".format(
-                sta.station, sta.channel.upper(), tloc)))
 
-            # Get waveforms, with extra 1 second to avoid
-            # traces cropped too short - traces are trimmed later
-            try:
-                st = client.get_waveforms(
-                    network=sta.network,
-                    station=sta.station, location=loc,
-                    channel=channelsZNE, starttime=start,
-                    endtime=end+1., attach_response=False)
-                if len(st) == 3:
-                    print("*              - ZNE Data Downloaded")
-
-                # It's possible if len(st)==1 that data is Z12
+            # Construct Channel List
+            chaZNE = (sta.channel.upper() + 'Z,' +
+                      sta.channel.upper() + 'N,' +
+                      sta.channel.upper() + 'E')
+            msgZNE = "*          {1:2s}[ZNE].{2:2s} - Checking Network".format(
+                sta.station, sta.channel.upper(), tloc)
+
+            chaZ12 = (sta.channel.upper() + 'Z,' +
+                      sta.channel.upper() + '1,' +
+                      sta.channel.upper() + '2')
+            msgZ12 = "*          {1:2s}[Z12].{2:2s} - Checking Network".format(
+                sta.station, sta.channel.upper(), tloc)
+
+            # Loop over possible channels
+            for channel, msg in zip([chaZNE, chaZ12], [msgZNE, msgZ12]):
+                print(msg)
+
+                kwargs = {}
+                if remove_response:
+                    kwargs = dict(attach_response=remove_response)
+
+                # Get waveforms, with extra 1 second to avoid
+                # traces cropped too short - traces are trimmed later
+                try:
+                    st = client.get_waveforms(
+                        network=sta.network,
+                        station=sta.station, location=loc,
+                        channel=channel, starttime=start,
+                        endtime=end+1., **kwargs)
+                except Exception as e:
+                    if verbose:
+                        print("* Met exception:")
+                        print("* " + e.__repr__())
+                    st = None
                 else:
-                    # Construct Channel List
-                    channelsZ12 = sta.channel.upper() + 'Z,' + \
-                        sta.channel.upper() + '1,' + \
-                        sta.channel.upper() + '2'
-                    msg = "*          {1:2s}[Z12].{2:2s} - Checking Network".format(
-                        sta.station, sta.channel.upper(), tloc)
-                    print(msg)
-                    try:
-                        st = client.get_waveforms(
-                            network=sta.network,
-                            station=sta.station, location=loc,
-                            channel=channelsZ12, starttime=start,
-                            endtime=end+1., attach_response=False)
-                        if len(st) == 3:
-                            print("*              - Z12 Data Downloaded")
-                        else:
-                            st = None
-                    except:
-                        st = None
-            except:
-                st = None
+                    if len(st) == 3:
+                        # It's possible if len(st)==1 that data is Z12
+                        print("*              - ZNE Data Downloaded")
+                        break
 
             # Break if we successfully obtained 3 components in st
             if not erd:
-
                 break
 
     # Check the correct 3 components exist
-    if st is None:
+    if st is None or len(st) < 3:
         print("* Error retrieving waveforms")
         print("**************************************************")
         return True, None
 
     # Three components successfully retrieved
+    if remove_response:
+        st.remove_response()
+        print("*")
+        print("* Restituted stream to true ground velocity.")
+        print("*")
+
+    # Detrend and apply taper
+    st.detrend('demean').detrend('linear').taper(
+        max_percentage=0.05, max_length=5.)
+
+    # Check start times
+    if not np.all([tr.stats.starttime == start for tr in st]):
+        print("* Start times are not all close to true start: ")
+        [print("*   "+tr.stats.channel+" " +
+               str(tr.stats.starttime)+" " +
+               str(tr.stats.endtime)) for tr in st]
+        print("*   True start: "+str(start))
+        print("* -> Shifting traces to true start")
+        delay = [tr.stats.starttime - start for tr in st]
+        st_shifted = Stream(
+            traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)])
+        st = st_shifted.copy()
+
+    # Check sampling rate
+    sr = st[0].stats.sampling_rate
+    sr_round = float(floor_decimal(sr, 0))
+    if not sr == sr_round:
+        print("* Sampling rate is not an integer value: ", sr)
+        print("* -> Resampling")
+        st.resample(sr_round, no_filter=False)
+
+    # Try trimming
+    try:
+        st.trim(start, end)
+    except:
+        print("* Unable to trim")
+        print("* -> Skipping")
+        print("**************************************************")
+        return True, None
+
+    # Check final lengths - they should all be equal if start times
+    # and sampling rates are all equal and traces have been trimmed
+    if not np.allclose([tr.stats.npts for tr in st[1:]], st[0].stats.npts):
+        print("* Lengths are incompatible: ")
+        [print("*     "+str(tr.stats.npts)) for tr in st]
+        print("* -> Skipping")
+        print("**************************************************")
+
+        return True, None
+
+    elif not np.allclose([st[0].stats.npts], int((end - start)*sr),
+                         atol=1):
+        print("* Length is too short: ")
+        print("*    "+str(st[0].stats.npts) +
+              " ~= "+str(int((end - start)*sr)))
+        print("* -> Skipping")
+        print("**************************************************")
+
+        return True, None
+
     else:
+        print("* Waveforms Retrieved...")
+        return False, st
+
+
+def optimal_wlevel(rfdatas, wlevels=10**(np.arange(-5, 0.25, 0.25)), iplot=False):
+    """
+    Optimal water level for 'water2' deconvolution using Generalized Cross
+    Validation
+
+    Accepts
+    -------
+
+    rfdatas : list of rfpy.RFData
+        rfdata elements hold rotated seismograms in rfdata.data attribute
+    wlevels : iterable
+        Range of water levels (fraction of maximum stacked amplitude) to try.
+        Default: 1e-5 to 10
+
+    Returns:
+    --------
+    wlevel : float
+        Optimal water level (fraction of maximum stacked amplitude)
+    """
+
+    nt = len(rfdatas)
+    nf = len(rfdatas[0].specs[0].data)
+
+    wft = np.zeros((nt, nf), dtype=complex)
+    vft = np.zeros((nt, nf), dtype=complex)
+
+    for it, rfdata in enumerate(rfdatas):
+
+        if not rfdata.meta.rotated:
+            msg = 'Element {:} of rfdatas must be rotated.'.format(it)
+            raise ValueError(msg)
+
+        wft[it, :] = np.fft.fft(rfdata.data[0].data, n=nf)
+        vft[it, :] = np.fft.fft(rfdata.data[1].data, n=nf)
+
+    # Auto and cross spectra
+    wwft = wft*np.conjugate(wft)
+    vwft = vft*np.conjugate(wft)
+
+    if nt > 0:
+        wwft = np.sum(wwft)
+        vwft = np.sum(vwft)
+
+    betas = wlevels*np.amax(wwft)
+    misfits = np.zeros_like(betas)
+    modnorm = np.zeros_like(betas)
+    gcvf = np.zeros_like(betas)
+
+    for ib, beta in enumerate(betas):
+        # Define operator W W* / (W W* + B) and deconvolve to get
+        # impulse response in frequency domain.
+        wwft2 = wwft + beta
+        rft = vwft / wwft2
+        xft = wwft / wwft2
 
-        # Detrend and apply taper
-        st.detrend('demean').detrend('linear').taper(
-            max_percentage=0.05, max_length=5.)
-
-        # Check start times
-        if not np.all([tr.stats.starttime == start for tr in st]):
-            print("* Start times are not all close to true start: ")
-            [print("*   "+tr.stats.channel+" " +
-                   str(tr.stats.starttime)+" " +
-                   str(tr.stats.endtime)) for tr in st]
-            print("*   True start: "+str(start))
-            print("* -> Shifting traces to true start")
-            delay = [tr.stats.starttime - start for tr in st]
-            st_shifted = Stream(
-                traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)])
-            st = st_shifted.copy()
-
-        # Check sampling rate
-        sr = st[0].stats.sampling_rate
-        sr_round = float(floor_decimal(sr, 0))
-        if not sr == sr_round:
-            print("* Sampling rate is not an integer value: ", sr)
-            print("* -> Resampling")
-            st.resample(sr_round, no_filter=False)
-
-        # Try trimming
-        try:
-            st.trim(start, end)
-        except:
-            print("* Unable to trim")
-            print("* -> Skipping")
-            print("**************************************************")
-            return True, None
-
-        # Check final lengths - they should all be equal if start times
-        # and sampling rates are all equal and traces have been trimmed
-        if not np.allclose([tr.stats.npts for tr in st[1:]], st[0].stats.npts):
-            print("* Lengths are incompatible: ")
-            [print("*     "+str(tr.stats.npts)) for tr in st]
-            print("* -> Skipping")
-            print("**************************************************")
-
-            return True, None
-
-        elif not np.allclose([st[0].stats.npts], int((end - start)*sr),
-                             atol=1):
-            print("* Length is too short: ")
-            print("*    "+str(st[0].stats.npts) +
-                  " ~= "+str(int((end - start)*sr)))
-            print("* -> Skipping")
-            print("**************************************************")
-
-            return True, None
+        # Compute model norm.
+        modnorm[ib] = np.linalg.norm(rft)**2
 
+        # Compute data misfit. Note misfit is numerator of GCV function
+        # Note also earlier mistake where norm(nft)^2 was norm(nft).
+        if nt == 1:
+            nft = vft - wft*rft  # this should be SP
+            misfits[ib] = np.linalg.norm(nft)**2
         else:
-            print("* Waveforms Retrieved...")
-            return False, st
+            for it in range(nt):
+                nft = vft[it, :] - wft[it, :]*rft
+                misfits[ib] = misfits[ib] + np.linalg.norm(nft)**2
+
+        # Compute denominator and GCV function.
+        den = (nf*nt - np.real(np.sum(xft)))**2
+        gcvf[ib] = misfits[ib]/den
+
+    # Compute best beta.
+    ibest = np.argmin(gcvf)
+
+    if iplot:
+        import matplotlib.pyplot as mp
+
+        fig, axs = mp.subplots(ncols=2)
+        ax = axs[0]
+        ax.plot(modnorm, misfits, color='black')
+        ax.plot(modnorm, misfits, marker='+', color='red')
+        ax.plot(modnorm[ibest], misfits[ibest], marker='o', color='green')
+        ax.set_xlabel('Model Norm')
+        ax.set_ylabel('Data Misfit')
+        ax = axs[1]
+        ax.plot(betas, gcvf)
+        ax.plot(betas, gcvf, marker='+', color='red')
+        ax.plot(betas[ibest], gcvf[ibest], color='green', marker='o')
+        ax.set_xscale('log')
+        ax.set_xlabel('Regularization Parameter')
+        ax.set_ylabel('GCV Function')
+        fig.show()
+        input('Press key to continue')
+        mp.close(fig)
+
+    if ibest == 0 or ibest == len(betas):
+        # Minimum not found
+        msg = 'No minimum found. Try extending wlevels.'
+        raise ValueError(msg)
+
+    wlevel = float(betas[ibest]/np.amax(wwft))
+    return wlevel
diff --git a/setup.py b/setup.py
index c9f78dc..25d6287 100644
--- a/setup.py
+++ b/setup.py
@@ -35,10 +35,14 @@ def find_version(*paths):
     install_requires=['numpy', 'obspy', 'stdb>=0.2.0'],
     python_requires='>=3.6',
     packages=setuptools.find_packages(),
-    entry_points={'console_scripts':[
-    'rfpy_calc=rfpy.scripts.rfpy_calc:main',
-    'rfpy_recalc=rfpy.scripts.rfpy_recalc:main',
-    'rfpy_plot=rfpy.scripts.rfpy_plot:main',
-    'rfpy_harmonics=rfpy.scripts.rfpy_harmonics:main',
-    'rfpy_hk=rfpy.scripts.rfpy_hk:main',
-    'rfpy_ccp=rfpy.scripts.rfpy_ccp:main']})
+    include_package_data=True,
+    package_data={'': ['examples/data/*.pkl']},
+    entry_points={
+        'console_scripts':
+            ['rfpy_get_data=rfpy.scripts.rfpy_get_data:main',
+             'rfpy_calc_simdec=rfpy.scripts.rfpy_calc_simdec:main',
+             'rfpy_calc=rfpy.scripts.rfpy_calc:main',
+             'rfpy_plot=rfpy.scripts.rfpy_plot:main',
+             'rfpy_harmonics=rfpy.scripts.rfpy_harmonics:main',
+             'rfpy_hk=rfpy.scripts.rfpy_hk:main',
+             'rfpy_ccp=rfpy.scripts.rfpy_ccp:main']})