diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 00000000..3d47d0b7
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,4 @@
+v0.1.0 (2018-07-18)
+    * Initial release; implements core functionality, 
+      documentation, and tests.
+    * Includes basic pipeline run scripts.
diff --git a/README.md b/README.md
index 3df6af03..d5cc3abe 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@ For usage examples and documentation, see http://hera-pspec.readthedocs.io/en/la
 
 ### Code Dependencies
 
-* numpy >= 1.10
+* numpy >= 1.14
 * pyuvdata (`pip install pyuvdata` or use https://github.com/HERA-Team/pyuvdata.git)
 * aipy (```conda install -c conda-forge aipy```)
 * scipy >= 0.19
diff --git a/hera_pspec/VERSION b/hera_pspec/VERSION
index 77d6f4ca..6e8bf73a 100644
--- a/hera_pspec/VERSION
+++ b/hera_pspec/VERSION
@@ -1 +1 @@
-0.0.0
+0.1.0
diff --git a/hera_pspec/grouping.py b/hera_pspec/grouping.py
index fa1b3203..f91ab55f 100644
--- a/hera_pspec/grouping.py
+++ b/hera_pspec/grouping.py
@@ -838,7 +838,7 @@ def bootstrap_average_blpairs(uvp_list, blpair_groups, time_avg=False,
 
 
 def bootstrap_resampled_error(uvp, blpair_groups=None, time_avg=False, Nsamples=1000, seed=None,
-                              normal_std=True, robust_std=True, cintervals=None, bl_error_tol=1.0,
+                              normal_std=True, robust_std=False, cintervals=None, bl_error_tol=1.0,
                               add_to_history='', verbose=False):
     """
     Given a UVPSpec object, generate bootstrap resamples of its average
diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py
index ae75aa79..13bdacfa 100644
--- a/hera_pspec/plot.py
+++ b/hera_pspec/plot.py
@@ -10,7 +10,8 @@
 def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False, 
                    average_times=False, fold=False, plot_noise=False, 
                    delay=True, deltasq=False, legend=False, ax=None,
-                   component='real'):
+                   component='real', lines=True, markers=False, error=None,
+                   **kwargs):
     """
     Plot a 1D delay spectrum (or spectra) for a group of baselines.
     
@@ -56,10 +57,25 @@ def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False,
         this case, even if the existing plot has completely different axis 
         labels etc.) If None, a new Axes object will be created. Default: None.
 
-    component : str
+    component : str, optional
         Component of complex spectra to plot, options=['abs', 'real', 'imag'].
         Default: 'real'. 
 
+    lines : bool, optional
+        If True, plot lines between bandpowers for a given pspectrum.
+        Default: True.
+
+    markers : bool, optional
+        If True, plot circles at bandpowers. Filled for positive, empty
+        for negative. Default: False.
+
+    error : str, optional
+        If not None and if error exists in uvp stats_array, plot errors
+        on bandpowers. Default: None.
+
+    kwargs : dict, optional
+        Extra kwargs to pass to _all_ ax.plot calls.
+
     Returns
     -------
     fig : matplotlib.pyplot.Figure
@@ -112,20 +128,61 @@ def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False,
     for blgrp in blpairs:
         # Loop over blpairs in group and plot power spectrum for each one
         for blp in blgrp:
+            # setup key and casting function
             key = (spw, blp, pol)
             if component == 'real':
-                power = np.abs(np.real(uvp_plt.get_data(key))).T
+                cast = np.real
             elif component == 'imag':
-                power = np.abs(np.imag(uvp_plt.get_data(key))).T
+                cast = np.imag
             elif component == 'abs':
-                power = np.abs(uvp_plt.get_data(key)).T
+                cast = np.abs
+
+            # get power array and repeat x array
+            power = cast(uvp_plt.get_data(key))
 
-            # flag spectra that have zero integration
+            # flag records that have zero integration
             flags = np.isclose(uvp_plt.get_integrations(key), 0.0)
-            power[:, flags] = np.nan
+            power[flags] = np.nan
+
+            # get errs if requessted
+            if error is not None and hasattr(uvp_plt, 'stats_array'):
+                if error in uvp_plt.stats_array:
+                    errs = uvp_plt.get_stats(error, key)
+                    errs[flags] = np.nan
+
+            # iterate over integrations per blp
+            for i in range(power.shape[0]):
+                # get y data
+                y = power[i]
+
+                # plot elements
+                cax = None
+                if lines:
+                    cax, = ax.plot(x, np.abs(y), label="%s" % str(key), marker='None', **kwargs)
+
+                if markers:
+                    if cax is None:
+                        c = None
+                    else:
+                        c = cax.get_color()
+                    # plot positive w/ filled circles
+                    cax, = ax.plot(x[y >= 0], np.abs(y[y >= 0]), c=c, ls='None', marker='o', 
+                                  markerfacecolor=c, markeredgecolor=c, **kwargs)
+                    # plot negative w/ unfilled circles
+                    c = cax.get_color()
+                    cax, = ax.plot(x[y < 0], np.abs(y[y < 0]), c=c, ls='None', marker='o',
+                                   markerfacecolor='None', markeredgecolor=c, **kwargs)
+
+                if error is not None and hasattr(uvp_plt, 'stats_array'):
+                    if error in uvp_plt.stats_array:
+                        if cax is None:
+                            c = None
+                        else:
+                            c = cax.get_color()
+                        cax = ax.errorbar(x, np.abs(y), fmt='none', ecolor=c, yerr=cast(errs[i]), **kwargs)
+                    else:
+                        raise KeyError("Error variable '%s' not found in stats_array of UVPSpec object." % error)
 
-            ax.plot(x, power, label="%s" % str(key))
-            
             # If blpairs were averaged, only the first blpair in the group 
             # exists any more (so skip the rest)
             if average_blpairs: break
@@ -329,18 +386,20 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa
         Nkeys = len(waterfall)
         Nside = int(np.ceil(np.sqrt(Nkeys)))
         fig, axes = plt.subplots(Nside, Nside, figsize=figsize)
+    
     # Ensure axes is an ndarray
     if isinstance(axes, matplotlib.axes._subplots.Axes):
         axes = np.array([[axes]])
     if isinstance(axes, list):
         axes = np.array(axes)
-    # ensure its 2D and get side lengths
+    
+    # Ensure its 2D and get side lengths
     if axes.ndim == 1:
         axes = axes[:, None]
     assert axes.ndim == 2, "input axes must have ndim == 2"
     Nvert, Nhorz = axes.shape
 
-    # get LST range: setting y-ticks is tricky due to LST wrapping...
+    # Get LST range: setting y-ticks is tricky due to LST wrapping...
     y = uvp_plt.lst_avg_array[uvp_plt.key_to_indices(waterfall.keys()[0])[1]]
     if lst_in_hrs:
         lst_units = "Hr"
@@ -369,7 +428,7 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa
         psunits = psunits.replace("beam normalization not specified", 
                                  r"{\rm unnormed}")
 
-    # iterate over waterfall keys
+    # Iterate over waterfall keys
     keys = waterfall.keys()
     k = 0
     for i in range(Nvert):
@@ -402,7 +461,7 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa
             # configure left-column plots
             if j == 0:
                 # set yticks
-                ax.set_yticks(np.arange(Ny//Ny_thin+1))
+                ax.set_yticks(np.arange(Ny)[::Ny_thin])
                 ax.set_yticklabels(y[::Ny_thin])
                 ax.set_ylabel(r"LST [{}]".format(lst_units), fontsize=16)
             else:
@@ -428,7 +487,8 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa
             units = "$%sP(k_\parallel)$ $[%s]$" % (logunits, psunits)
 
         spwrange = np.around(np.array(uvp_plt.get_spw_ranges()[spw][:2]) / 1e6, 2)
-        axes[0][0].get_figure().suptitle("{}\n{} polarization | {} -- {} MHz".format(units, pol, *spwrange), y=1.03, fontsize=14)
+        axes[0][0].get_figure().suptitle("{}\n{} polarization | {} -- {} MHz".format(units, pol, *spwrange), 
+                                         y=1.03, fontsize=14)
 
     # Return Axes
     if new_plot:
diff --git a/hera_pspec/tests/test_plot.py b/hera_pspec/tests/test_plot.py
index 44566d70..2f678fc7 100644
--- a/hera_pspec/tests/test_plot.py
+++ b/hera_pspec/tests/test_plot.py
@@ -4,7 +4,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import os, copy, sys
-from hera_pspec import pspecdata, pspecbeam, conversions, plot, utils
+from hera_pspec import pspecdata, pspecbeam, conversions, plot, utils, grouping
 from hera_pspec.data import DATA_PATH
 from pyuvdata import UVData
 
@@ -142,6 +142,19 @@ def test_plot_average(self):
         self.assertTrue( axes_contains(f4.axes[0], elements) )
         plt.close(f5)
         
+        # test errorbar plotting w/ markers
+
+        # bootstrap resample
+        (uvp_avg, _,
+         _) = grouping.bootstrap_resampled_error(self.uvp, time_avg=True,
+                                                 Nsamples=100, normal_std=True,
+                                                 robust_std=False, verbose=False)
+
+        f6 = plot.delay_spectrum(uvp_avg, uvp_avg.get_blpairs(), spw=0,
+                                pol='xx', average_blpairs=False, average_times=False,
+                                component='real', error='normal_std', lines=False,
+                                markers=True)
+        plt.close(f6)
 
     def test_plot_cosmo(self):
         """