Repository: lanl/fiducia (language: python)

Issue:
Initial release of code
Release merge fiducia code into repo.


Discussion / hints:

Fix patch:
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..0ed0081
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,36 @@
+© 2020. Triad National Security, LLC. All rights reserved.
+This program was produced under U.S. Government contract 89233218CNA000001 
+for Los Alamos National Laboratory (LANL), which is operated by Triad National 
+Security, LLC for the U.S. Department of Energy/National Nuclear 
+Security Administration. All rights in the program are reserved by Triad 
+National Security, LLC, and the U.S. Department of Energy/National 
+Nuclear Security Administration. The Government is granted for itself and 
+others acting on its behalf a nonexclusive, paid-up, irrevocable 
+worldwide license in this material to reproduce, prepare derivative 
+works, distribute copies to the public, perform publicly and display 
+publicly, and to permit others to do so.
+
+This program is open source under the BSD-3 License.
+Redistribution and use in source and binary forms, with or without 
+modification, are permitted provided that the following conditions are met:
+1. Redistributions of source code must retain the above copyright notice, 
+this list of conditions and the following disclaimer.
+
+2.Redistributions in binary form must reproduce the above copyright notice, 
+this list of conditions and the following disclaimer in the documentation 
+and/or other materials provided with the distribution.
+
+3.Neither the name of the copyright holder nor the names of its contributors 
+may be used to endorse or promote products derived from this software 
+without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE 
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
\ No newline at end of file
diff --git a/fiducia/.coverage b/fiducia/.coverage
new file mode 100644
index 0000000..4f43693
--- /dev/null
+++ b/fiducia/.coverage
@@ -0,0 +1,1 @@
+!coverage.py: This is a private format, don't read it directly!{"lines":{"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/examples/plot_example_test.py":[6,9,10,11,18,19,25,26,27,28,29],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/__init__.py":[7],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/__init__.py":[7],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_cspline.py":[9,12,13,14,15,18,19,20,21,23,26,27,35,45,54,57,58,84,94,105,106,107,127,136,145,146,147,155,164,167,168,184,195,198,199,250,259,29,30,31,32,33,39,40,41,43,47,48,49,50,52,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,86,87,88,89,90,92,96,97,98,99,100,102,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,129,130,131,132,134,138,139,140,141,143,148,149,150,151,152,153,157,158,159,160,162,170,171,172,175,176,178,179,180,181,182,186,187,188,189,190,191,193,201,202,203,204,205,206,207,208,210,211,212,213,214,215,216,217,218,219,220,221,223,224,225,226,227,228,229,230,232,233,234,235,236,237,238,239,241,242,243,244,245,246,247,248,252,253,254,255,257,261,262,263,264,266],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/cspline.py":[7,10,11,12,13,14,16,17,21,22,23,24,25,26,27,28,29,30,31,32,33,37,81,116,159,201,232,263,299,363,412,456,529,567,664,742,77,78,112,113,145,146,147,148,149,150,152,153,154,155,156,188,189,190,191,192,193,195,196,197,198,223,224,225,226,227,228,229,254,255,256,257,258,259,260,291,292,293,295,296,340,341,346,349,351,354,355,356,357,358,359,360,404,406,408,409,452,453],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/pltDefaults.py":[9,12,13,14,15,20,21,25,26,27,28,29,30,34,36,37,38,39,40,41,44,45,46,47,48,65,124,212,215,218,222,223],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/stats.py":[8,10,14,15,16,17,18,23,59,117,181,274,101,102,103,104,107,109,110,112,113,114,106,320,321,322,325,328,330,332,339,346,347,349,351,358,360,334,337,352,356,341,344],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/loader.py":[9,12,13,14,15,19,20,25,26,27,28,29,30,31,35,65,104,143,203,247,308,371,415,176,178,180,129,130,131,132,134,135,136,137,139,140,183,184,186,190,191,193,194,196,198,200],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/misc.py":[9,12,13,14,18,19,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,46,48,85,79,80,81],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_loader.py":[9,12,13,16],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_main.py":[9,12,13,16,18,19,29,21,22,23,24,25,26,27,31,32,33,34,35,36,37,38,40],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/main.py":[47,50,51,52,53,54,57,61,62,63,64,65,66,71,72,73,74,75,76,80,104,176,235,429,572,158,159,161,163,164,165,166,168,169,170,171,172,173],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/rawProcess.py":[13,16,17,18,19,23,24,25,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,55,90,125,158,200,250,308,380,442,473,529,577,616,679,772,891,956,1000,1059,1165,1239,1293,1294,1295,1469,1544],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/visualization.py":[9,12,13,16,21,22,23,24,31,94,150,193,222,259],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/response.py":[10,13,18,23,24,25],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/error.py":[9,12,13,14,15,16,18,19,20,21,22,25,26,27,28,29,30,37,143,144,145,146,226,300,392,481],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_misc.py":[9,12,13,16,19,23,24,25],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_response.py":[9,12,13,16],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_stats.py":[8,11,12,15,17,18,25,31,39,20,21,23,27,29,33,34,35,36,42,43],"/home/pkozlowski/LANL_postdoc/codes/python/fiducia/fiducia/tests/test_visualization.py":[9,12,13,16]}}
\ No newline at end of file
diff --git a/fiducia/__init__.py b/fiducia/__init__.py
new file mode 100644
index 0000000..f97aad1
--- /dev/null
+++ b/fiducia/__init__.py
@@ -0,0 +1,8 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 25 13:14:52 2019
+
+@author: Pawel M. Kozlowski
+"""
+
diff --git a/fiducia/cspline.py b/fiducia/cspline.py
new file mode 100644
index 0000000..dcb67c6
--- /dev/null
+++ b/fiducia/cspline.py
@@ -0,0 +1,836 @@
+"""
+Created on Fri Mar  8 09:41:36 2019
+
+Functions for working with cubic spline equation in matrix form.
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import scipy.sparse as sparse
+from scipy import integrate
+import matplotlib.pyplot as plt
+import xarray as xr
+# custom modules
+import fiducia.pltDefaults
+from fiducia.stats import dotVariance
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["splineCoords",
+           "splineCoordsInv",
+           "yCoeffArr",
+           "dCoeffArr",
+           "dToyArr",
+           "responseInterp",
+           "yChiCoeffArr",
+           "yChiCoeffArrEnergies",
+           "fancyTrapz2",
+           "segmentsArr",
+           "detectorArr",
+           "knotSolve",
+           "reconstructSpectrum",
+           ]
+
+
+def splineCoords(energy, energyStart, energyEnd):
+    r"""
+    Convert photon energy value into normalized coordinates for a particular
+    spline region.
+    
+    Parameters
+    ----------
+    energy: float, numpy.ndarray
+        Energy value to be converted into normalized spline coordinate.
+        
+    energyStart: float
+        Lower bound energy for the spline region based on knot points.
+        
+    energyEnd: float
+        Upper bound energy for the spline region based on knot points.
+    
+    normCoord: float
+        Return value of energy converted into normalized spline coordinates.
+    
+    Returns
+    -------
+    normCoord: float, numpy.ndarray
+        Normalized energy coordinate(s).
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+#    if energyStart > energyEnd:
+#        raise ValueError(f"Lower bound energy {energyStart} must be less than "
+#                        f"upper bound energy {energyEnd}.")
+#    if np.all(energy < energyStart) or np.all(energy > energyEnd):
+#        raise ValueError(f"Given energy value {energy} is outside spline "
+#                        f"coordinate bounds {energyStart} to {energyEnd}.")
+    normCoord = (energy - energyStart) / (energyEnd - energyStart)
+    return normCoord
+
+
+def splineCoordsInv(energyNorm, energyStart, energyEnd):
+    r"""
+    Given a normalized energy value and the bounds of a spline segment, 
+    return the un-normalized photon energy value. This is the inverse of
+    splineCoords().
+    
+    Parameters
+    ----------
+    energyNorm: float, numpy.ndarray
+        Normalized photon energy
+        
+    energyStart: float
+        Lower bound energy for the spline region based on knot points.
+        
+    energyEnd: float
+        Upper bound energy for the spline region based on knot points.
+    
+    Returns
+    -------
+    energy: float, numpy.ndarray
+        Absolute photon energy (un-normalized).
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    energy = (energyEnd - energyStart) * energyNorm + energyStart
+    return energy
+
+
+def yCoeffArr(energyNorm, chLen):
+    r"""
+    Returns the matrix M_y(t) for a given value of t in:
+        
+    .. math::
+        Y_i(t) = M_y(t) y_i + M_D(t) D_i
+        
+    Parameters
+    ----------
+    energyNorm: float
+        normalized photon energy for a spline section.
+        
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+    
+    Returns
+    -------
+    mArr: scipy.sparse.lil.lil_matrix
+        Sparse matrix :math:`M_y(t)`.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    coeff1 = 1 - 3 * energyNorm ** 2 + 2 * energyNorm **3 
+    coeff2 = 3 * energyNorm ** 2 - 2 * energyNorm ** 3
+    mArr = sparse.diags([coeff1, coeff2],
+                        [0,1],
+                        shape=(chLen, chLen + 1),
+                        format='lil')
+    #First and last row elements in y array are different (added by DHB 3/25/19)
+    mArr[0,0] = 1 - energyNorm ** 3
+    mArr[0,1] = energyNorm ** 3
+    mArr[chLen-1,chLen-1] = 1 - energyNorm ** 3
+    mArr[chLen-1,chLen] = energyNorm ** 3
+    return mArr.tocsc()
+
+
+def dCoeffArr(energyNorm, chLen):
+    r"""
+    Returns the matrix M_D(t) for a given value of t in:
+        
+    .. math::
+        Y_i(t) = M_y(t) y_i + M_D(t) D_i
+    
+    Parameters
+    ----------
+    energyNorm: float
+        normalized photon energy for a spline section.
+        
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+    
+    Returns
+    -------
+    dArr: scipy.sparse.lil.lil_matrix
+        Sparse matrix :math:`M_D(t)`.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    coeff1 = energyNorm - 2 * energyNorm ** 2 + energyNorm ** 3
+    coeff2 = -1 * energyNorm ** 2 + energyNorm ** 3
+    dArr = sparse.diags([coeff1, coeff2],
+                        [0,1],
+                        shape=(chLen, chLen + 1), 
+                        format='lil')
+    #First and last row elements in D array are different (added by DHB 3/25/19)
+    dArr[0,0] = energyNorm - energyNorm ** 3
+    dArr[chLen-1,chLen-1] = energyNorm - energyNorm ** 2
+    dArr[chLen-1,chLen] = energyNorm ** 2 - energyNorm ** 3
+    return dArr.tocsc()
+
+
+def __chi1Arr__(chLen):
+    r"""
+    Parameters
+    ----------
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+    
+    Returns
+    -------
+    chi1: scipy.sparse.lil.lil_matrix
+        Sparse matrix :math:`\chi_1`.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # chi1array in Dan Barnak's notation = TDy'' in Jim's paper 
+    chi1 = sparse.diags([1, 4, 1], 
+                        [-1,0,1],
+                        shape=(chLen + 1, chLen + 1),
+                        format='lil')
+    chi1[0, 0] = 2
+    chi1[-1, -1] = 2
+    return chi1.tocsc()
+
+
+def __chi3Arr__(chLen):
+    r"""
+    Parameters
+    ----------
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+        
+    Returns
+    -------
+    chi3: scipy.sparse.lil.lil_matrix
+        Sparse matrix :math:`\chi_3`.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # chi2array in Dan Barnak's notation = TDy in Jim's paper
+    chi3 = sparse.diags([-1, 0, 1],
+                        [-1,0,1],
+                        shape=(chLen + 1, chLen + 1),
+                        format='lil')
+    chi3[0, 0] = -1
+    chi3[-1, -1] = 1
+    return chi3.tocsc()
+
+
+def dToyArr(chLen):
+    r"""
+    Construct matrix for converting from :math:`D_i` to :math:`y_i` vector.
+    
+    Parameters
+    ----------
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+        
+    Returns
+    -------
+    diToyi: numpy.ndarray
+        Matrix for converting from :math:`D_i` to :math:`y_i` vector.
+        
+    Notes
+    -----
+    The matrix is given by:
+        
+    .. math::
+        D_i = 3 \chi_1^{-1} \chi_3 y_i
+    
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    chi1Arr = __chi1Arr__(chLen)
+    chi3Arr = __chi3Arr__(chLen)
+    chi1ArrInv = np.linalg.inv(chi1Arr.toarray())
+    # constructing array for converting values from Di to yi
+    diToyi = 3 * np.dot(chi1ArrInv, chi3Arr.toarray())
+    return diToyi
+
+
+def responseInterp(energyNorm, energyMin, energyMax, responseFrame, channels):
+    r"""
+    Given a DANTE detector response as a function of energy, convert the
+    response to normalized photon energy, t, over a given spline segment, and
+    return interpolated response values for a given value of t. Returns an
+    array of interpolated responses corresponding to the number of channels.
+    
+    Parameters
+    ----------
+    energyNorm: float, numpy.ndarray
+        normalized photon energy
+        
+    energyMin: float
+        Lower bound photon energy of the spline segment over which we are
+        normalizing.
+        
+    eneryMax: float
+        Upper bound photon energy of the spline segment over which we are
+        normalizing.
+        
+    responseFrame: pandas.core.frame.DataFrame
+        DANTE channel responses as a function of photon energy (not normalized).
+
+    channels: numpy.ndarray
+        numpy array of DANTE channel numbers.
+    
+    Returns
+    -------
+    responsesInterpd: numpy.ndarray
+        Returns a matrix of (energyNorms, channels) of response functions.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+
+    """
+    chLen = len(channels)
+    if np.shape(energyNorm):
+        # if energyNorm is a vector, then get the length of that vector
+        energyLen = len(energyNorm)
+    else:
+        # otherwise the length is just 1
+        energyLen = 1
+    # construct an array of interpolated response functions based on
+    # number of photon energy points we want and number of DANTE channels
+    responsesInterpd = np.zeros((energyLen, chLen))
+    # fetch the original array of energy values (not normalized)
+    energyArr = responseFrame['Energy(eV)']
+    # converting normalized energy values to un-normalized energy values,
+    # since our response functions are in terms of absolute photon energy
+    energyReg = splineCoordsInv(energyNorm, energyMin, energyMax)
+    for idx, channel in enumerate(channels):
+        responseArr = responseFrame[channel]
+        responsesInterpd[:,idx] = np.interp(x=energyReg,
+                                            xp=energyArr,
+                                            fp=responseArr)
+    return responsesInterpd
+
+
+def yChiCoeffArr(energyNorm, chLen, dToY):
+    r"""
+    This is the matrix corresponding to:
+        
+    .. math::
+        M_y(t) + 3 M_D(t) \chi_1^{-1} \chi_3
+    
+    which describes the cubic spline interpolation of the x-ray
+    spectrum.
+    
+    Parameters
+    ----------
+    energyNorm: float
+        normalized photon energy
+        
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+        
+    dToY: numpy.ndarray
+        Matrix for converting from :math:`D_i` to :math:`y_i` values in
+        cubic spline interpolation. See dToyArr().
+    
+    Returns
+    -------
+    yChiArr: numpy.ndarray
+        Returns a 2D matrix for a particular value of energyNorm.
+        
+    Notes
+    -----
+    The matrix is given by:
+        
+    .. math::
+        M_y(t) + 3 M_D(t) \chi_1^{-1} \chi_3
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # M_y(t)
+    yCoeff = yCoeffArr(energyNorm, chLen)
+    # M_D(t)
+    dCoeff = dCoeffArr(energyNorm, chLen)
+    # folded x-ray spectrum and detector response array
+    yChiArr = yCoeff.toarray() + np.dot(dCoeff.toarray(), dToY)
+    return yChiArr
+
+
+def yChiCoeffArrEnergies(energyNorms, chLen, dToY):
+    r"""
+    This is the matrix corresponding to:
+        
+    .. math::
+        M_y(t) + 3 M_D(t) \chi_1^{-1} \chi_3
+    
+    which describes the cubic spline interpolation of the x-ray
+    spectrum.
+    
+    energyNorms: numpy.ndarray
+        Vector of normalized photon energies.
+        
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+        
+    dToY: numpy.ndarray
+        Matrix for converting from :math:`D_i` to :math:`y_i` values in cubic
+        spline interpolation. See dToyArr().
+    
+    Returns
+    -------
+    yChiArrEnergies: numpy.ndarray
+        Returns a 3D matrix composed of a series of 2D yChiCoeff matrices
+        corresponding to the given energyNorm values. This matrix is
+        indexed as (energyNorms, knotIndex, knotIndex).
+        
+    Notes
+    -----
+    The matrix is given by:
+        
+    .. math::
+        M_y(t) + 3 M_D(t) \chi_1^{-1} \chi_3
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    yChiArrEnergies = np.array([yChiCoeffArr(energyNorm, chLen, dToY) for energyNorm in energyNorms])
+    return yChiArrEnergies
+    
+
+def fancyTrapz2(energyNorms, yChis, segments, responseFrame, channels):
+    r"""
+    Trap rule integration of the folding between our :math:`M_{y \chi}` matrix
+    and response function matrix, with respect to normalized photon energy,
+    for each channel. The result should be a matrix with shape
+    (`len(channels)`, `len(segments)`, `len(knotIndex)`).
+    
+    Parameters
+    ----------
+    energyNorms: numpy.ndarray
+        1D array of normalized photon energies
+        
+    yChis: numpy.ndarray
+        3D array of :math:`M_{y \chi}` values corresponding to
+        (energyNorms, segments, knotIndex).
+        
+    responses: numpy.ndarray
+        2D array of DANTE channel response functions corresponding to
+        (energyNorms, channels).
+        
+    channels: numpy.ndarray
+        Array of DANTE channel numbers.      
+        
+    Returns
+    -------
+    integArr : xarray.DataArray
+        A matrix containing the folded integration of the :math:`M_{y \chi}` 
+        matrix and response function matrix, with respect to normalized photon 
+        energy. Has shape (`len(channels)`, `len(segments)`, `len(knotIndex)`).
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    shape = np.shape(yChis)
+    segmentsLen = np.shape(segments)[0]
+    knotsLen = shape[2]
+    chLen = segmentsLen
+    # initialize integArr for storing photon energy integrated values
+    integArr = xr.DataArray(np.zeros((chLen, segmentsLen, knotsLen)),
+                            dims=['channel', 'segment', 'knot_point'],
+                            coords={'channel':channels})
+    # loop over relevant DANTE channels for analysis
+    for channelIdx in np.arange(integArr['channel'].size):
+        # loop over photon energy segments (between knot points)
+#        print(f"Calculating for channel {channels[channelIdx]}.")
+        for segmentNum, segment in enumerate(segments):
+            energyMin, energyMax = segment
+#            print(f"Calculating for segment {segments[segmentNum]}.")
+            # loop over knot points
+            for knotNum in np.arange(integArr['knot_point'].size):
+                # multiplication of response by spline matrix
+                responses = responseInterp(energyNorms,
+                                           energyMin,
+                                           energyMax,
+                                           responseFrame,
+                                           channels)
+#                plt.plot(energyNorms, responses[:, channelIdx])
+#                plt.title(f"Response channel {channels[channelIdx]}, segment {segmentNum}, knot {knotNum}")
+#                plt.show()
+                multArr = yChis[:, segmentNum, knotNum] * responses[:, channelIdx]
+#                integArr[channelIdx, segmentNum, knotNum] = np.trapz(y=multArr,
+#                                                                     x=energyNorms)
+                integVal = (energyMax - energyMin) * integrate.simps(y=multArr, x=energyNorms)
+                integArr[channelIdx, segmentNum, knotNum] = integVal
+    return integArr
+
+
+def segmentsArr(knots):
+    r"""
+    Returns the bounds of each spline segment, given the spline knot points.
+    
+    Returns an array of tuples of (energyMin, energyMax) describing the bounds
+    of each spline segment, given an array of spline knots (photon energies
+    corresponding to K-edges).
+    
+    Parameters
+    ----------
+    knots: numpy.ndarray
+        numpy array of photon energies describing positions of spline knots.
+        
+    Returns
+    -------
+    segments : numpy.ndarray
+        A 1D array of tuples of (energyMin, energyMax), corresponding to the 
+        bounds of each spline segment. Has length of 'len(knots) - 1'.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # initializing array of segments
+    segments = np.zeros(len(knots) - 1, dtype=object)
+    for idx, knot in enumerate(knots):
+        if idx == len(segments):
+            # skipping the last knot since there are 1 fewer segments than knots
+            continue
+        segments[idx] = knots[idx], knots[idx + 1]
+    return segments
+    
+
+def detectorArr(channels, knots, responseFrame, boundary="y0", npts=1000):
+    r"""
+    Matrix representing the spectrally integrated folding of the detector
+    response with a cubic spline interpolation of the x-ray spectrum. This
+    is applied to the measured DANTE channel signals to recover knot points
+    :math:`y_i` of the cubic spline, which can then be used to reconstruct the
+    inferred x-ray spectrum.
+    
+    Parameters
+    ----------
+    channels: numpy.ndarray
+        Array of DANTE channel numbers.
+    
+    knots: numpy.ndarray
+        Array of photon energies describing positions of spline knots.
+    
+    responseFrame: pandas.core.frame.DataFrame
+        DANTE channel responses as a function of photon energy (not 
+        normalized).
+    
+    boundary: str, optional
+        Choose whether yGuess corresponds to :math:`y_0` (lowest photon 
+        energy) or :math:`y_{n+1}` (highest photon energy) boundary condition.
+        This should correspond to the photon energy value given in knots.
+        Options are 'y0' or 'yn+1'. Default 'y0'.
+        
+    npts: int
+        Number of points used in computing the integral
+        
+    Returns
+    -------
+    detArr: numpy.ndarray
+        Matrix representing the spectrally integrated folding of the detector
+        response with a cubic spline interpolation of the x-ray spectrum.
+        2D array of channels and knot points of shape (n, n).
+        
+    Notes
+    -----
+    For each spline segment we have:
+        
+    .. math::
+        M_{stuff} = \int_0^{1} (M_y(t) + 3 M_D(t) \chi_1^{-1} \chi_3) R_d(t) dt
+    
+    Each spline is then summed to form the full detector matrix for recovering
+    the knot points.
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # number of DANTE channels where we have useful measurements
+    chLen = len(channels)
+    # initialize normalized energies array
+    # array of normalized energies over which we do the integral
+    energyNorms = np.linspace(0, 1, num=npts)
+    # producing segments from knots
+    segments = segmentsArr(knots)
+    # calculating array for converting from values of D_i to y_i. This
+    # is an optimization as this array is constant!
+    dToY = dToyArr(chLen)
+    # M_{y \chi} coefficients array corresponding to given normalized
+    # energies. Array shape is (energyNorms, segments, knotIndex).
+    yChis = yChiCoeffArrEnergies(energyNorms, chLen, dToY)
+
+    integFoldArr = fancyTrapz2(energyNorms,
+                               yChis,
+                               segments,
+                               responseFrame,
+                               channels)
+    # sum along segment axis, as each segment must contribute to the
+    # overall signal.
+    #detArr = np.sum(integFoldArr, axis=1)
+    detArr = integFoldArr.sum(dim="segment")
+    detArr.attrs['boundary'] = boundary
+    if boundary == "y0":
+        # extracting column corresponding to y0
+        detArrBoundaryCol = detArr.isel(knot_point=0)
+        detArr = detArr.isel(knot_point=slice(1, None))
+    elif boundary == "yn+1":
+        # extracting column corresponding to y_{n+1}
+        detArrBoundaryCol = detArr.isel(knot_point=-1)
+        detArr = detArr.isel(knot_point=slice(None, -1))
+    else:
+        raise Exception(f"No method found for boundary {boundary}.")
+    return detArr, detArrBoundaryCol
+
+
+def knotSolve(signals,
+              detArr,
+              detArrBoundaryCol,
+              detArrVarianceBoundaryCol,
+              detArrInv,
+              stdDetArrInv,
+              signalsUncertainty=None,
+              yGuess=1e-77,
+              npts=1000):
+    r"""
+    Get knot points :math:`y_i` from measured DANTE signals :math:`S_d`.
+    
+    Parameters
+    ----------
+    signals: numpy.ndarray
+        numpy array of DANTE measured signal for each channel at a particular
+        point in time.
+        
+    detArr : numpy.ndarray
+        Matrix representing the spectrally integrated folding of the detector
+        response with a cubic spline interpolation of the x-ray spectrum.
+        2D array of channels and knot points of shape (n, n).
+    
+    detArrBoundaryCol : xarray.DataArray
+        Column of cublic spline matrix corresponding to the knots at the 
+        boundary chosen with `boundary`.
+    
+    detArrVarianceBoundaryCol: xarray.DataArray
+        Column of variances in the cublic spline matrix corresponding to the 
+        knots at the boundary chosen with `boundary`.
+        
+    detArrInv : xarray.DataArray
+        Inversion of detArr, with the column corresponding to boundary removed 
+        so detArr is invertible.
+    
+    stdDetArrInv: xarray.DataArray
+        Array of the standard deviation of each element in detArrInv based on
+        variance using the 'responseUncertaintyFrame' propagated with Monte Carlo. 
+
+    signalsUncertainty: xarray.DataArray, optional
+        numpy array of the uncertainty of the DANTE measured signal for each 
+        channel at a particular point in time. The default is None.
+
+    yGuess: float, optional
+        Guess for position of boundary knot point. Default is 1e-77.
+    
+    npts: int, optional
+        Number of points used in computing the integral. Default is 1000.
+        
+    Returns
+    -------
+    knotsY : numpy.ndarray
+        Array of knot point intensity values with yGuess appended.
+        
+    knotsYVariance: numpy.ndarray
+        Array with each element corresponding to the variance of the same
+        element in 'knotsY'.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """ 
+    #subtract boundary col from signals
+    signalsy0 = signals - yGuess * detArrBoundaryCol
+    signalsy0Variance = signalsUncertainty**2  + yGuess**2 * detArrVarianceBoundaryCol
+    signalsy0Uncertainty = np.sqrt(signalsy0Variance)
+    # applying inverted array to signals to recover knot points
+    #huge difference between xarray.dot and np.dot. See issue in Gitlab
+    knotsY = np.dot(detArrInv, signalsy0)
+    #knotsYother = detArrInv.dot(signalsy0, {subscripts:'i,i'})
+    #print("knot diff", knotsYother-knotsY)
+    knotsYVariance = dotVariance(detArrInv, signalsy0, stdDetArrInv, signalsy0Uncertainty)
+    return knotsY, knotsYVariance
+
+
+def reconstructSpectrum(chLen,
+                        knots,
+                        knotsY,
+                        knotsYUncertainty=None,
+                        npts=1000,
+                        plot=False):
+    r"""
+    Reconstruct the inferred DANTE spectrum given the knot points
+    :math:`y_i` obtained from knotSolve().
+    
+    Parameters
+    ----------
+    chLen: int
+        Number of DANTE channels (equal to number of spline knots).
+        
+    knots: list, numpy.ndarray
+        List or array of knot point photon energy value. See knotFind().
+        
+    knotsY: numpy.ndarray
+        Array of knot point intensity values with yGuess appended.
+        See knotSolve() and analyzeSpectrum().
+
+    knotsYUncertainty: numpy.ndarray
+        Array of knot point intensity uncertainty values with yGuess appended.
+        See knotSolve() and analyzeSpectrum(). The default is None.
+    
+    npts: int
+        Number of points used in computing the integral. The default is 1000.
+        
+    plot: Bool
+        Flag for plotting unfolded spectrum. The default is False.
+        
+    Returns
+    -------
+    photonEnergies: numpy.ndarray
+        Photon energy axis of unfolded spectrum.
+        
+    intensities: numpy.ndarray
+        Spectral intensity axis of unfolded spectrum.
+        
+    intensitiesVariance: numpy.ndarray
+        Uncertaitny (1 :math:`\sigma`) on spectral intensity values.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    #check for no uncertainty provided.
+    if knotsYUncertainty is None:
+        #give them uncertainty of 0
+        knotsYUncertainty = np.zeros(knotsY.shape)
+    # initialize normalized energies array
+    energyNorms = np.linspace(0, 1, num=npts)
+    # producing segments from knots
+    segments = segmentsArr(knots)
+    # getting cubic spline matrices for reconstructing spline
+    dToY = dToyArr(chLen)
+    yChis = yChiCoeffArrEnergies(energyNorms, chLen, dToY)
+    #set a constant simple uncertainty for the yChis.
+    #TOOD replace with something physical
+    yChisUncertainty = np.zeros(yChis.shape)
+    # recovering the spectral values
+    spectrum = np.dot(yChis, knotsY)
+    spectrumVariance = dotVariance(yChis,
+                                   knotsY,
+                                   yChisUncertainty,
+                                   knotsYUncertainty)
+    # print(spectrumVariance)
+    
+    # initialize arrays for holding photon energies and corresponding
+    # intensities of stitched spectrum.
+    photonEnergies = np.array([])
+    intensities = np.array([])
+    intensitiesVariance = np.array([])
+    # loop over segments to reconstruct spectrum
+    for segNum, segment in enumerate(segments):
+        energyMin, energyMax = segment
+        energyRegs = splineCoordsInv(energyNorms, energyMin, energyMax)
+        # appending each segment to form a single array of spectral
+        # values. We need to cut the last element from each array as the
+        # last element overlaps the first element of the next segment.
+        photonEnergies = np.append(photonEnergies, energyRegs[:-1])
+        intensities = np.append(intensities, spectrum[:-1, segNum])
+        intensitiesVariance = np.append(intensitiesVariance,
+                                        spectrumVariance[:-1, segNum])
+        if plot:
+            plt.plot(energyRegs, spectrum[:, segNum])
+            
+    if plot:
+        plt.ylabel('Spectrum (GW/sr/eV)')
+        plt.xlabel('Photon Energy (eV)')
+        plt.title('Spectrum from cubic spline knots')
+        plt.show()
+    return photonEnergies, intensities, intensitiesVariance
\ No newline at end of file
diff --git a/fiducia/error.py b/fiducia/error.py
new file mode 100644
index 0000000..757c549
--- /dev/null
+++ b/fiducia/error.py
@@ -0,0 +1,597 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tues June 16 13:48:21 2020
+
+Utilities for calculating response uncertainty
+
+@author: Myles Brophy
+"""
+
+# python modules
+import numpy as np
+import xarray as xr
+import pandas as pd
+import matplotlib.pyplot as plt
+import datetime
+#custom modules
+from fiducia.cspline import splineCoordsInv, segmentsArr, dToyArr, yChiCoeffArrEnergies, responseInterp, detectorArr
+from fiducia.misc import areDataFramesCompatible
+from fiducia.stats import trapzVariance, gradientVariance, interpVariance
+from fiducia.response import knotFind
+from fiducia.loader import loadResponses, loadResponseUncertainty
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["detectorErrMC",
+           "knotVarianceFind",
+           "responseInterpVariance",
+           "fancyTrapz2Variance",
+           "detectorArrVariance",
+           "detectorUncertainty"
+           ]
+
+
+
+    
+def detectorErrMC(detArr, detArrVariance, samples=10000,
+                  boundary="y0", MChistogram=False):
+    r"""
+    Monte Carlo simulation and statistics to determine cubic spline uncertainty.
+    
+    Calculate the cubic spline matrix uncertainty using a Monte Carlo simulation
+    and statistics on the MC's output.
+    
+    Parameters
+    ----------
+    detArr: numpy.ndarray
+        Matrix representing the spectrally integrated folding of the detector
+        response with a cubic spline interpolation of the x-ray spectrum. See
+        'cspline.detectorArr()'.
+    
+    detArrVariance:  numpy.ndarray
+        A DataFrame containing the uncertainty for each Dante channel for the
+        photon energy range that detArr spans.
+    
+    samples : int, optional
+        Number of MCs amples to run. Default is `10000`.
+    
+    boundary : str, optional
+        Choose whether yGuess corresponds to :math:`y_0` (lowest photon energy) or
+        :math:`y_{n+1}` (highest photon energy) boundary condition. This should
+        correspond to the photon energy value given in knots. Options are `y0`
+        or `yn+1`. Default is 'y0'.
+    
+    MChistogram : bool, optional
+        Plot histograms corresponding to each variant of detArr generated with
+        Monte Carlo uncertainty propagation.
+        Default is `False`.
+    
+    Returns
+    -------
+    stdErrorMatrix : numpy.ndarray
+         A numpy.ndarray with the standard deviation of the inverted matrices
+         generated using random weights based on the channel uncertainty.
+        
+    Raises
+    ------
+    Exception
+        If `boundary` doesn't equal `y0` or `yn+1`.
+    ValueError
+        If the shapes of `detArr` and `detUncertaintyArr` aren't equal.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    
+    if detArr.shape != detArrVariance.shape:
+        raise ValueError(f"Shape of the cubic spline matrix ({detArr.shape}) is not equal to"+
+                         f" the shape of the detector uncertainty matrix ({detArrVariance.shape})")
+    #square root the uncertainty matrix because it is sigma^2 and we want sigma
+    detArrUncertainty = np.sqrt(detArrVariance)
+    #make a 3d matrix, where each 2d slice is an error matrix generated by MC
+    print("Generating "+str(samples) + " MC samples...")
+    error = xr.DataArray(np.zeros((samples,)+detArr.shape), dims=['sample', 'channel', 'knot_point'])
+    for s in range(samples):
+        error[s] = detArr + np.random.normal(loc=0, scale=1, size=detArr.shape)*detArrUncertainty
+
+    print("Finished generating MC samples.")
+    
+    detArrInv = xr.DataArray(np.linalg.inv(detArr), dims=['channel', 'knot_point'], attrs={'boundary':boundary})  
+    invError = xr.DataArray(np.zeros(error.shape), dims=['sample', 'channel', 'knot_point'])
+    
+    # #invert each 2d slice of the error matrix
+    for sample in range(samples):
+        invErrSlice = np.linalg.inv(error[sample])
+        invError[sample] = invErrSlice
+            
+    # #find std of each element compared to other samples
+    stdDetArrInv = invError.std(dim="sample")
+    stdDetArrInv.attrs['boundary'] = boundary
+    
+    if MChistogram:
+        for ch in invError['channel']:
+            for knot in invError['knot_point']:
+                # print(len(invError.sel(channel=ch, knot_point=knot).values))
+                plt.hist(invError.sel(channel=ch, knot_point=knot).values,
+                         samples,
+                         alpha=1,
+                         label=f'Channel {ch.values} knot {knot.values}')
+                plt.axvline(detArrInv.sel(channel=ch, knot_point=knot).values,
+                            label='DetArrInv',
+                            color='red')
+                plt.axvline(detArrInv.sel(channel=ch,knot_point=knot).values + stdDetArrInv.sel(channel=ch, knot_point=knot).values,
+                            label='DetArrInv+sigma',
+                            color='purple')
+                plt.axvline(detArrInv.sel(channel=ch, knot_point=knot).values - stdDetArrInv.sel(channel=ch, knot_point=knot).values,
+                            label='DetArrInv-sigma',
+                            color='green')
+                plt.legend(loc='upper right')
+                plt.title("Detector MC")
+                plt.show()
+                
+    return detArrInv, stdDetArrInv
+    
+
+#TODO check to see if this has any effect on final error bars. Pawel thinks not.
+def knotVarianceFind(channels,
+                     responseUncertaintyFrame=None,
+                     forceKnot=np.array([]), 
+                     knotBoundaryY=1e-77,
+                     boundary='y0'):
+    r"""
+    Modification of response.knotFind()
+
+    Parameters
+    ----------
+    channels: numpy.ndarray
+        Array of DANTE channel numbers. 
+
+    responseUncertaintyFrame: pandas.core.frame.DataFrame, optional
+        DataFrame holding percent uncertainties of DANTE channel responses as
+        a function of photon energy (not normalized). The default is `None`.
+
+    forceKnot : TYPE, optional
+        DESCRIPTION. The default is `np.array([])`.
+    
+    knotBoundaryY : float, optional
+        Guess for position of y_0 or y_{n+1} knot point. Default is 1e-77.
+    
+    boundary: str, optional
+        Choose whether yGuess corresponds to :math:`y_0` (lowest photon
+        energy) or :math:`y_{n+1}` (highest photon energy) boundary condition.
+        This should correspond to the photon energy value given in knots.
+        Options are `y0` or `yn+1`. Default `y0`.
+
+    Returns
+    -------
+     knotUncertainty: numpy.ndarray
+         An array of uncertainty in knot points, with each element
+         corresponding to a channel or boundary condition. 
+         See :func:`response.knotFind`.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #modified from response.knotFind
+    knotsUncertainty = np.zeros(len(channels))
+    
+    #check if uncertainty frame is passed.
+    if responseUncertaintyFrame is None:
+        #return ndarray of length channels + 1 for boundary condition
+        return  np.zeros(len(channels)+1)
+    
+    for idx, channel in enumerate(channels):
+        # if forceKnot isn't an empty list, then we go about forcing the user
+        # provided values.
+        if forceKnot.size != 0:
+            if channel in forceKnot[:,0]:
+                # if the channel is forced by the user then put the user provided
+                # photon energy into the knots array
+                forceIdx = np.where(forceKnot[:,0] == channel)
+                knotsUncertainty[idx] = forceKnot[forceIdx, 1]
+            else:
+                # finding largest negative gradient, which should correpond to the
+                # K-edge of the DANTE channel filter.
+                grad = -gradientVariance(responseUncertaintyFrame[channel])
+                maxIndex = np.argmax(grad)
+                knotsUncertainty[idx] = responseUncertaintyFrame['Energy(eV)'][maxIndex]
+        else:
+            # finding largest negative gradient, which should correpond to the
+            # K-edge of the DANTE channel filter.
+            grad = -gradientVariance(responseUncertaintyFrame[channel])
+            maxIndex = np.argmax(grad)
+            knotsUncertainty[idx] = responseUncertaintyFrame['Energy(eV)'][maxIndex]
+        
+    if boundary == 'y0':
+        knotsUncertaintyAppend = np.append([knotBoundaryY], knotsUncertainty)
+    elif boundary == 'yn+1':
+        knotsUncertaintyAppend = np.append(knotsUncertainty, [knotBoundaryY])
+    else:
+        raise Exception(f"No method found for boundary {boundary}.")
+    return knotsUncertaintyAppend
+
+
+def responseInterpVariance(energyNorm,
+                           energyMin,
+                           energyMax,
+                           responseUncertaintyFrame,
+                           channels):
+    r"""
+    Given a DANTE detector response as a function of energy, convert the
+    response to normalized photon energy, t, over a given spline segment, and
+    return interpolated response values for a given value of t. Returns an
+    array of interpolated responses corresponding to the number of channels.
+    
+    Parameters
+    ----------
+    energyNorm: float, numpy.ndarray
+        normalized photon energy
+        
+    energyMin: float
+        Lower bound photon energy of the spline segment over which we are
+        normalizing.
+        
+    eneryMax: float
+        Upper bound photon energy of the spline segment over which we are
+        normalizing.
+        
+    responseFrame: pandas.core.frame.DataFrame
+        DANTE channel responses as a function of photon energy (not normalized).
+
+    channels: numpy.ndarray
+        numpy array of DANTE channel numbers.
+    
+    Returns
+    -------
+    responsesInterpdVariance: numpy.ndarray
+        Returns a matrix of (energyNorms, channels) of response functions.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    cspline.repsonseInterp
+    
+    Examples
+    --------
+
+    """
+    chLen = len(channels)
+    if np.shape(energyNorm):
+        # if energyNorm is a vector, then get the length of that vector
+        energyLen = len(energyNorm)
+    else:
+        # otherwise the length is just 1
+        energyLen = 1
+    # construct an array of interpolated response functions based on
+    # number of photon energy points we want and number of DANTE channels
+    responsesInterpdVariance = np.zeros((energyLen, chLen))
+    # fetch the original array of energy values (not normalized)
+    energyArr = responseUncertaintyFrame['Energy(eV)']
+    # converting normalized energy values to un-normalized energy values,
+    # since our response functions are in terms of absolute photon energy
+    energyReg = splineCoordsInv(energyNorm, energyMin, energyMax)
+    for idx, channel in enumerate(channels):
+        responseUncertaintyArr = responseUncertaintyFrame[channel]
+        responsesInterpdVariance[:,idx] = interpVariance(x=energyReg,
+                                                         xp=energyArr,
+                                                         fpUnc=responseUncertaintyArr)
+    return responsesInterpdVariance
+
+
+def fancyTrapz2Variance(energyNorms,
+                        yChis,
+                        segments,
+                        responseUncertaintyFrame,
+                        channels,
+                        interpProp=True):
+    r"""
+    Calculate the variance when propogating uncertainties
+    through :func:`fiducia.cspline.fancyTrapz2`.
+
+    Parameters
+    ----------
+    energyNorms : numpy.ndarray
+        Array of normalized energies over which the integral is computed.
+    
+    yChis : numpy.ndarray
+        3D array corresponding to the :math:`M_{y \chi}` coefficients.
+        Array shape corresponds to (`energyNorms`, `chLen`, `dToY`). 
+        See :func:`fiducia.error.detectorArrVariance`
+    
+    segments : numpy.ndarray
+        Array of segments produced by :func:`segmentsArr` with the knots
+   
+    responseUncertaintyFrame : pandas.core.frame.DataFrame
+        DataFrame holding uncertainty percentages of DANTE channel responses
+        as a function of photon energy (not normalized).
+        
+    channels : numpy.ndarray
+        Array of DANTE channel numbers.
+        
+    interpProp : bool, optional
+        Boolean to decide if :func:`error.responseInterpVariance` should be
+        used. If `False, :func:`cspline.responseInterp()` is used, speeding
+        up the calculation. Note that the uncertainty is would not be
+        propagated correctly if `False`. With future optimizations, this
+        option to choose may be removed. Default is `True`.
+    
+    Returns
+    -------
+    integArrVariance : xarray.Dat
+        A matrix containing the folded integration of the :math:`M_{y \chi}` 
+        matrix and response function uncertainty matrix, with respect to
+        normalized photon  energy.
+        Has shape (`len(channels)`, `len(segments)`, `len(knotIndex)`).
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    cspline.fancyTrapz2()
+    
+    Examples
+    --------
+        
+    """
+    shape = np.shape(yChis)
+    segmentsLen = np.shape(segments)[0]
+    knotsLen = shape[2]
+    chLen = segmentsLen
+    # initialize integArr for storing photon energy integrated values
+    integArrVariance = xr.DataArray(np.zeros((chLen, segmentsLen, knotsLen)),
+                                    dims=['channel', 'segment', 'knot_point'],
+                                    coords={'channel':channels})  
+
+    # loop over relevant DANTE channels for analysis
+    for channelIdx in np.arange(chLen):
+        # loop over photon energy segments (between knot points)
+        for segmentNum, segment in enumerate(segments):
+            energyMin, energyMax = segment
+            # loop over knot points
+            for knotNum in np.arange(knotsLen):
+                # multiplication of response by spline matrix
+                #use interpolation propagation by default, but runs slower
+                #might not have much an effect. For now we'll let the user decide
+                #need to optimize reponseInterpVariance
+                if interpProp:
+                    responsesVariance = responseInterpVariance(energyNorms,
+                                                               energyMin,
+                                                               energyMax,
+                                                               responseUncertaintyFrame,
+                                                               channels)
+                    responsesUncertainty = np.sqrt(responsesVariance)
+                else:
+                    responsesVariance = responseInterp(energyNorms,
+                                                       energyMin,
+                                                       energyMax,
+                                                       responseUncertaintyFrame,
+                                                       channels)
+                    
+                multArr = np.abs(yChis[:, segmentNum, knotNum] * responsesUncertainty[:, channelIdx])
+                integVar = (energyMax - energyMin) ** 2 * trapzVariance(multArr,
+                                                                        x=energyNorms)
+                integArrVariance[channelIdx, segmentNum, knotNum] = integVar
+    return integArrVariance
+
+
+def detectorArrVariance(channels, knots, responseUncertaintyFrame, boundary="y0", npts=1000):
+    r"""
+    Propagates uncertanity through :func:`cspline.detectorArr() to find the variance in :math:`M_{int}`.
+    
+    Parameters
+    ----------
+    channels : numpy.ndarray
+        Array of DANTE channel numbers.
+    
+    knots : numpy.ndarray
+        Array of photon energies describing positions of spline knots.
+    
+    responseUncertaintyFrame : pandas.core.frame.DataFrame
+        DataFrame holding uncertainty percentages of DANTE channel responses
+        as a function of photon energy (not normalized).
+        
+    npts : int, optional
+        Number of points used in computing the integral. The default is 1000.
+        
+    Returns
+    -------
+    detArrVariance : xarray.DataArray
+        2D array of channels and knot points uncertainties of shape `(n, n+1)`.
+
+    detArrVarianceBoundaryCol : xarray.DataArray
+        Column of variances in the cublic spline matrix corresponding to the 
+        knots at the boundary chosen with `boundary`.
+        
+    Notes
+    -----
+    Covariances between segments is not currently accounted for. This
+    covariance should be small compared to the other uncertainties, but should
+    be noted.
+    
+    See also
+    --------
+    cspline.detectorArr()
+    
+    Examples
+    --------
+    
+    """
+    # number of DANTE channels where we have useful measurements
+    chLen = len(channels)
+    # initialize normalized energies array
+    # array of normalized energies over which we do the integral
+    energyNorms = np.linspace(0, 1, num=npts)
+    # producing segments from knotsUncertainty. Use the same segmentsArr
+    #because there is no calculation done
+    segments = segmentsArr(knots)
+    # calculating array for converting from values of D_i to y_i. This
+    # is an optimization as this array is constant!
+    dToY = dToyArr(chLen)
+    # M_{y \chi} coefficients array corresponding to given normalized
+    # energies. Array shape is (energyNorms, segments, knotIndex).
+    yChis = yChiCoeffArrEnergies(energyNorms, chLen, dToY)
+    print("running fancytrapz2Variance")
+    integFoldArrVariance = fancyTrapz2Variance(energyNorms,
+                               yChis,
+                               segments,
+                               responseUncertaintyFrame,
+                               channels)
+    print("finished fancytrapz2")
+    # sum along segment axis, as each segment must contribute to the
+    # overall signal.
+    #detArrVariance = np.sum(integFoldArrVariance, axis=1)
+    detArrVariance = integFoldArrVariance.sum(dim="segment")
+    detArrVariance.attrs['boundary'] = boundary
+    if boundary == "y0":
+        # extracting column corresponding to y0
+        detArrVarianceBoundaryCol = detArrVariance.isel(knot_point=0)
+        detArrVariance = detArrVariance.isel(knot_point=slice(1, None))
+    elif boundary == "yn+1":
+        # extracting column corresponding to y_{n+1}
+        detArrVarianceBoundaryCol = detArrVariance.isel(knot_point=-1)
+        detArrVariance = detArrVariance.isel(knot_point=slice(None, -1))
+    else:
+        raise Exception(f"No method found for boundary {boundary}.")
+    return detArrVariance, detArrVarianceBoundaryCol
+
+
+def detectorUncertainty(channels,
+                        responseFile,
+                        responseUncertaintyFile=None,
+                        boundary="y0",
+                        npts=1000,
+                        samples=1000,
+                        MChistogram=False,
+                        saveDataset=True,
+                        csplineDatasetFile=''):
+    r"""
+    Finds the cspline detector matrix, it`s inverse matrix and std matrix
+    using Monte Carlo uncertainty propagation.
+
+    Propagates response uncertainties through 
+
+    Parameters
+    ----------
+    channels: numpy.ndarray
+        Array of DANTE channel numbers.
+        
+    responseFile: str
+        Path to the `.csv` holding DANTE channel responses as a function of 
+        photon energy (not normalized).
+    
+    responseUncertaintyFile: str, optional
+        Path to the `.csv` holding DANTE channel response uncertainties as a 
+        function of photon energy. Uncertainty values provided as percentages.
+
+    boundary: str, optional
+        Choose whether yGuess corresponds to :math:`y_0` (lowest photon
+        energy) or :math:`y_{n+1}` (highest photon energy) boundary condition.
+        This should correspond to the photon energy value given in knots.
+        Options are `y0` or `yn+1`. Default 'y0'.
+    
+    npts: int, optional
+        Number of points used in computing the integral. Default is 1000.
+    
+    samples: int, optional
+        Number of samples to generate during Monte Carlo propagation.
+        See :func:`error.detectorErrMC`. Default is 1000.
+        
+    Returns
+    -------
+    detArr : xarray.DataArray
+        Matrix representing the spectrally integrated folding of the detector
+        response with a cubic spline interpolation of the x-ray spectrum.
+        2D array of channels and knot points of shape (n, n).
+ 
+    detArrBoundaryCol : xarray.DataArray
+        Column of cublic spline matrix corresponding to the knots at the
+        boundary chosen with `boundary`.
+        
+    detArrVarianceBoundaryCol: xarray.DataArray
+        Column of variances in the cublic spline matrix corresponding to the 
+        knots at the boundary chosen with `boundary`.
+    
+    detArrInv : xarray.DataArray
+        Inversion of detArr, with the column corresponding to boundary removed
+        so detArr is invertible.
+   
+    stdDetArrInv : xarray.DataArray
+        Array of the standard deviation of each element in detArrInv based on
+        variance using the `responseUncertaintyFrame` propagated with Monte
+        Carlo. 
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #load response functions
+    responseFrame = loadResponses(channels, responseFile)
+    
+     #create responseUncertainty with 0s if not given.
+    if responseUncertaintyFile is None:
+        responseUncertaintyFrame = responseFrame.copy()
+        cols = list(responseFrame.columns.values)
+        cols.remove('Energy(eV)')
+        for col in cols:
+            responseUncertaintyFrame.loc[:,col] = 0
+    else:
+        responseUncertaintyFrame = loadResponseUncertainty(responseFrame, responseUncertaintyFile)
+
+    knots = knotFind(channels, responseFrame)
+    
+    #knotsVariance = knotVarianceFind(channels, responseUncertaintyFrame)
+    if not areDataFramesCompatible(channels, responseFrame, responseUncertaintyFrame):
+        raise ValueError("Response frame and response uncertainty frame are not compatible. Check formats.")
+    detArr, detArrBoundaryCol = detectorArr(channels,
+                                            knots,
+                                            responseFrame,
+                                            boundary,
+                                            npts)
+    detArrVariance, detArrVarianceBoundaryCol = detectorArrVariance(channels,
+                                                                    knots,
+                                                                    responseUncertaintyFrame,
+                                                                    boundary,
+                                                                    npts)
+    detArrInv, stdDetArrInv = detectorErrMC(detArr,
+                                            detArrVariance,
+                                            samples,
+                                            boundary,
+                                            MChistogram)
+    #combine into Dataset
+    csplineDataset = xr.Dataset(data_vars={'detArr' : detArr,
+                                           'detArrBoundaryCol' : detArrBoundaryCol,
+                                           'detArrVarianceBoundaryCol' : detArrVarianceBoundaryCol,
+                                           'detArrInv' : detArrInv,
+                                           'stdDetArrInv' : stdDetArrInv}, 
+                                attrs={'boundary' : boundary,
+                                       'channels' : np.asarray(channels), 
+                                       'yGuess' : 0,
+                                       'responseFile' : responseFile, 
+                                       'responseUncertaintyFile' : responseUncertaintyFile,
+                                       'generatedDatetime' : str(datetime.datetime.now())})
+    if saveDataset:
+        #save as netCDF as recommended by http://xarray.pydata.org/en/stable/io.html#netcdf
+        if not csplineDatasetFile:
+            csplineDatasetFile = 'csplineDataset_' + str(datetime.datetime.now().date()) + '.nc'
+        csplineDataset.to_netcdf(csplineDatasetFile)
+    return csplineDataset 
diff --git a/fiducia/examples/README.txt b/fiducia/examples/README.txt
new file mode 100644
index 0000000..4ce11dc
--- /dev/null
+++ b/fiducia/examples/README.txt
@@ -0,0 +1,13 @@
+.. _examples-index:
+
+Examples
+========
+
+
+.. _general_examples:
+
+General examples
+----------------
+
+General-purpose and introductory examples from Fiducia
+ 
diff --git a/fiducia/examples/plot_example.py b/fiducia/examples/plot_example.py
new file mode 100644
index 0000000..c6b219e
--- /dev/null
+++ b/fiducia/examples/plot_example.py
@@ -0,0 +1,31 @@
+"""
+Test Example for Sphinx Docs
+===================================
+
+We create a test plot to see if sphinx gallery works
+"""
+
+
+import numpy as np
+import fiducia
+import matplotlib.pyplot as plt
+
+
+############################################################
+# Generating some data
+
+
+dataX = np.arange(10)
+dataY = dataX ** 2
+
+############################################################
+# Plotting the data
+
+
+plt.plot(dataX, dataY)
+plt.xlabel('X')
+plt.ylabel('Y')
+plt.title('Example data')
+plt.show()
+
+
diff --git a/fiducia/loader.py b/fiducia/loader.py
new file mode 100644
index 0000000..1153aab
--- /dev/null
+++ b/fiducia/loader.py
@@ -0,0 +1,476 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:20:37 2019
+
+Utilities for loading DANTE measurement and response function data.
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import scipy.sparse as sparse
+
+
+# custom modules
+from fiducia.misc import find_nearest
+import fiducia.pltDefaults
+
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["cleanupHeader",
+           "loadResponses",
+           "loadResponseUncertainty",
+           "readDanProcessed",
+           "signalsAtTime",
+           "signalInt",
+           "readDanteData",
+           ]
+
+
+def __responseName__(channelNum):
+    r"""
+    Convenience function for generating response function file name
+    given the DANTE channel number.
+    
+    Parameters
+    ----------
+    channelNum: int
+        DANTE channel number
+    
+    Returns
+    -------
+    fileName: str
+        The file name of the response file for channel 'channelNum'
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    baseName = 'DanteDataS_1_'
+    extension = '.dat'
+    fileName = baseName + str(channelNum) + extension
+    return fileName
+
+
+def __readResponse__(channelNum, directory):
+    r"""
+    Read a single DANTE channel response function file given the channel
+    number and path to the directory containing the response function files.
+    
+    Parameters
+    ----------
+    channelNum: int
+        DANTE channel number
+        
+    directory: str
+        Path to channel response function files
+        
+    Raises
+    ------
+    Exception
+        If file does not exist.
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    try:
+        fileName = __responseName__(channelNum)
+        responseArr = np.loadtxt(directory + fileName)
+        return responseArr
+    except IOError:
+        # file doesnt exist, return dummy data
+        Exception(f"File {fileName} does not exist!")
+
+
+def cleanupHeader(dataFrame):
+    r"""
+    Strip whitespace and rename DataFrame headers.
+
+    Parameters
+    ----------
+    dataFrame : pandas.core.frame.DataFrame
+        DataFrame to be cleaned.
+
+    Returns
+    -------
+    cleanedDataFrame : pandas.core.frame.DataFrame
+        DataFrame with stripped and renamed channel headers.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+
+    """
+    # stripping whitespace from header names
+    colNames = dataFrame.columns
+    colNamesStripped = [header.strip() for header in colNames]
+    renameStripDict = {colName: colNamesStripped[idx] for idx, colName in enumerate(colNames)}
+    df1 = dataFrame.rename(columns=renameStripDict)
+    # renaming channel headers from e.g. Ch2 to just 2
+    allChannels = np.arange(18) + 1
+    renameDict = {}
+    for channel in allChannels:
+        renameDict['Ch' + str(channel)] =  channel
+
+    cleanedDataFrame = df1.rename(columns=renameDict)
+    return cleanedDataFrame
+
+
+def loadResponses(channels, fileName, solid=True):
+    r"""
+    Load DANTE measurement data from files given the channels and path to the
+    directory containing the response function files. Returns a dataframe
+    with the data.
+    
+    Parameters
+    ----------
+    channels: list, numpy.ndarray
+        List or array of relevant channels
+        
+    fileName: str
+        Full path and filename of .csv file containing DANTE respones
+        functions.
+    
+    solid: Bool, optional
+        Includes solid angle in response function value if true. The default is true.
+        
+    Returns
+    -------
+    responseFrame:  pandas.core.frame.DataFrame
+        DataFrame with the response function data for the 'channels' requested
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+
+    solidAngles = fiducia.misc.solidAngles
+    # loading all the response functions
+    dataFrame = pd.read_csv(fileName)
+    #clean headers
+    cleanedFrame = cleanupHeader(dataFrame)
+    # filtering for channels we care about
+    # for this particular shot
+    colFilter = ['Energy(eV)'] + channels
+    responseFrame = cleanedFrame[colFilter].copy()
+    # convert energy column from strings to floats (if necessary)
+    if type(responseFrame['Energy(eV)'][0]) == str:
+        energyFloats = responseFrame['Energy(eV)'].str.replace(',', '').astype(float)
+        responseFrame.loc[:,'Energy(eV)'] = energyFloats
+    else:
+        energyFloats = responseFrame['Energy(eV)'].astype(float)
+        responseFrame.loc[:,'Energy(eV)'] = energyFloats
+
+    if solid:
+        for chan in channels:
+            #multiply each element by the corresponding channel solid angle
+            responseFrame.loc[:, chan] *= solidAngles[chan-1]
+            #save metadata that we already include solid angle
+            responseFrame.solid = True
+
+    return responseFrame
+
+
+def loadResponseUncertainty(responseFrame, fileName):
+    r"""
+    Load uncertainty percentages into a DataFrame.
+    
+    Parameters
+    ----------
+    responseFrame: pandas.core.frame.DataFrame
+        DataFrame to base the respones uncertainty frame on. 
+        
+    fileName: str
+        Full path and filename of .csv file containing DANTE response uncertainty
+        percentages functions.
+
+    Returns
+    -------
+    responseUncertaintyFrame : pandas.core.frame.DataFrame
+        DataFrame with each column being a channel and each element being the
+        channel's uncertainty percentage. Extended to match the photon energy
+        range in the response frame.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #read raw data
+    channelUncertaintyFrame = pd.read_csv(fileName)
+    #clean headers
+    cleanedFrame = cleanupHeader(channelUncertaintyFrame)
+    # filtering for columns we care about
+    channels = list(responseFrame.columns.values)
+    channels.remove('Energy(eV)')
+    responseUncertaintyFrame = responseFrame.copy()
+    #go through each column and fill each element with corresponding uncertainty
+    for chan in channels:
+        #multiply each responseFrame element by the percent uncertainty
+        responseUncertaintyFrame.loc[:, chan] *= cleanedFrame.loc[0, chan]/100
+    return responseUncertaintyFrame
+
+
+def readDanProcessed(channels, directory):
+    r"""
+    Loads DANTE measurement data from files given the channels and path to the
+    directory containing the reduced and aligned DANTE data. Returns a dataframe
+    with the data.
+    Note that this is *not* for raw data. It is for reading DANTE signals
+    that have already been processed by Dan Barnak's scripts.
+    
+    Parameters
+    ----------
+    channels: list, numpy.ndarray
+        List or array of relevant channels
+        
+    directory: str
+        Path to channel response function files
+    
+    Returns
+    -------
+    dataFrame : pandas.core.frame.DataFrame
+        Dataframe of aligned signals from Dan's analysis.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    # initialize finder for length of data array for longest channel
+    longestChLen = 0
+    # load all the arrays into a dict
+    responseDict = {}
+    for channel in channels:
+        responseData = __readResponse__(channel, directory)
+        responseDict[channel] = responseData
+        # getting length of array for longest channel
+        responseLen = len(responseData)
+        if responseLen > longestChLen:
+            longestChLen = responseLen
+    # initialize pandas dataframe
+    indices = np.arange(longestChLen)
+    dataFrame = pd.DataFrame(index=indices)
+    # add index axis label
+    dataFrame.index.name = 'indices'
+    # write channel response data into dataframe
+    for channel in channels:
+        chData = responseDict[channel]
+        # photon energy
+        dataFrame['Time' + str(channel)] = chData[:,0]
+        # channel response
+        dataFrame['Signal' + str(channel)] = chData[:,1]
+    return dataFrame
+
+
+def signalsAtTime(time,
+                  measurementFrame,
+                  channels,
+                  plot=False,
+                  method="interp"):
+    r"""
+    Get DANTE signals from each channel at a particular time. Default is
+    to return an interpolated value of the signal at the given time.
+    Alternatively, this function can return the nearest value in
+    the signal data array for the given time.
+    
+    
+    Parameters
+    ----------
+    time: float
+        Time for which we want DANTE signals (in ns).
+        
+    measurementFrame: pandas.core.frame.DataFrame
+        Pandas dataframe containing DANTE measurement data. See
+        readDanteData() and readDanProcessed().
+        
+    plot: Bool
+        When True, plots DANTE signals vs channel index at a particular time.
+        
+    method: str
+        Either 'nearest' or 'interp'. 'nearest' finds the nearest point in the
+        DANTE signal to the given time. 'interp' returns an interpolated
+        signal value for the given time. Default is 'interp'.
+    
+    Returns
+    -------
+    signals : numpy.ndarray
+        Dante signals for each channel at a particular time step.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    chLen = len(channels)
+    signals = np.zeros(chLen)
+    for idx, channel in enumerate(channels):
+        if method == "nearest":
+            timeIdx, _ = find_nearest(array=measurementFrame['Time' + str(channel)],
+                                      value=time)
+            signals[idx] = measurementFrame['Signal' + str(channel)][timeIdx]
+        elif method == "interp":
+            signals[idx] = np.interp(x=time,
+                                     xp=measurementFrame['Time' + str(channel)],
+                                     fp=measurementFrame['Signal' + str(channel)])
+        else:
+            raise Exception(f"Method {method} not found!")
+        
+    if plot:
+        plt.scatter(channels, signals)
+        plt.xticks(channels)
+        plt.xlabel('DANTE channel')
+        plt.ylabel('Signal (V)')
+        plt.title(f'DANTE signals @ t = {time} ns')
+        plt.show()
+    return signals
+
+
+def signalInt(channels, measurementFrame, tStart, tEnd):
+    r"""
+    Get time-integrated Dante signals for a specified time interval. Used in
+    getting time-integrated spectrum from the unfold.
+    
+    Parameters
+    ----------
+    measurementFrame: pandas.core.frame.DataFrame
+        Pandas dataframe containing DANTE measurement data. See
+        loadDanteData().
+    tStart: float
+        Lower bound for time integration.
+        
+    tEnd: float
+        Upper bound for time integration
+    
+    Returns
+    -------
+    signalInt: numpy.ndarray
+        Time integrated Dante signals for each channel.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    from scipy import integrate
+    
+    chLen = len(channels)
+    signalInt = np.zeros(chLen)
+    for idx, channel in enumerate(channels):
+        timeseries = measurementFrame['Time' + str(channel)]
+        chanseries = measurementFrame['Signal' + str(channel)]
+        timeIdx1, _ = find_nearest(array=timeseries, value=tStart)
+        timeIdx2, _ = find_nearest(array=timeseries, value=tEnd)
+        signalInt[idx] = integrate.simps(y=chanseries[timeIdx1:timeIdx2],
+                                         x=timeseries[timeIdx1:timeIdx2])
+    return signalInt
+
+
+def readDanteData(filePath):
+    r"""
+    Reads Dante .dat file and returns header info and channel signals
+    as two separate pandas dataframes.
+    
+    Parameters
+    ----------
+    filePath: str
+        Full path to the Dante .dat file.
+    
+    Returns
+    -------
+    headerFrame: pandas.core.frame.DataFrame
+        Header of Dante data file. This typically include information
+        about the various components used in each Dante channel, such
+        as oscilloscopes, XRDs, etc.
+    
+    dataFrame: pandas.core.frame.DataFrame
+        Dante data.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # reading the entire dante file
+    dataAndHeaderFrame = pd.read_csv(filePath, sep='\t', header=None)
+    # generating dante dataframe header names
+    headerNames1st = [str(num) for num in np.arange(18) + 1]
+    headerNames2nd = [str(num) + ' bkg' for num in np.arange(18) + 1]
+    headerName = headerNames1st + headerNames2nd
+    # replacing dataframe column header names for more intuitive access
+    dataAndHeaderFrame.columns = headerName
+    # splitting into header frame and measurement data frame
+    headerLen = 18
+    headerFrame = dataAndHeaderFrame[:][:headerLen]
+    dataFrame = dataAndHeaderFrame[:][headerLen:]
+    # replacing row names for header frame
+    indexNamesReplace = {0:'Signal Cable',
+                         1:'Attenuator 1',
+                         2:'Attenuator 2',
+                         3:'Attenuator 3',
+                         4:'Attenuator 4',
+                         5:'Jumper Cable',
+                         6:'XRD SN',
+                         7:'Mirror SN',
+                         8:'Filter 1 SN',
+                         9:'Filter 2 SN',
+                         10:'Filter 3 SN',
+                         11:'Fiducial T',
+                         12:'Scope type',
+                         13:'Full scale Hor time',
+                         14:'#Hor Pts',
+                         15:'Full Scale Vert mV',
+                         16:'HV bias for XRDs',
+                         17:'(unused field)'}
+    headerFrame.rename(index=indexNamesReplace, inplace=True)
+    return headerFrame, dataFrame
diff --git a/fiducia/main.py b/fiducia/main.py
new file mode 100644
index 0000000..27c3529
--- /dev/null
+++ b/fiducia/main.py
@@ -0,0 +1,781 @@
+
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 25 12:22:01 2019
+
+FIDUCIA: Filtered Diode Unfolder (using) Cubic Spline Algorithm
+
+DANTE spectrum deconvolver based on cubic splines method [1]. Translated
+from Dan Barnak's Mathematica code.
+
+DANTE channels are bounded by edge absorption feature (knot point) due to
+filter for the respective channel. Cubic splines representing the estimated
+spectrum are fitted in each spectral region bounded by knot points. The
+detector signal for each channel is then equal to the response function of
+the detector folded with the matrix representation of the cubic spline. A
+triadiagonal matrix representation of the cubic spline equation is used to
+make the problem numerically tractable. This way a matrix inversion can be
+used to solve for the unknown coefficients in the cubic spline equation,
+using the measured signals. These coefficients are then plugged back into
+the cubic spline equation over each interval (between knot points) to
+make a piecewise reconstruction of the x-ray spectrum at each time step.
+    
+
+References
+----------
+
+Cubic spline deconvolution method
+[1] J. P. Knauer and N. C. Gindele. Temporal and spectral deconvolution of
+data from diamond, photoconductive devices. Rev. Sci. Instrum. 75, 3714 (2004)
+https://doi.org/10.1063/1.1785274
+
+Error propagation for cubic spline deconvolution method
+[2] D. L. Fehl and F. Briggs. Verification of unfold error estimates in the
+unfold operator code. Rev. Sci. Instrum. 68, 890 (1997)
+https://doi.org/10.1063/1.1147713
+    
+Useful description of cubic spline matrix representation
+[3] http://mathworld.wolfram.com/CubicSpline.html
+
+Paper comparing cubic splines unfolds to other methods
+[4] D. H. Barnak, J. R. Davies, J. P. Knauer, and P. M. Kozlowski. Soft
+x-ray spectrum unfold of K-edge filtered x-ray diode arrays using
+cubic splines. Submitted to Review of Scientific Instruments in 2020.
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import matplotlib.pyplot as plt
+import astropy.units as u
+from astropy.constants import sigma_sb, k_B
+import xarray as xr
+
+# custom modules
+from fiducia.rawProcess import (loadCorrected,
+                                hysteresisCorrect,
+                                align,
+                                constructMeasurementFrame)
+from fiducia.loader import signalsAtTime, loadResponses
+from fiducia.cspline import knotSolve, reconstructSpectrum
+from fiducia.visualization import plotStreak, plotTraces, plotResponse
+from fiducia.response import knotFind
+from fiducia.error import trapzVariance
+import fiducia.pltDefaults
+
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["simulateSignal",
+           "inferRadTemp",
+           "inferPower",
+           "analyzeSpectrum",
+           "analyzeStreak",
+           "feelingLucky",
+           ]
+
+
+def simulateSignal():
+    """
+    Takes the inferred spectrum and folds it with the instrument
+    function to retrieve the forward propagated signal for each
+    DANTE channel.
+    
+    Parameters
+    ----------
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    return
+
+
+def inferRadTemp(power, area, angle, powerUncertainty=None):
+    r"""
+    Gets the inferred radiation temperature by calculating in from radiated
+    power through the Stefan-Boltzmann Law.
+    
+    Parameters
+    ----------
+    power: float, np.ndarray
+        Total radiated power as a function of time calculated from unfolded
+        spectra. See :func:`main.inferPower()`.
+
+    area: float
+        Area of emitting surface in units of mm^2. For hohlraums/halfraums,
+        this is the area of the LEH.
+
+    angle: float
+        Angle between the surface area normal and the Dante line of sight in 
+        degrees. Usually 37.4 degrees for hohlraums/halfraums. Must be
+        between 0 and 90 degrees.
+
+    powerUncertainty: float, np.ndarray, optional
+        Uncertainty in total radiated power as a function of time calculated from unfolded
+        spectra. See :func:`main.inferPower()`. The default is None.
+
+    Returns
+    -------
+    tRad : numpy.ndarray
+        Radiation temperature of the blackbody emitter.
+        
+    tRadVariance: numpy.ndarray
+        Variance :math:`\sigma^2` on the radiation temperature.
+
+    Notes
+    -----
+    Total x-ray flux (power) from a black body emitter is given by:
+        
+    .. math::
+        P = \sigma_{SB} A \cos(\theta) T^4
+        
+        
+    Where P = power, :math:`\sigma_{SB}` = Stefan-Boltzmann constant, A is the
+    area of radiating surface, :math:`\theta` is the viewing angle between the
+    surface area normal and the Dante line-of-sight, T is the radiation
+    temperature of the black body emitter.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # limit viewing angles to 90 degrees.
+    if not 0 <= angle <= 90:
+        raise Exception("Angle must be between 0 and 90 degrees.")
+    
+    # when power uncertainty is not given, it is set to zero.
+    if powerUncertainty is None:
+        if type(power) is np.ndarray:
+            powerUncertainty = np.zeros(power.shape)
+        else:
+            powerUncertainty = 0
+    
+    cosAngle = np.cos(np.deg2rad(angle))
+    #Stefan-Boltzmann constant
+    sb = sigma_sb.to(u.GW/(u.mm ** 2 * u.K ** 4))
+    sigma_sb_unitless = sb.value
+    kb = k_B.to(u.eV/u.K).value
+    denom = area * cosAngle * sigma_sb_unitless
+    # value to be quadratic rooted for getting radiation temperature
+    val2Root = power / denom
+    # quadratic rooting to get Trad in Kelvins, and converting from
+    # Kelvins to electronvolts.
+    tRad = np.power(val2Root, 0.25) * kb
+    tRadVariance = 0.0
+    if np.shape(powerUncertainty) == np.shape(power):
+#        tRadVariance = kb ** 2 / np.sqrt(denom) * (0.25 * power ** (-0.75) * powerUncertainty) ** 2
+        tRadVariance = (0.25 * tRad * powerUncertainty / power) ** 2
+        # cleaning up cases where we get nan because input power is zero
+        tRadVariance = np.nan_to_num(tRadVariance)
+    return tRad, tRadVariance
+
+
+def inferPower(energies, spectra, spectraUncertainty=None):
+    r"""
+    Gets the inferred total radiation power as a function of time.
+    
+    Parameters
+    ----------
+    energies: numpy.ndarray
+        Photon energies corresponding to input spectrum
+        
+    spectra: numpy.ndarray
+        Spectral Flux values as a function of photon energy in units of
+        (GW/sr/eV)
+    
+    Returns
+    -------
+    power: numpy.ndarray
+        Total x-ray power (flux) as a function of time.
+        
+    powerVariance: numpy.ndarray
+        Variance :math:`\sigma^2` on total x-ray power.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    if spectraUncertainty is None:
+        spectraUncertainty = np.zeros(spectra.shape)
+
+    timesLength = np.shape(spectra)[1]
+    power = np.zeros(timesLength)
+    powerVariance = np.zeros(timesLength)
+    for idx in range (0, timesLength):
+        power[idx] = np.trapz(y = spectra[:, idx], x=energies[:, idx])
+        powerVariance[idx] = trapzVariance(spectraUncertainty[:, idx],
+                                           x=energies[:, idx])
+    return power, powerVariance
+
+
+def analyzeSpectrum(channels,
+                    knots,
+                    detArr,
+                    detArrBoundaryCol,
+                    detArrVarianceBoundaryCol,
+                    detArrInv,
+                    stdDetArrInv,
+                    measurementFrame,
+                    time,
+                    signalsUncertainty = None,
+                    yGuess=0,
+                    boundary="y0",
+                    nPtsIntegral=100,
+                    nPtsSpectrum=100,
+                    plotSignal=False,
+                    plotKnots=False,
+                    plotSpectrum=True):
+    r"""
+    Given the response function file and the DANTE measurement data file,
+    run cubic spline analysis to reconstruct spectrum for a given time.
+    
+    
+    Parameters
+    ----------
+    channels: list, numpy.ndarray
+        List or array of relevant DANTE channel numbers.
+    
+    responseFrame: pandas.core.frame.DataFrame
+        Pandas dataFrame containing response functions for each DANTE
+        channel. See loadResponses().
+        
+    knots: list, numpy.ndarray
+        List or array of knot point photon energy value. See knotFind().
+    
+    detArr : xarray.DataArray
+        Matrix representing the spectrally integrated folding of the detector
+        response with a cubic spline interpolation of the x-ray spectrum.
+        2D array of channels and knot points of shape (n, n).
+ 
+    detArrBoundaryCol : xarray.DataArray
+        Column of cublic spline matrix corresponding to the knots at the boundary
+        chosen with `boundary`.
+ 
+    detArrVarianceBoundaryCol: xarray.DataArray
+        Column of variances in the cublic spline matrix corresponding to the 
+        knots at the boundary chosen with `boundary`.   
+ 
+    detArrInv : xarray.DataArray
+        Inversion of detArr, with the column corresponding to boundary removed so detArr is invertible.
+   
+    stdDetArrInv : xarray.DataArray
+        Array of the standard deviation of each element in detArrInv based on variance
+        using the `responseUncertaintyFrame` propagated with Monte Carlo. 
+    
+    measurementFrame: pandas.core.frame.DataFrame
+        Pandas dataframe containing DANTE measurement data. See
+        readDanteData() and readDanProcessed().
+        
+    time: float
+        Time for which we want DANTE signals (in ns).
+
+    signalsUncertainty: numpy.ndarray, optional
+        One dimensional array with each element corresponding to the uncertainty
+        each signal. The default is None.        
+        
+    yGuess: float, optional
+        Guess for position of boundary knot point. Default 0.
+        
+    boundary: str, optional
+        Choose whether yGuess corresponds to :math:`y_0` (lowest photon energy) or
+        :math:`y_{n+1}` (highest photon energy) boundary condition. This should
+        correspond to the photon energy value given in knots. Options are `y0`
+        or `yn+1`. Default 'y0'.
+    
+    nPtsIntegral: int, optional
+        Number of points used in computing the integral. Default is 100.
+        
+    nPtsSpectrum: int, optional
+        Number of points to use in reconstructing the spectrum. Default is 100.
+    
+    plotKnots: Bool, optional
+        Flag for plotting the Dante signal at the given time across all
+        channels. Default is False.
+    
+    plotKnots: Bool, optional
+        Flag for plotting just the solved knot points. Default is False.
+        
+    plotSpectrum: Bool, optional
+        Flag for plotting the unfolded spectrum. Default is True.
+    
+    Returns
+    -------    
+    knotsYAll: numpy.ndarray
+        Array of knot point intensity values with yGuess appended.
+        See knotSolve() and analyzeSpectrum().
+
+    knotsYVariance: numpy.ndarray
+        Array of knot point intensity uncertainty values with yGuess appended.
+        See knotSolve() and analyzeSpectrum().
+    
+    photonEnergies: numpy.ndarray
+        Photon energy axis of unfolded spectrum.
+        
+    intensities: numpy.ndarray
+        Spectral intensity axis of unfolded spectrum.
+        
+    intensitiesVariance: numpy.ndarray
+        Uncertaitny (1 :math:`\sigma`) on spectral intensity values.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # Select a particular time step for generating x-ray spectrum from DANTE
+    signals = signalsAtTime(time=time,
+                            measurementFrame=measurementFrame,
+                            channels=channels,
+                            plot=plotSignal,
+                            method="interp")
+    # solve knots and plot
+    # if signalsUncertainty is None:
+    #     signalsUncertainty = np.ones(signals.shape)
+        
+    knotsY, knotsYVariance = knotSolve(signals,
+                                       detArr, detArrBoundaryCol,
+                                       detArrVarianceBoundaryCol,
+                                       detArrInv, stdDetArrInv,
+                                       signalsUncertainty,
+                                       yGuess=yGuess,
+                                       npts=nPtsIntegral)
+    #Root the variance to get the uncertainty (sigma)
+    knotsYUncertainty = np.sqrt(knotsYVariance)
+    
+    # prepending y0 guess to form full array of y knot values equal in
+    # length to knot energy values
+    if boundary == "y0":
+        knotsYAll = np.append([yGuess], knotsY)
+        knotsYUncertaintyAll = np.append([yGuess], knotsYUncertainty)
+
+    elif boundary == "yn+1":
+        knotsYAll = np.append(knotsY, [yGuess])
+        knotsYUncertaintyAll = np.append(knotsYUncertainty, [yGuess])
+    else:
+        raise Exception(f"No method found for boundary {boundary}.")
+    
+    # making simple plot at this intermediate step to get rough idea
+    # of the shape of the spectrum
+    if plotKnots:
+            plt.plot(knots, knotsYAll)
+            plt.xlabel("Energy (eV)")
+            plt.ylabel("X-ray spectrum")
+            plt.title("cubic spline knots")
+            plt.show()
+    
+    # reconstruct spectrum using inferred knot points
+    chLen = len(channels)
+    photonEnergies, intensities, intensitiesVariance = reconstructSpectrum(chLen,
+                                                                           knots,
+                                                                           knotsYAll,
+                                                                           knotsYUncertaintyAll,
+                                                                           nPtsSpectrum,
+                                                                           plot=plotSpectrum)
+    photonEnergiesPlus, intensitiesPlus, _ = reconstructSpectrum(chLen,
+                                                                 knots,
+                                                                 knotsYAll+knotsYUncertaintyAll,
+                                                                 knotsYUncertaintyAll,
+                                                                 nPtsSpectrum,
+                                                                 plot=False)
+    photonEnergiesMinus, intensitiesMinus, _ = reconstructSpectrum(chLen,
+                                                                   knots,
+                                                                   knotsYAll-knotsYUncertaintyAll,
+                                                                   knotsYUncertaintyAll,
+                                                                   nPtsSpectrum,
+                                                                   plot=False)
+    #subtract to just get the delta from intensities. plot_line_shaded just wants the difference
+    intensitiesPosDiff = intensitiesPlus-intensities
+    intensitiesNegDiff = intensities - intensitiesMinus
+    if plotSpectrum:
+        fiducia.pltDefaults.plot_line_shaded(photonEnergies, 
+                                             intensities, 
+                                             intensitiesPosDiff, 
+                                             intensitiesNegDiff,)
+        plt.xlabel("Photon Energy (eV)")
+        plt.ylabel("Spectrum (GW/sr/eV)")
+        plt.title("Spectra with error bars")
+        plt.show()    
+    return knotsYAll, knotsYVariance, photonEnergies, intensities, intensitiesVariance
+
+
+def analyzeStreak(channels,
+                  responseFrame,
+                  knots,
+                  detArr,
+                  detArrBoundaryCol,
+                  detArrVarianceBoundaryCol,
+                  detArrInv, 
+                  stdDetArrInv,
+                  measurementFrame,
+                  timeStart,
+                  timeStop,
+                  timeStep,
+                  signalsUncertainty=None,
+                  yGuess=0,
+                  boundary="y0",
+                  nPtsIntegral=100,
+                  nPtsSpectrum=100):
+    r"""
+    Given the response function file and the DANTE measurement data file,
+    run cubic spline analysis to reconstruct spectrum for a given time.
+    
+    Parameters
+    ----------
+    channels: list, numpy.ndarray
+        List or array of relevant DANTE channel numbers.
+    
+    responseFrame: pandas.core.frame.DataFrame
+        Pandas dataFrame containing response functions for each DANTE
+        channel. See loadResponses().
+        
+    knots: list, numpy.ndarray
+        List or array of knot point photon energy value. See knotFind().
+    
+    detArr : xarray.DataArray
+        Matrix representing the spectrally integrated folding of the detector
+        response with a cubic spline interpolation of the x-ray spectrum.
+        2D array of channels and knot points of shape (n, n).
+ 
+    detArrBoundaryCol : xarray.DataArray
+        Column of cublic spline matrix corresponding to the knots at the boundary
+        chosen with `boundary`.
+ 
+    detArrVarianceBoundaryCol: xarray.DataArray
+        Column of variances in the cublic spline matrix corresponding to the 
+        knots at the boundary chosen with `boundary`.   
+        
+    detArrInv : xarray.DataArray
+        Inversion of detArr, with the column corresponding to boundary removed so detArr is invertible.
+   
+    stdDetArrInv : xarray.DataArray
+        Array of the standard deviation of each element in detArrInv based on variance
+        using the `responseUncertaintyFrame` propagated with Monte Carlo. 
+    
+    measurementFrame: pandas.core.frame.DataFrame
+        Pandas dataframe containing DANTE measurement data. See
+        :func:`loader.readDanteData` and :func:`readDanProcessed`.
+        
+    timeStart: float
+        Start time for producing temporally streaked DANTE spectra (in ns).
+        
+    timeStop: float
+        End time for producing temporally streaked DANTE spectra (in ns).
+        
+    timeStep: float
+        Time step size for producing temporally streaked DANTE spectra (in ns).
+        
+    signalsUncertainty: numpy.ndarray, optional
+        One dimensional array with each element corresponding to the uncertainty
+        each signal. The default is None. 
+    
+    yGuess: float, optional
+        Guess for position of boundary knot point. Default is 1e-77.
+        
+    boundary: str, optional
+        Choose whether yGuess corresponds to :math:`y_0` (lowest photon energy) or
+        :math:`y_{n+1}` (highest photon energy) boundary condition. This should
+        correspond to the photon energy value given in knots. Options are `y0`
+        or `yn+1`. Default 'y0'.
+    
+    nPtsIntegral: int, optional
+        Number of points used in computing the integral. Default is 100.
+        
+    nPtsSpectrum: int, optional
+        Number of points to use in reconstructing the spectrum. Default is 100.
+    
+    Returns
+    -------
+    times
+    
+    
+    energies
+    
+    
+    spectra
+    
+    
+    spectraVariance
+        
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+        
+    """
+    chLen = len(channels)
+    # generating evenly spaced timesteps
+    times = np.arange(timeStart, timeStop + timeStep, timeStep)
+    timesLen = len(times)
+    # initializing array for holding spectral values
+    spectralPts = chLen * nPtsSpectrum - chLen
+    energies = np.zeros((spectralPts, timesLen))
+    spectra = np.zeros((spectralPts, timesLen))
+    #TODO implement spectraVariance once DotVariance works on N-D arrays
+    spectraVariance = np.zeros((spectralPts, timesLen))
+    
+    for idt, time in enumerate(times):
+        # process the spectrum for the given time
+        results = analyzeSpectrum(channels,
+                                  knots,
+                                  detArr, 
+                                  detArrBoundaryCol,
+                                  detArrVarianceBoundaryCol,
+                                  detArrInv, 
+                                  stdDetArrInv,
+                                  measurementFrame,
+                                  time,
+                                  signalsUncertainty=signalsUncertainty,
+                                  yGuess=yGuess,
+                                  boundary=boundary,
+                                  nPtsIntegral=100,
+                                  nPtsSpectrum=100,
+                                  plotKnots=False,
+                                  plotSpectrum=False)
+        knotsYAll, knotsYVariance, photonEnergies, intensities, intensitiesVariance = results
+        # including unfolding spectrum for time step into array of spectra
+        # for all time steps
+        energies[:, idt] = photonEnergies
+        spectra[:, idt] = intensities
+        spectraVariance[:, idt] = intensitiesVariance
+        print(f"Completed time step {time} ns.")    
+    # plotting streaked spectrum
+    plotStreak(times, energies, spectra)
+    return times, energies, spectra, spectraVariance
+
+
+def feelingLucky(dataFile,
+                 attenuatorsFile,
+                 offsetsFile,
+                 responseFile,
+                 csplineDatasetFile,
+                 channels,
+                 area,
+                 angle,
+                 signalsUncertainty=None,
+                 peaksNum=2):
+    r"""
+    Attempt processing dante signals given dante data file and calibration
+    files using sensible defaults.
+    
+    
+    Parameters
+    ----------
+    dataFile: str
+        Full path to the Dante .dat file containing dante signals from LLE
+        site.
+        
+    attenuatorsFile: str
+        Full path to file containing attenuator serial numbers and
+        corresponding attenuation factors.
+        
+    offsetsFile: str
+        Full path to file containing oscilloscope channel offsets in
+        time and voltage.
+        
+    responseFile: str
+        Full path and filename of .csv file containing DANTE respones
+        functions corresponding to dataFile.
+    
+    csplineDatasetFile : str
+        File pointing to the path of the saved dataset containing ''detArr'',
+        ''detArrBoundaryCol'',  ''detArrInv'', and ''stdDetArrInv''. 
+        See :func:'error.analyzeSpectrumUncertainty()'.
+  
+    channels: list, numpy.ndarray
+        List or array of relevant channels for which to apply analysis.
+        
+    area: float
+        Area of emitting surface in units of mm^2. For hohlraums/halfraums,
+        this is the area of the LEH. Used in Trad calculation.
+        
+    angle: float
+        Angle between the surface area normal and the Dante line of sight in 
+        degrees. Usually 37.4 degrees for hohlraums/halfraums. Used in Trad 
+        calculation.
+    
+    signalsUncertainty: numpy.ndarray, optional
+        One dimensional array with each element corresponding to the uncertainty
+        each signal. The default is None.
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #load cspline matrix related files
+    csplineDataset =  xr.open_dataset(csplineDatasetFile).load()
+    detArr = csplineDataset['detArr']
+    detArrBoundaryCol = csplineDataset['detArrBoundaryCol']
+    detArrVarianceBoundaryCol = csplineDataset['detArrVarianceBoundaryCol']
+    detArrInv = csplineDataset['detArrInv']
+    stdDetArrInv = csplineDataset['stdDetArrInv']
+    #close out of the files
+    csplineDataset.close()
+    
+    #find what boundary was used when calculating cspline matrices
+    boundary = csplineDataset.attrs['boundary']
+    
+    # loading data and applying corrections.
+    timesFrame, dfAtten, onChList, hf, dfVolt = loadCorrected(danteFile=dataFile,
+                                                              attenuatorsFile=attenuatorsFile,
+                                                              offsetsFile=offsetsFile,
+                                                              plot=True)
+    
+    # removing hysteresis and background with a polynomial fit
+    dfPoly = hysteresisCorrect(timesFrame=timesFrame,
+                               df=dfAtten,
+                               channels=channels,
+                               order=5,
+                               prominence=0.2,
+                               width=10,
+                               avgMult=1)
+    
+    # aligning signals to peak
+    # aligning to 1e-9 seconds by default.
+    # looking for 2 peaks and aligning to the first (0th index) peak.
+    timesAligned = align(timesFrame=timesFrame,
+                         df=dfPoly,
+                         channels=channels,
+                         peaksNum=peaksNum,
+                         peakAlignIdx=0,
+                         referenceTime=1e-9,
+                         prominence=0.01,
+                         width=10,
+                         avgMult=1.5)
+    
+    # constructing dataframe for passing to analyzeStreak()
+    measurementFrame = constructMeasurementFrame(timesFrame=timesAligned,
+                                                 df=dfPoly,
+                                                 channels=channels)
+        
+    # testing plot traces
+    plotTraces(channels, measurementFrame, scale='log')
+    
+    # loading dante responses
+    responseFrame = loadResponses(channels, responseFile)
+
+    # getting knots from response function edges
+    # currently not forcing any knots
+    knots = knotFind(channels=channels,
+                      responseFrame=responseFrame,
+                      forceKnot=np.array([]),
+                      knotBoundary=0,
+                      boundary=boundary)
+    
+    # print(f"Channels: {channels}")
+    # print(f"Channels length: {len(channels)}")
+    # print(f"Knot energies (eV): {knots}")
+    # print(f"Knot energies length: {len(knots)}")
+    
+    # plotting all relevant response functions
+    # plotResponse(channels=channels, 
+    #              responseFrame=responseFrame,
+    #              knots=knots)
+    
+    # unfolding the spectrum for a particular time step
+    # currently unfolding at 1 ns
+    time = 1.0
+    spectrumResults = analyzeSpectrum(channels, 
+                                      knots, 
+                                      detArr,
+                                      detArrBoundaryCol,
+                                      detArrVarianceBoundaryCol,
+                                      detArrInv,
+                                      stdDetArrInv,
+                                      measurementFrame,
+                                      signalsUncertainty=signalsUncertainty,
+                                      time=time,
+                                      yGuess=0,
+                                      boundary=boundary,
+                                      plotKnots=True)
+    knotsYAll, knotsYVariance, photonEnergies, intensities, intensitiesVariance = spectrumResults
+    
+    # unfolding time resolved spectra
+    timeStart = -1
+    timeStop = 4
+    timeStep = 0.1
+    streakResults = analyzeStreak(channels,
+                                  responseFrame,
+                                  knots,
+                                  detArr, 
+                                  detArrBoundaryCol,
+                                  detArrVarianceBoundaryCol,
+                                  detArrInv, 
+                                  stdDetArrInv,
+                                  measurementFrame,
+                                  timeStart,
+                                  timeStop,
+                                  timeStep,
+                                  signalsUncertainty=signalsUncertainty,
+                                  yGuess=0,
+                                  boundary=boundary,
+                                  nPtsIntegral=100,
+                                  nPtsSpectrum=100)
+    times, energies, spectra, spectraVariance = streakResults
+
+
+    spectraUncertainty = np.sqrt(spectraVariance)
+    power, powerVariance = inferPower(energies, spectra, spectraUncertainty)
+    plt.scatter(times, power)
+    plt.xlabel("Time (ns)")
+    plt.ylabel("Radiated Power (GW/sr)")
+    plt.show()
+    powerUncertainty = np.sqrt(powerVariance)
+    fiducia.pltDefaults.plot_line_shaded(times, power, powerUncertainty)
+    plt.xlabel("Time (ns)")
+    plt.ylabel("Radiated Power (GW/sr)")
+    axes = plt.gca()
+    axes.set_ylim([0,None])
+    plt.show()
+    
+    tRad, tRadVariance = inferRadTemp(power, area, angle, powerUncertainty)
+    plt.scatter(times, tRad)
+    plt.xlabel("Time (ns)")
+    plt.ylabel("Radiation Temperature (eV)")
+    plt.show()
+    tRadUncertainty = np.sqrt(tRadVariance)
+    fiducia.pltDefaults.plot_line_shaded(times, tRad, tRadUncertainty)
+    plt.xlabel("Time (ns)")
+    plt.ylabel("Radiation Temperature (eV)")
+    axes = plt.gca()
+    axes.set_ylim([0,None])
+    plt.show()
+#    return None
+    return times, energies, spectra, power, tRad
diff --git a/fiducia/misc.py b/fiducia/misc.py
new file mode 100644
index 0000000..a890fcc
--- /dev/null
+++ b/fiducia/misc.py
@@ -0,0 +1,152 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:25:05 2019
+
+Miscellaneous utilities
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import astropy.units as u
+import copy
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["find_nearest",
+           "areDataFramesCompatible",
+           ]
+
+
+# solid angles of channels 1 thru 18 (from C Sorce)
+# This is in steradians (added by DHB 3/22/2019)
+solidAngles = [4.48e-6,
+               2.76e-6,
+               2.76e-6,
+               2.76e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6,
+               7.07e-6]
+
+# radius of OMEGA chamber
+chamberRadius = 180 * u.cm
+
+def find_nearest(array, value):
+    r"""
+    Find nearest value in array and return index, and value as a tuple.
+    
+    Parameters
+    ----------
+    array: list, numpy.ndarray
+        Array of values to be searched.
+        
+    value: int, float
+        Value for which this function will find the nearest value
+        in the array.
+    
+    Returns
+    -------
+    idx: int
+        Index at which nearest value to input value occurs in the array.
+    
+    array[idx]: int, float
+        The nearest value to the input value.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return idx, array[idx]
+
+
+#TODO Add support for channels of type ndarray
+def areDataFramesCompatible(channels, *frames):
+    r"""
+    Check DataFrame compatibility for specified channels.
+    
+    Checks if multiple pandas.core.frame.DataFrame objects are compatible and have
+    the channels that are requested. Checks that the DataFrames span the same
+    energy range. Returns true if the frames pass all checks, false otherwise.
+    
+    Parameters
+    ----------
+    channels: list
+        List of relevant channels
+    
+    *frames: pandas.core.frame.DataFrame
+        The DataFrames that you want to check for compatiblity with the relevant
+        channels
+        
+    Returns
+    -------
+    bool
+        True if frames are compatible with the requested channels, and False
+        otherwise.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #deep copy so we don't modify the passed channel list
+    ch = copy.deepcopy(channels)
+    ch.append('Energy(eV)')
+    cols = []
+    frameShapes = []
+    energyCols = []
+    for frame in frames:
+        cols.append(set(frame.columns.values))
+        frameShapes.append(frame.shape)
+        #check if the frame has an Energy(eV) col
+        if not 'Energy(eV)' in frame.columns:
+            #if not stop here, it definitely wont pass the next tests
+             return False
+        energyCols.append(frame['Energy(eV)'].values)
+     
+    #check if all frames have 2d shape
+    for shape in frameShapes:
+        if len(shape) != 2:
+            return False
+    
+    #check if the frames have the channels we need
+    overlap = set.intersection(*cols)
+    if not set(ch).issubset(overlap):
+        return False
+    
+    #check if the energy columns are equal
+    energyColToChangeAgainst = energyCols[0]
+    for energyCol in energyCols:
+        #check if the energy column lengths are equal
+        if len(energyColToChangeAgainst) != len(energyCol):
+            return False
+        #check if the energy column elements are equal
+        if not np.array_equal(energyColToChangeAgainst, energyCol):
+            return False
+    
+    #no more tests to run, DataFrames are compatible with the relevant channels
+    return True
diff --git a/fiducia/pltDefaults.py b/fiducia/pltDefaults.py
new file mode 100755
index 0000000..8563bdc
--- /dev/null
+++ b/fiducia/pltDefaults.py
@@ -0,0 +1,250 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Oct 27 02:37:12 2017
+
+Default plotting parameters
+
+@author: Pawel M. Kozlowski
+"""
+
+# importing Python modules
+import numpy as np
+import matplotlib 
+import matplotlib.pyplot as plt
+from itertools import cycle
+
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["plot_line_shaded",
+           "plot_scatter_bars"]
+
+
+# setting default plot properties
+plotFont = {'family' : 'serif',
+            'weight' : 'normal',
+            'size'   : 16}
+plt.rc('font', **plotFont)
+plt.rc('lines', linewidth=2)
+plt.rc('lines', markersize=8)
+#plt.rc('text', usetex=True)
+#plt.rc('text', usetex=False)
+# setting default figure size
+matplotlib.rcParams['figure.figsize'] = 10, 6
+
+validMathTextFonts = ['dejavusans',
+                      'dejavuserif',
+                      'cm',
+                      'stix',
+                      'stixsans',
+                      'custom']
+
+# setting latex rendered fonts to be same as regular fonts
+try:
+    matplotlib.rcParams['mathtext.fontset'] = 'dejavuserif'
+    matplotlib.rcParams['mathtext.rm'] = 'DejaVu Serif'
+    matplotlib.rcParams['mathtext.it'] = 'DejaVu Serif:italic'
+    matplotlib.rcParams['mathtext.bf'] = 'DejaVu Serif:bold'
+except:
+    print("Couldn't load dejavuserif fonts for plot defaults."
+          "Falling back to stix fonts.")
+    try:
+        matplotlib.rcParams['mathtext.fontset'] = 'stix'
+        matplotlib.rcParams['font.family'] = 'STIXGeneral'
+    except:
+        print("Couldn't load stix fonts for plot defaults.")
+
+# for testing whether matplotlib and python fonts match
+#plt.title(r'cm$\rm cm^{-3}$')
+
+
+#%% convenience functions for plotting
+
+# shaded error bars for line plot
+def plot_line_shaded(xData, yData, yErrsPos, yErrsNeg=[], label="", **kwargs):
+    """
+    Generate a line plot with shaded region representing y-error bars.
+    Can be run multiple times before plt.show(), to plot multiple data
+    sets on the same axes.
+    
+    Parameters
+    ----------
+    xData: numpy.ndarray
+        X-axis data to be plotted.
+        
+    yData: numpy.ndarray
+        Y-axis data to be plotted.
+        
+    yErrsPos: numpy.ndarray
+        Errors on yData.
+        
+    yErrsPos: numpy.ndarray
+        When errors on yData are asymmetric, these are the positive
+        side errors.
+        
+    yErrsNeg: numpy.ndarray
+        When errors on yData are asymmetric, these are the negative
+        side errors.
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #make errors symmetric unless specified otherwise
+    if yErrsNeg == []:
+        yErrsNeg = yErrsPos
+    # check that arrays are 1D
+    
+    # check that arrays are equal length
+    if not len(xData) == len(yData) == len(yErrsPos) ==len(yErrsNeg):
+        raise ValueError("Arrays must of equal length!")
+
+    # line plot
+    ax = kwargs.pop('ax', plt.gca())
+    base_line, = ax.plot(xData, yData, label=label, **kwargs)
+    # shaded error region
+    plt.fill_between(xData,
+                     yData - yErrsNeg,
+                     yData + yErrsPos,
+                     alpha=0.5,
+                     edgecolor=base_line.get_color(),
+                     facecolor=base_line.get_color())
+    return
+
+
+def plot_scatter_bars(xData, yData, yErrsPos, yErrsNeg=[], label="", **kwargs):
+    """
+    Generate a scatter plot with y-error bars.
+    Can be run multiple times before plt.show(), to plot multiple data
+    sets on the same axes.
+    
+    Parameters
+    ----------
+    xData: numpy.ndarray
+        X-axis data to be plotted.
+        
+    yData: numpy.ndarray
+        Y-axis data to be plotted.
+        
+    yErrsPos: numpy.ndarray
+        Errors on yData.
+        
+    yErrsPos: numpy.ndarray
+        When errors on yData are asymmetric, these are the positive
+        side errors.
+        
+    yErrsNeg: numpy.ndarray
+        When errors on yData are asymmetric, these are the negative
+        side errors.
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #make errors symmetric unless specified otherwise
+    if yErrsNeg == []:
+        yErrsNeg = yErrsPos
+    # check that arrays are 1D
+    
+    # check that arrays are equal length
+    if not len(xData) == len(yData) == len(yErrsPos) == len(yErrsNeg):
+        raise ValueError("Arrays must of equal length!")
+        
+    # line plot
+    ax = kwargs.pop('ax', plt.gca())
+    ax.errorbar(xData,
+                yData,
+                yerr=[yErrsNeg, yErrsPos],
+                label=label,
+                fmt='o',
+                fillstyle='none',
+                capsize=4,
+                elinewidth=1,
+                **kwargs)
+    return
+
+#%% testing custom plotting tools
+
+#xData = np.arange(10)
+#yData = xData ** 2
+#yData2 = xData ** 3
+#yErrs = np.ones_like(xData) * 30
+#plot_line_shaded(xData, yData, yErrs, label="#1", color='r')
+#plot_line_shaded(xData, yData2, yErrs, label="#2")
+#plt.ylabel("y")
+#plt.xlabel("x")
+#plt.legend(loc='upper left', frameon=False, labelspacing=0.001,
+#           fontsize=14, borderaxespad=0.4)
+#plt.show()
+#
+#
+#xData = np.arange(10)
+#yData = xData ** 2
+#yData2 = xData ** 3
+#yErrs = np.ones_like(xData) * 30
+#plot_scatter_bars(xData, yData, yErrs, label="#1", color='r')
+#plot_scatter_bars(xData, yData2, yErrs, label="#2")
+#plt.ylabel("y")
+#plt.xlabel("x")
+#plt.legend(loc='upper left', frameon=False, labelspacing=0.001,
+#           fontsize=14, borderaxespad=0.4)
+#plt.show()
+
+#%% style cycling
+# printing default colors
+default_colors = matplotlib.colors.cnames.keys()
+#print(f"Colors: {default_colors}")
+# printing default linestyles
+default_lines = matplotlib.lines.lineStyles.keys()
+#print(f"Lines: {default_lines}")
+# printing default marker styles
+default_markers = matplotlib.markers.MarkerStyle.markers.keys()
+#print(f"Markers: {default_markers}")
+
+# custom list of linestyles (excluding blank line styles)
+lines = ["-","--","-.",":"]
+linecycler = cycle(lines)
+#plt.figure()
+#for i in range(10):
+#    x = range(i,i+10)
+#    plt.plot(range(10),x,next(linecycler))
+#plt.show()
+    
+
+#%% Custom colors based on https://towardsdatascience.com/making-matplotlib-beautiful-by-default-d0d41e3534fd
+#CB91_Blue = '#2CBDFE'
+#CB91_Green = '#47DBCD'
+#CB91_Pink = '#F3A0F2'
+#CB91_Purple = '#9D2EC5'
+#CB91_Violet = '#661D98'
+#CB91_Amber = '#F5B14C'
+#
+#color_list = [CB91_Blue,
+#              CB91_Pink,
+#              CB91_Green,
+#              CB91_Amber,
+#              CB91_Purple,
+#              CB91_Violet]
+#plt.rcParams['axes.prop_cycle'] = plt.cycler(color="seaborn-colorblind")
+
+
+#%% built-in custom color cycles
+#print(plt.style.available)
+#plt.style.use('seaborn')
\ No newline at end of file
diff --git a/fiducia/rawProcess.py b/fiducia/rawProcess.py
new file mode 100644
index 0000000..6cbf320
--- /dev/null
+++ b/fiducia/rawProcess.py
@@ -0,0 +1,1700 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 13 16:43:39 2019
+
+Utilities for processing raw DANTE data. Typical steps include:
+    - attenuator correction
+    - background shot subtraction
+    - channel alignment (via e.g. peak finding)
+    - temporal axis calibration
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import find_peaks
+#from skimage import feature
+
+# custom modules
+from fiducia.misc import find_nearest
+from fiducia.loader import readDanteData
+import fiducia.pltDefaults
+
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["noScope",
+           "noXRD",
+           "onChannels",
+           "timesScope",
+           "voltageScale",
+           "bkgCorrect",
+           "offsetCorrect",
+           "attenuationFactors",
+           "attenuationCorrect",
+           "timeAvgBkg",
+           "avgBkgCorrect",
+           "polyBkg",
+           "signalEdges",
+           "polyBkgFrame",
+           "highestPeak",
+           "highestN",
+           "getPeaks",
+           "alignPeaks",
+           "constructMeasurementFrame",
+           "loadCorrected",
+           "hysteresisCorrect",
+           "align",
+           ]
+
+
+def noScope(hf):
+    r"""
+    Given a header frame, return a list of channels with no scope. These
+    are the channels that are off.
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Header frame from DANTE measurement data. See readDanteData().
+    
+    Returns
+    -------
+    set
+        Set of channels corresponding to oscilloscopes marked as "off"
+        in the Dante data file header.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    totalChannels = int(np.shape(hf)[1] / 2)
+    offScope = []
+    for idx in range(totalChannels):
+        ch = idx + 1
+        chStr = str(ch)
+        if hf[chStr]['Scope type'] == 9999:
+            offScope.append(ch)
+    return set(offScope)
+
+
+def noXRD(hf):
+    r"""
+    Given a header frame, return a list of channels with no XRD. If there is
+    a scope, then these channels may still register a signal!
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Header frame from DANTE measurement data. See readDanteData().
+    
+    Returns
+    -------
+    set
+        Set of channels corresponding to no XRDs marked in the Dante data file
+        header.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    totalChannels = int(np.shape(hf)[1] / 2)
+    offXRD = []
+    for idx in range(totalChannels):
+        ch = idx + 1
+        chStr = str(ch)
+        if hf[chStr]['XRD SN'] == 0:
+            offXRD.append(ch)
+    return set(offXRD)
+
+
+def onChannels(hf):
+    r"""
+    Given a header frame, return a list of which dante channels
+    were on for the shot.
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Header frame from DANTE measurement data. See readDanteData().
+    
+    Returns
+    -------
+    set
+        Set of channels corresponding oscilloscopes and XRDs marked as on
+        within the Dante data file header.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    offXRDSet = noXRD(hf)
+    offScopeSet = noScope(hf)
+    allOff = offXRDSet.union(offScopeSet)
+    allChSet = set(np.arange(18) + 1)
+    onChSet = allChSet - allOff
+    return onChSet
+
+
+def chStatus(hf):
+    r"""
+    Generates a list of statuses for each channel (on, no scope, or no XRD).
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Header frame from DANTE measurement data. See readDanteData().
+    
+    Returns
+    -------
+    status: dict
+        Statuses of Dante channels according to data file header.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    offXRDSet = noXRD(hf)
+    offScopeSet = noScope(hf)
+    allOff = offXRDSet.union(offScopeSet)
+    allChArr = np.arange(18) + 1
+    allChSet = set(allChArr)
+    onChSet = allChSet - allOff
+    status = {}
+    for chNum in allChArr:
+        if chNum in onChSet:
+            status[chNum] = "On"
+        elif chNum in offXRDSet:
+            status[chNum] = "No XRD"
+        elif chNum in offScopeSet:
+            status[chNum] = "No Scope"
+        else:
+            status[chNum] = "Off?"
+    return status
+
+
+def timesScope(hf):
+    r"""
+    Given a headerFrame, returns a timesFrame containing an array of
+    oscilloscope times for each channel and background shot in the
+    headerFrame.
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Header frame from DANTE measurement data. See readDanteData().
+    
+    Returns
+    -------
+    timesFrame: pandas.core.frame.DataFrame
+        Returns a timesFrame containing the corresponding times for each
+        oscilloscope trace contained in the header frame.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    timesFrame = pd.DataFrame(index=np.arange(1024), columns=hf.keys())
+    for key in hf.keys():
+        timeCode = str(hf[key]['Full scale Hor time'])
+        # if the time code is 99, it means the scope was off
+        if timeCode == '99':
+            print(f"No timespan for channel {key}.")
+        # extracting total time spanned by the oscilloscope trace based on the
+        # time code. For example, a time code of 208 means 2.0e-8 seconds.
+        else:
+            # assuming decimal point comes after first number, and all but the last
+            # number are part of the significand.
+            timeSpan = float(timeCode[0] + '.' + timeCode[1:-1] + 'e-' + timeCode[-1])
+            # getting total number of points recorded by oscilloscope
+            scopePts = hf[key]['#Hor Pts']
+            # getting length of time steps between points by dividing the total time
+            # over which oscilloscope trace is recorded by the total number of
+            # time steps. There is one fewer step than number of points.
+            timeStep = timeSpan / (scopePts - 1)
+            # constructing array of times for convience in analysis/plotting
+            timesFrame[key] = np.arange(scopePts) * timeStep
+
+    return timesFrame
+
+
+def voltageScale(hf, df):
+    r"""
+    Scales voltage (vertical) axis of dante signals based on information
+    contained in the header. Returns a dataframe with the dante signals
+    in units of volts. Also returns an errors/uncertainties frame
+    in units of volts, where the uncertainty due to the 11-bit ADC
+    converter has been calculated.
+    
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Header dataframe from dante .dat file. See readDanteData().
+        
+    df: pandas.core.frame.DataFrame
+        Dante dataframe. See readDanteData().
+    
+    Returns
+    -------
+    dfScaled: pandas.core.frame.DataFrame
+        Dante dataframe with signals in units of volts.
+        
+    errFrame: pandas.core.frame.DataFrame
+        Correspoding errors for dfScaled. Also in units of volts.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # number of bits in ADC converter
+    bits = 11
+    # number of bins/segments that voltage range is divided into. This
+    # defines the uncertainty in the measurements due to the oscilloscope
+    # and is one less than the total number of points spanning the range
+    # of the oscilloscope's ADC.
+    bitRange = 2 ** bits - 1
+    # initializing scaled dataframe and error frame with same structure
+    # as input dataframe.
+    dfScaled = pd.DataFrame().reindex_like(df)
+    errFrame = pd.DataFrame().reindex_like(df)
+    for key in hf.keys():
+        # channel voltage scale in units of volts
+        maxVoltage = hf[key]['Full Scale Vert mV'] / 1e3
+        # scale conversion factor
+        countsToVolts = maxVoltage / bitRange
+        # converting oscilloscope traces from counts to volts
+        dfScaled[key] = df[key] * countsToVolts
+        # uncertainties are constant and based on the ADC resolution
+        # for each channel
+        errFrame[key] = np.ones_like(df[key]) * countsToVolts
+    return dfScaled, errFrame
+
+
+def bkgCorrect(df, timesFrame):
+    r"""
+    Give a Dante data frame containing measurement data and background
+    shot data, remove the background from the data and return the corrected
+    data as a dataframe. Note that the returned dataframe is different
+    in a few ways from the input dataframe. First, the returned dataframe
+    is assumed to have strings as column headers, whereas the returned
+    dataframe will have integers (corresponding to dante channel
+    number) as the column headers. In addition, the input dataframe
+    will start indexing at some number above 0 (usually 18, due to the
+    header length), whereas the returned dataframe is re-indexed to
+    begin at 0.
+    
+    A dataframe with corresponding time scales to df is also passed to this
+    function for reindexing from strings to integers. This also acts as a
+    placeholder in case it is necessary to interpolate values if the
+    background shot and measurement shot timescales are not the same. Though
+    this type of interpolation is not currently implemented.
+    
+    
+    Parameters
+    ----------
+    df: pandas.core.frame.DataFrame
+        Dataframe of raw dante data. This should contain both the shot
+        measurement and the shot background as columns. See readDanteData().
+        The columns in this dataframe are assumed to be strings.
+        
+    timesFrame: pandas.core.frame.DataFrame
+        Dataframe containing time axis corresponding to dante signals in
+        df dataframe. See timesScope().
+    
+    Returns
+    -------
+    timesBkg: pandas.core.frame.DataFrame
+        Returns a dataframe of times corresponding to dfCorrected signals.
+        The columns in this dataframe are integers corresponding to Dante
+        channel number.
+    
+    dfCorrected: pandas.core.frame.DataFrame
+        Returns a dataframe of background subtracted dante signals. The
+        columns in this dataframe are integers corresponding to Dante
+        channel number.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # reassign index (dataframe currently starts at 18 or so due to header)
+    dfnew = df.set_index(np.arange(np.shape(df)[0]))
+    # get total number of channels in dataframe
+    totalChannels = int(np.shape(dfnew)[1] / 2)
+    # produce a list of all channels
+    channels = np.arange(totalChannels) + 1
+    # get the total number of datapoints in the dataframe
+    dataRange = np.arange(np.shape(dfnew)[0])
+    # initialize an empty dataframe for storing background corrected values.
+    dfCorrected = pd.DataFrame(index=dataRange, columns=channels)
+    # initializing a new timesframe
+    timesBkg = pd.DataFrame(index=dataRange, columns=channels)
+    for ch in channels:
+        # subtract shot background from shot measurement.
+        dfCorrected[ch] = dfnew[str(ch)] - dfnew[str(ch) + ' bkg']
+        # passing through time values to new frame
+        timesBkg[ch] = timesFrame[str(ch)]
+    return timesBkg, dfCorrected
+
+
+def offsetCorrect(df, timesFrame, offsetsFile):
+    r"""
+    Reads given offset correction file (.xls) and applies offsets to
+    dante measurement data given in dataframe. The input dataframe
+    should already be background corrected and scaled to units of volts, 
+    see bkgCorrect() and voltageScale().
+    Note that although timing offsets are also applied, they are not as
+    relevant since timing should be realigned to a fiducial peak
+    anyway.
+    Additional attenuation is not implemented and an error will be raised if
+    the offsets file contains attenuation values other than 1.
+    Returns a dataframe with applied offsets.
+    
+    Parameters
+    ----------
+    df: pandas.core.frame.DataFrame
+        Dante dataframe with background corrected values and scaled
+        to units of volts. See readDanteData(), bkgCorrect() and
+        voltageScale().
+        
+    timesFrame: pandas.core.frame.DataFrame
+        Dataframe containing time axis corresponding to dante signals in
+        df dataframe. See timesScope() and bkgCorrect().
+    
+    
+    Returns
+    -------
+    offsetsFile: str
+        Full path to .xls file containing dante channel offsets.
+    
+    dfOffset: pandas.core.frame.DataFrame
+        Dante dataframe with applied offset corrections
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # read excel file with offsets
+    offsetsFrame = pd.read_excel(offsetsFile)
+    offsets = offsetsFrame.set_index(np.arange(np.shape(offsetsFrame)[0]) + 1)
+    # initialize dataframe for storing offset corrected values
+    dfOffset = pd.DataFrame().reindex_like(df)
+    timeOffset = pd.DataFrame().reindex_like(timesFrame)
+    for key in df.keys():
+        # correct for channel offset
+        # Offsets are added to signal
+        dfOffset[key] = df[key] + offsets['OffsetV (V)'][key]
+        if offsets['additional attenuation'][key] != 1:
+            raise NotImplementedError("Non-zero attenuation in offsets "
+                                      "file not implemented!")
+        # applying time offsets
+        # need to scale by 1e-9 since timesFrame is given in seconds, whereas
+        # offsets are given in nanoseconds.
+        timeOffset[key] = timesFrame[key] + 1e-9 * offsets['OffsetT (ns)'][key]
+    return timeOffset, dfOffset
+
+
+def __getAttenuation__(attenuatorsFrame, serialNum):
+    r"""
+    Given dataframe with attenuator info and a serial number,
+    returns the attenuation factor
+    
+    Parameters
+    ----------
+    attenuatorsFrame: 
+        
+        
+    serialNum: 
+    
+        
+    Returns
+    -------
+    factorVal: 
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    factor = attenuatorsFrame[attenuatorsFrame['#'] == serialNum]['Factor']
+    factorVal = factor.values
+    return factorVal
+
+
+def attenuationFactors(hf, channels, attenuatorsPath):
+    r"""
+    Given a header frame, return the attenuation factors applied to each
+    channel.
+    
+    Parameters
+    ----------
+    hf: pandas.core.frame.DataFrame
+        Pandas dataframe containing dante header information. See
+        readDanteData().
+    
+    channels: set
+        Set of channels to be analyzed.
+        
+    attenuatorsPath: str
+        Full path to excel file containing attenuator serial numbers and
+        corresponding attenuation factors.
+    
+    Returns
+    -------
+    chFactors
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # load attenuators dataframe
+    attenuatorsFrame = pd.read_excel(attenuatorsPath)
+    # list of pandas dataframe indices containing attenuator serial numbers.
+    attenuatorIdxs = ['Attenuator 1',
+                      'Attenuator 2',
+                      'Attenuator 3',
+                      'Attenuator 4']
+    # initializing array to store attenuation factors corresponding to each
+    # channel.
+    chFactors = pd.DataFrame(1, index=['Attenuation Factor'], columns=channels)
+    # cycling over channels
+    for ch in channels:
+        # getting attenuator serial numbers
+        serials = hf[str(ch)][attenuatorIdxs]
+        # filtering null values (no attenuator)
+        realSerials = serials[serials != 0]
+        # fetching attenuation factors
+        factors = [__getAttenuation__(attenuatorsFrame, sn) for sn in realSerials]
+        # product sum of attenuation factors and saving to pandas
+        # dataframe. If the list is empty then prod will return a value
+        # of 1, which is exactly what we want!
+        chFactors[ch] = np.prod(factors)
+    return chFactors
+
+
+def attenuationCorrect(attenuatorsFile, hf, df, channels):
+    r"""
+    Given a Dante data frame and header frame, return a data frame with
+    attenuation corrections applied to each channel.
+    
+    Parameters
+    ----------
+    attenuatorsFile: str
+        Full path to .xls file containing attenuator serial numbers
+        and corresponding attenuation factors. See attenuationFactors().
+    
+    hf: pandas.core.frame.DataFrame
+        Header dataframe from dante .dat file. See readDanteData().
+        
+    df: pandas.core.frame.DataFrame
+        Dante dataframe. This frame should already be voltage scaled,
+        background corrected, and offset corrected. See readDanteData().
+        
+    channels: list
+        List of dante channels in df to be analyzed.
+    
+    Returns
+    -------
+    dfAtten: pandas.core.frame.DataFrame
+        Returns dataframe with attenuation corrected signal values for
+        the given channels.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    # getting attenuation factors
+    attFactors = attenuationFactors(hf, channels, attenuatorsFile)
+    # initializing dataframe for attenuation corrected values
+    dfAtten = pd.DataFrame(index=np.arange(len(df)),
+                           columns=channels)
+    for ch in channels:
+        # applying attenuation factors
+        dfAtten[ch] = df[ch] * attFactors[ch]['Attenuation Factor']
+    return dfAtten
+
+
+def timeAvgBkg(times, signals, timeStart, timeEnd):
+    r"""
+    Calculates time averaged background for given data.
+    
+    Parameters
+    ----------
+    times:
+        
+        
+    signals:
+        
+        
+    timeStart:
+        
+        
+    timeEnd:
+        
+    
+    Returns
+    -------
+    avg: 
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # getting nearest indices to start and stop times
+    idxStart, valueStart = find_nearest(times, timeStart)
+    idxEnd, valueEnd = find_nearest(times, timeEnd)
+    # averaging
+    avg = np.mean(signals[idxStart:idxEnd])
+    return avg
+
+
+def avgBkgCorrect(timesFrame, df, channels, timeLength=1e-9):
+    r"""
+    Applies background correction to bring the signal down to zero, based
+    on averaging the signal background over a section of time from earliest
+    time contained in timesFrame to earliest time plus timeLength.
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        Dataframe of time axis values corresponding to signals in df. These
+        should be in units of seconds.
+        
+    df: pandas.core.frame.DataFrame
+        Dataframe of dante signals. These should already be attenuation
+        corrected and in units of volts.
+        
+    channels: set
+        Set of channels to be analyzed.
+    
+    timeLength: float
+        Duration of time from initial time over which to take the average. In
+        units of seconds.
+        
+    dfAvg: pandas.core.frame.DataFrame
+        Returns a dataframe containing average background corrected
+        signals.
+    
+    Returns
+    -------
+    dfAvg:
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    dfAvg = pd.DataFrame().reindex_like(df)
+    for ch in channels:
+        # getting earliest time in frame.
+        timeStart = timesFrame[ch][0]
+        # added duration to earliest time to get end time for average.
+        timeEnd = timeStart + timeLength
+        # calculating average background by grabbing a section of the signal
+        avgBkg = timeAvgBkg(times=timesFrame[ch].values,
+                            signals=df[ch].values,
+                            timeStart=timeStart,
+                            timeEnd=timeEnd)
+        # applying correction
+        dfAvg[ch] = df[ch] - avgBkg
+    return dfAvg
+
+
+def polyBkg(time,
+            signal,
+            lowerEdge,
+            upperEdge,
+            order=3,
+            lowerLength=None,
+            upperLength=None,
+            plot=False):
+    r"""
+    Fit polynomial function to ends of the signal as an estimate of the
+    background signal + hyesteresis. Default is cubic fit.
+    
+    Parameters
+    ----------
+    time: numpy.ndarray
+        array of times corresponding to signal
+        
+    signal: numpy.ndarray
+        array of signal values for a single dante channel
+        
+    lowerEdge: int
+         Index of time array corresponding to lower edge of detected
+         signal. See signalEdges().
+        
+    upperEdge: int
+        Index of time array corresponding to upper edge of detected
+        signal. See signalEdges().
+        
+    order: int
+        Order of polynomial to be fitted to estimated background/hysteresis.
+        
+    lowerLength: int
+        Length over which to take the polynomial background fit on the
+        lower end (earlier in time) segment of the signal, with respect
+        to lowerEdge. Default is None, which just takes the first point
+        in the signal.
+        
+    upperLength: int
+        Length over which to take the polynomial background fit on the
+        upper end (later in time) segment of the signal, with respect
+        to upperEdge. Defualt is None, which then just picks the second
+        to last point in the signal.
+        
+    plot: bool
+        Flag for plotting polynomial fitted background signal. Default
+        is False.
+        
+    Returns
+    -------
+    time:
+        
+        
+    fitSignal: 
+        
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    if lowerLength == None:
+        lowerMin = 0
+    else:
+        lowerMin = lowerEdge - lowerLength
+    lowerMax = lowerEdge
+    upperMin = upperEdge
+    if upperLength == None:
+        upperMax = -1
+    else:
+        upperMax = upperEdge + upperLength
+    dataEdgeSignal = np.concatenate((signal[lowerMin:lowerMax], signal[upperMin:upperMax]))
+    dataEdgeTime = np.concatenate((time[lowerMin:lowerMax], time[upperMin:upperMax]))
+    # fit polynomial
+    fit = np.polyfit(dataEdgeTime, dataEdgeSignal, order)
+    polyObj = np.poly1d(fit)
+    fitSignal = polyObj(time)
+    
+    # plotting
+    if plot:
+        plt.plot(time, signal, label='signal')
+        plt.plot(time, fitSignal, '--', label='fit')
+        plt.scatter(dataEdgeTime, dataEdgeSignal, label='bkg')
+        plt.legend()
+        plt.show()
+    
+    return time, fitSignal
+
+
+
+def signalEdges(timesFrame,
+                df,
+                channels,
+                sigmaMult=3,
+                plot=False,
+                prominence=0.1,
+                width=10,
+                avgMult=1):
+    r"""
+    Determines locations and widths of peaks above the mean of the signal
+    for each dante channel. Edges of the signal containing region are then
+    obtained by moving sigmaMult peak widths away from the earliest and
+    latest peaks. Returns these lower and upper bound edges of the signal
+    containing region as a dataframe. These edges are useful for fitting
+    and removing the background/hysteresis.
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        A dataframe containing time axis values corresponding to signals in
+        df.
+    
+    df: pandas.core.frame.DataFrame
+        A dataframe of corrected/calibrated dante signal measurements.
+    
+    channels: list
+        List of channels in df for which edges will be determined.
+    
+    sigmaMult: float
+        Multiplier factor by which the lower and upper bounds of the signal
+        containing region are determined. The lower bound is determined
+        by sigmaMult times the width of the earliest peak away from the
+        earliest peak. The upper bound is determined by sigmaMult times the
+        width of the latest peak away from the latest peak. Default is 3
+        for approximately 3*sigma away from each peak.
+        
+    plot: bool
+        Flag for plotting peak locations and widths. Default is False.
+        
+    edgesFrame: pandas.core.frame.DataFrame
+        Lower and upper bound edges of the signal containing region for
+        each dante channel. The lower bound is in 0 index and the upper
+        bound is in 1 index. The bounds are given in index coordinates
+        and they have been rounded to the nearest point.
+        
+    prominence: float
+        Prominence threshold for identifying peaks in scipy's find_peaks().
+        
+    width: int
+        Width in index units for identifying peaks in scipy's find_peaks().
+        
+    avgMult: float
+        Multiplicative factor for setting minimum intensity threshold
+        for indentifying peaks in scipy's find_peaks(). This is a multiple
+        of the signal average.
+    
+    Returns
+    -------
+    edgesFrame: 
+        
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    edgesFrame = pd.DataFrame(index=[0,1], columns=df.keys())
+    for ch in channels:
+        times = timesFrame[ch]
+        signal = df[ch]
+        # finding peaks which are above the mean of the signal
+        peaks, properties = find_peaks(signal,
+                                       height=avgMult * np.mean(signal),
+                                       prominence=prominence,
+                                       width=width)
+        # determining lower and upper bounds of signal region by taking first
+        # and last peaks and shifting by distance of peak width
+        width1 = properties["widths"][0]
+        # width / 2 is kind of like sigma, and then that multiplied by 3 is
+        # like 3 sigma away.
+        lowerBnd = peaks[0] - sigmaMult * width1 / 2
+        width2 = properties["widths"][-1]
+        upperBnd = peaks[-1] + sigmaMult * width2 / 2
+        # saving these edges to the dataframe which is to be returned
+        edgesFrame[ch][0] = int(round(lowerBnd))
+        edgesFrame[ch][1] = int(round(upperBnd))
+        if plot:
+            # plotting
+            plt.plot(times, signal)
+            plt.plot(times[peaks], signal[peaks], "x")
+            plt.vlines(x=times[peaks],
+                       ymin=signal[peaks] - properties["prominences"],
+                       ymax = signal[peaks],
+                       color = "C1")
+            plt.hlines(y=properties["width_heights"],
+                       xmin=times[properties["left_ips"]],
+                       xmax=times[properties["right_ips"]],
+                       color = "C1")
+            plt.hlines(y=avgMult * np.mean(signal),
+                       xmin=times[0],
+                       xmax=times[len(signal) - 1],
+                       color = "C2",
+                       linestyles='--')
+            plt.scatter(times[edgesFrame[ch][0]],
+                        signal[edgesFrame[ch][0]],
+                        color="C3")
+            plt.scatter(times[edgesFrame[ch][1]],
+                        signal[edgesFrame[ch][1]],
+                        color="C3")
+            plt.title(f"Ch {ch}")
+            plt.xlabel("Time (s)")
+            plt.ylabel("Signal (V)")
+            plt.show()
+    
+    return edgesFrame
+
+
+def polyBkgFrame(timesFrame,
+                 df,
+                 edgesFrame,
+                 channels,
+                 order=3,
+                 plot=False):
+    r"""
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        A dataframe containing time axis values corresponding to signals in
+        df.
+        
+    df: pandas.core.frame.DataFrame
+        A dataframe of corrected/calibrated dante signal measurements.
+        
+    edgesFrame: pandas.core.frame.DataFrame
+        Dataframe describing edges of the signal containing region, outside
+        of which should be just background. This function will fit to these
+        two early time and late time background containing regions. See
+        signalEdges().
+        
+    channels: list
+        A list of channels for which to apply analyis.
+        
+    order: int
+        Order of polynomial to be fitted to estimated background/hysteresis.
+        
+    plot: bool
+        Flag for plotting fitted background and background subtracted
+        signal. Default is False.
+        
+    Returns
+    -------
+    dfPoly: 
+        
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    dfPoly = pd.DataFrame().reindex_like(df)
+    for ch in channels:
+        # getting polynomial fit of background/hysteresis
+        fitTimes, fitSignals = polyBkg(time=timesFrame[ch],
+                                       signal=df[ch],
+                                       lowerEdge=edgesFrame[ch][0],
+                                       upperEdge=edgesFrame[ch][1],
+                                       order=order,
+                                       lowerLength=None,
+                                       upperLength=None,
+                                       plot=plot)
+        # subtracting hysteresis/background and saving to dataframe
+        dfPoly[ch] = df[ch] - fitSignals
+        
+        if plot:
+            plt.plot(timesFrame[ch], dfPoly[ch])
+            plt.title(f'Channel {ch}')
+            plt.xlabel('Time (s)')
+            plt.ylabel('Signal (V)')
+            plt.show()
+    return dfPoly
+
+
+def highestPeak(signal, peakIdxs):
+    r"""
+    Find the highest peak, and return list of peaks with the highest
+    peak removed from the list.
+    
+    Parameters
+    ----------
+    signal: pandas.core.series.Series
+        A data series consisting of signals from a single dante channel.
+        
+    peakIdxs: list
+        A list of indices corresponding to peaks identified in the signal
+        by using scipy's find_peaks() function.
+    
+    Returns
+    -------
+    peakHighestIdx: int
+        Returns the index corresponding to the highest peak.
+        
+    peakIdxs2: list
+        Returns a list of peak index locations with the highest peak
+        removed from the list. This makes it easier for highestN() to
+        apply highestPeak() iteratively to find the N highest peaks.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # search for first tallest peak
+    peakSignals = signal[peakIdxs]
+    peakHighestIdx = pd.Series.idxmax(peakSignals)
+    # finding index of tallest peak in peakIdxs array (since this is a
+    # numpy array and is indexed differently from peakSignals).
+    delIdx = np.where(peakIdxs==peakHighestIdx)
+    # remove first tallest peak
+    peakIdxs2 = np.delete(peakIdxs, delIdx)
+    return peakHighestIdx, peakIdxs2 
+
+
+def highestN(signal, peakIdxs, peaksNum=2):
+    r"""
+    Select the N tallest peaks.
+    
+    Parameters
+    ----------
+    signal: pandas.core.series.Series
+        A data series consisting of signals from a single dante channel.
+        
+    peakIdxs: list
+        A list of indices corresponding to peaks identified in the signal
+        by using scipy's find_peaks() function.
+        
+    peaksNum: int
+        Number of peaks to grab from peakIdxs. This function will grab
+        just the N tallest peaks where N=peaksNum.
+        
+    Returns
+    -------
+    highestPeaks: numpy.ndarray
+        Array of indices corresponding to the highest peaks. Peaks are
+        ordered from highest to lowest.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    if peaksNum > len(peakIdxs):
+        raise Exception("Number of peaks to grab should be less than or "
+                        "equal to total number of peaks provided!")
+    # initialize input to find highest peak first time around
+    peakIdxsInput = peakIdxs
+    # initialize array containing highest peaks
+    highestPeaks = np.zeros(peaksNum, dtype=int)
+    # loop over finding the N highest peaks
+    for idx in range(peaksNum):
+        # find the highest peak in the given list of peaks
+        peakHighestIdx, peakIdxs2 = highestPeak(signal=signal,
+                                                peakIdxs=peakIdxsInput)
+        # reassign input to be the list of peaks without the most
+        # recently found peak.
+        peakIdxsInput = peakIdxs2
+        # saving highest peak
+        highestPeaks[idx] = peakHighestIdx
+    return highestPeaks
+
+
+def getPeaks(timesFrame,
+             df,
+             channels,
+             peaksNum=2,
+             plot=False,
+             prominence=0.1,
+             width=10,
+             avgMult=1):
+    r"""
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        A dataframe containing time axis values corresponding to signals in
+        df.
+        
+    df: pandas.core.frame.DataFrame
+        A dataframe of corrected/calibrated dante signal measurements.
+    
+    channels: list
+        A list of channels for which to apply analyis.
+        
+    peaksNum: int
+        Number of peaks to grab from peakIdxs. This function will grab
+        just the N tallest peaks where N=peaksNum.
+    
+    peaksNum: int
+        Number of peaks to grab from peakIdxs. This function will grab
+        just the N tallest peaks where N=peaksNum.
+        
+    plot: bool
+        Flag for plotting identified peaks, with prominences, and widths,
+        overlaid with the corresponding dante signal, and the average
+        dante signal.
+        
+    peaksFrame: pandas.core.frame.DataFrame
+        Returns a dataframe containing indices of the identified peaks sorted
+        from the peak that occurs earliest in time to the latest in time.
+        
+    prominence: float
+        Prominence threshold for identifying peaks in scipy's find_peaks().
+        
+    width: int
+        Width in index units for identifying peaks in scipy's find_peaks().
+        
+    avgMult: float
+        Multiplicative factor for setting minimum intensity threshold
+        for indentifying peaks in scipy's find_peaks(). This is a multiple
+        of the signal average.
+        
+    
+    
+    Returns
+    -------
+    peaksFrame:
+        
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    peaksFrame = pd.DataFrame(index=np.arange(peaksNum), columns=df.keys())
+    for ch in channels:
+        times = timesFrame[ch]
+        signal = df[ch]
+        # finding peaks which are above the mean of the signal
+        peaks, properties = find_peaks(signal,
+                                       height=avgMult * np.mean(signal),
+                                       prominence=prominence,
+                                       width=width)
+        # selecting the N tallest peaks
+        highestPeakIdxs = highestN(signal=signal,
+                                   peakIdxs=peaks,
+                                   peaksNum=peaksNum)
+        # ordering peaks by which comes earliest
+        peaksOrdered = np.sort(highestPeakIdxs)
+        # saving these peaks to the dataframe which is to be returned
+        peaksFrame[ch] = peaksOrdered
+        if plot:
+            # plotting
+            plt.plot(times, signal)
+            plt.plot(times[peaks], signal[peaks], "x")
+            plt.vlines(x=times[peaks],
+                       ymin=signal[peaks] - properties["prominences"],
+                       ymax = signal[peaks],
+                       color = "C1")
+            plt.hlines(y=properties["width_heights"],
+                       xmin=times[properties["left_ips"]],
+                       xmax=times[properties["right_ips"]],
+                       color = "C1")
+            plt.hlines(y=avgMult * np.mean(signal),
+                       xmin=times[0],
+                       xmax=times[len(signal) - 1],
+                       color = "C2",
+                       linestyles='--')
+            plt.title(f"Ch {ch}")
+            plt.xlabel("Time (s)")
+            plt.ylabel("Signal (V)")
+            plt.show()
+    return peaksFrame
+
+
+def alignPeaks(timesFrame,
+               df,
+               peaksFrame,
+               channels,
+               peakAlignIdx=0,
+               referenceTime=1e-9,
+               plot=False):
+    r"""
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        A dataframe containing time axis values corresponding to signals in
+        df.
+        
+    df: pandas.core.frame.DataFrame
+        A dataframe of corrected/calibrated dante signal measurements.
+        
+    peaksFrame: pandas.core.frame.DataFrame
+        Dataframe containing positions of N highest peaks and sorted
+        from earliest in time to latest in time. See getPeaks().
+    
+    peakAlignIdx: int
+        Picks which peak to align to. 0 is first peak, 1 is second peak in
+        peaksFrame, etc.
+    
+    referenceTime: float
+        Time in s to which align peaks. Default is 1e-9 s or 1 ns.
+        
+    plot: bool
+        Flag for plotting aligned dante signals. Default is False.
+        
+    Returns
+    -------
+    timesAligned: pandas.core.frame.DataFrame
+        Returns a dataframe identical in shape to timesFrame, but with
+        the times for each dante channel offset such that the selected
+        peaks are temporally aligned.
+    
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    peakAlignIdx = 0 # 0 is first peak, 1 is 2nd peak, etc.
+    referenceTime = 1e-9 # set peaks to occur at 1 ns
+    timesAligned = pd.DataFrame().reindex_like(timesFrame)
+    for ch in channels:
+        time = timesFrame[ch]
+        # of the two highest peaks, grab the one with the lowest index (earliest
+        # in time)
+        # align all signals to reference time, by setting first peak to be at
+        # reference time.
+        # Get index of peak relative to signals dataframe
+        peakIdx = peaksFrame[ch][peakAlignIdx]
+        # get time at which peak occurs
+        peakTime = time[peakIdx]
+        # calculate necessary shift in peak to bring it to the reference time
+        offsetTime = peakTime - referenceTime
+        # generate new time scales for shifted peaks
+#        print(f'times: \n{time}')
+#        print(f'offset time: {offsetTime}')
+        timesAligned[ch] = time - offsetTime
+    # plotting
+    if plot:
+        for ch in channels:
+            plt.plot(timesAligned[ch], df[ch], label=ch)
+        plt.xlabel('Time (s)')
+        plt.ylabel('Signal (V)')
+        plt.title('Aligned')
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.show()
+    return timesAligned
+
+
+def constructMeasurementFrame(timesFrame, df, channels):
+    r"""
+    Takes out put timesFrame and dataFrame from rawProcess.py functions
+    and generates a measurementFrame that can be passed to analyzeStreak()
+    and other main.py functions.
+    
+    Converts units from seconds to nanoseconds.
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        A dataframe containing time axis values corresponding to signals in
+        df.
+        
+    df: pandas.core.frame.DataFrame
+        A dataframe of corrected/calibrated dante signal measurements.
+        
+    channels: list
+        A list of channels for which to apply analyis.
+        
+    Returns
+    -------
+    measurementFrame: pandas.core.frame.DataFrame
+        Returns a measurementFrame which can be passed to main.py functions
+        such as analyzeSpectrum() and analyzeStreak().
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # generating list of column names for measurementFrame
+    colsList= []
+    for ch in channels:
+        timeStr = 'Time' + str(ch)
+        signalStr = 'Signal' + str(ch)
+        colsList.append(timeStr)
+        colsList.append(signalStr)
+    # initialize dataframe
+    measurementFrame = pd.DataFrame(index=np.arange(len(df)), columns=colsList)
+    for ch in channels:
+        # writing data into measurement frame columns
+        measurementFrame['Time' + str(ch)] = timesFrame[ch] * 1e9
+        measurementFrame['Signal' + str(ch)] = df[ch]
+    return measurementFrame
+
+
+def loadCorrected(danteFile,
+                  attenuatorsFile,
+                  offsetsFile,
+                  cut=None,
+                  plot=False,
+                  addCh=[]):
+    r"""
+    Given a dante data file, an attenuators file, and an offsets file, reads
+    the file and applies background correction, attenuation correction, and
+    channel offset correction. Returns the corrected data traces as a pandas
+    dataframe. The row indices of this dataframe also contain the 
+    correct time scaling given the oscilloscope settings, but note that
+    the channels are not aligned. User must apply alignment correction
+    using some measured signal as a temporal fiducial.
+    
+    
+    Parameters
+    ----------
+    danteFile: str
+        Full path to .dat file containing raw dante traces.
+        
+    attenuatorsFile: str
+        Full path to the .xls file containing attenuator serial numbers and
+        corresponding attenuation factors.
+        
+    offsetsFile: str
+        Full path to .xls file containing dante channel offsets.
+        
+    cut: int
+        Number of points to cut from leading and trailing end of each
+        Dante channel trace. This is used to remove noise that occurs
+        at the edges of the signal. Default is None, which means no
+        cut is applied.
+        
+    plot: bool
+        Flag for plotting data after each calibration/correction step.
+        Default is False.
+        
+    addCh: list
+        Add channels to analyze. This is used to override which channels
+        are listed as on in the header of the data dante data file.
+    
+    Returns
+    -------
+    timeOffset: 
+        
+        
+    dfAvg:
+        
+        
+    onChList:
+        
+        
+    hf:
+        
+        
+    dfVolt:
+        
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # load data
+    hf, df = readDanteData(danteFile)
+    # get set of all dante channels with useful data
+    onCh = onChannels(hf)
+    print(f'Analyzing channels {onCh}')
+    # converting from set to list because some functions are implemented
+    # assuming a list of channels instead of a set.
+    if addCh:
+        onChList = list(onCh) + addCh
+    else:
+        onChList = list(onCh)
+    
+    # apply voltage scaling
+    dfVolt, errVolt = voltageScale(hf, df)
+    
+    # time axis calibration
+    timesFrame = timesScope(hf)
+    
+    # shot background correction
+    # Note: Here is where we switch from string indices to integer indices
+    # since we no longer need the bkg columns in our dataframe.
+    timesBkg, dfBkg = bkgCorrect(df=dfVolt,
+                                 timesFrame=timesFrame)
+    
+    # Applying offset correction
+    timeOffset, dfOffset = offsetCorrect(df=dfBkg,
+                                         timesFrame=timesBkg,
+                                         offsetsFile=offsetsFile)
+    
+    # Applying attenuation correction
+    dfAtten = attenuationCorrect(attenuatorsFile, hf, dfOffset, onChList)
+    
+    # applying an averaged background correction
+    dfAvg = avgBkgCorrect(timesFrame=timeOffset,
+                          df=dfAtten,
+                          channels=onCh,
+                          timeLength=1e-9)
+    
+    # plotting each calibration/correction step
+    if plot:
+        plt.plot(df)
+        plt.title('raw dante data, fresh off the load')
+        plt.xlabel('Time (counts)')
+        plt.ylabel('Signal (counts)')
+        plt.show()
+        
+        plt.plot(dfVolt)
+        plt.title('Voltage scaled')
+        plt.xlabel('Time (counts)')
+        plt.ylabel('Signal (V)')
+        plt.show()
+        
+        for ch in onChList:
+            plt.plot(timesFrame[str(ch)], dfVolt[str(ch)], label=ch)
+        plt.title('time and voltage scaled')
+        plt.xlabel('Time (s)')
+        plt.ylabel('Signal (V)')
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.show()
+        
+        for ch in onChList:
+            plt.plot(timesBkg[ch], dfBkg[ch], label=ch)
+        plt.title('bkg corrected')
+        plt.xlabel('Time (s)')
+        plt.ylabel('Signal (V)')
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.show()
+        
+        for ch in onChList:
+            plt.plot(timeOffset[ch], dfOffset[ch], label=ch)
+        plt.title('offset corrected')
+        plt.xlabel('Time (s)')
+        plt.ylabel('Signal (V)')
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.show()
+        
+        for ch in onChList:
+            plt.plot(timeOffset[ch], dfAtten[ch], label=ch)
+        plt.title('attenuation corrected')
+        plt.xlabel('Time (s)')
+        plt.ylabel('Signal (V)')
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.show()
+        
+        for ch in onChList:
+            plt.plot(timeOffset[ch], dfAvg[ch], label=ch)
+        plt.title('averaged background corrected')
+        plt.xlabel('Time (s)')
+        plt.ylabel('Signal (V)')
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.show()
+    
+    return timeOffset, dfAvg, onChList, hf, dfVolt
+
+
+def hysteresisCorrect(timesFrame,
+                      df,
+                      channels,
+                      order=5,
+                      prominence=0.2,
+                      width=10,
+                      avgMult=1):
+    r"""
+    Corrects for hysteresis by detecting edges of signal containing region
+    and fitting a polynomial background to regions that do not belong to 
+    signal. This background is then subtracted.
+    
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        Time corresponding to df
+        
+    df: pandas.core.frame.DataFrame
+        Dataframe of dante signals. See loadCorrected().
+        
+    channels: list
+        A list of channels for which to apply analyis.
+        
+    order: int
+        Polynomial order to be fitted to hysteresis/background.
+        
+    prominence: float
+        Prominence threshold for identifying peaks in scipy's find_peaks().
+        
+    width: int
+        Width in index units for identifying peaks in scipy's find_peaks().
+        
+    avgMult: float
+        Multiplicative factor for setting minimum intensity threshold
+        for indentifying peaks in scipy's find_peaks(). This is a multiple
+        of the signal average.
+        
+        
+    Returns
+    -------
+    dfPoly: pandas.core.frame.DataFrame
+        Returns a dataframe of hysteresis corrected dante signals.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    # getting edges of signal containing region
+    edgesFrame = signalEdges(timesFrame=timesFrame,
+                             df=df,
+                             channels=channels,
+                             plot=False,
+                             prominence=prominence,
+                             width=width,
+                             avgMult=avgMult)
+    
+    # removing hysteresis and background with a polynomial fit
+    dfPoly = polyBkgFrame(timesFrame=timesFrame,
+                          df=df,
+                          edgesFrame=edgesFrame,
+                          channels=channels,
+                          order=order,
+                          plot=False)
+    return dfPoly
+
+
+def align(timesFrame,
+          df,
+          channels,
+          peaksNum=1,
+          peakAlignIdx=0,
+          referenceTime=1e-9,
+          prominence=0.01,
+          width=10,
+          avgMult=1.5):
+    r"""
+    Aligns dante signals based on peak finding.
+    
+    Parameters
+    ----------
+    timesFrame: pandas.core.frame.DataFrame
+        Time corresponding to df
+        
+    df: pandas.core.frame.DataFrame
+        Dataframe of dante signals. See loadCorrected().
+        
+    channels: list
+        A list of channels for which to apply analyis.
+        
+    peaksNum: int
+        Number of peaks to grab from peakIdxs. This function will grab
+        just the N tallest peaks where N=peaksNum.
+        
+    peakAlignIdx: int
+        Picks which peak to align to. 0 is first peak, 1 is second peak in
+        peaksFrame, etc.
+        
+    referenceTime: float
+        Time in s to which align peaks. Default is 1e-9 s or 1 ns.
+        
+    prominence: float
+        Prominence threshold for identifying peaks in scipy's find_peaks().
+        
+    width: int
+        Width in index units for identifying peaks in scipy's find_peaks().
+        
+    avgMult: float
+        Multiplicative factor for setting minimum intensity threshold
+        for indentifying peaks in scipy's find_peaks(). This is a multiple
+        of the signal average.
+        
+    Returns
+    -------
+    timesAligned: pandas.core.frame.DataFrame
+        Returns a dataframe of times corresponding to signals in df, such
+        that the signals are now aligned to the given peak.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    if peakAlignIdx > peaksNum - 1:
+        raise Exception(f"Cannot align to {peakAlignIdx}th peak if there"
+                        f" are only {peaksNum} peaks.")
+    # finding tallest peaks in the signal
+    peaksFrame = getPeaks(timesFrame=timesFrame,
+                          df=df,
+                          channels=channels,
+                          peaksNum=peaksNum,
+                          plot=True,
+                          prominence=prominence,
+                          width=width,
+                          avgMult=avgMult)
+    
+    # aligning to selected peak
+    timesAligned = alignPeaks(timesFrame=timesFrame,
+                              df=df,
+                              peaksFrame=peaksFrame,
+                              channels=channels,
+                              peakAlignIdx=peakAlignIdx,
+                              referenceTime=referenceTime,
+                              plot=True)
+    return timesAligned
+
+#-----------------------------------------
+# PURGATORY FOR DEPRECATED METHOD ATTEMPTS
+#-----------------------------------------
+
+#def signalEdges(time,
+#                signal,
+#                sigma,
+#                windowMin=0,
+#                windowMax=1,
+#                title='',
+#                plot=False):
+#    """
+#    Gets the signal edges via Canny edge detection method
+#    
+#    time: numpy.ndarray
+#        Array of times corresponding to signal.
+#        
+#    signal: numpy.ndarray
+#        Array of signal intensities.
+#        
+#    sigma: float
+#        Width of Gaussian filter used in Canny edge detection. This should
+#        be close to the width of the features you want to detect.
+#        
+#    windowMin: int
+#        Lower bound index of window over which to search for edges in signal.
+#        Default is 0.
+#        
+#    windowMax: int
+#        Distance in index units from end of signal to upper bound of window
+#        over which to search for edges in signal. Default is 1.
+#        
+#    plot: bool
+#        Flag for plotting signal with detected edges. Default is False.
+#    """
+#    timeWin = time[windowMin:-windowMax]
+#    signalWin = signal[windowMin:-windowMax]
+#    # stacking signals to form an "image" so that scikit-image can process it.
+#    sigStack = np.vstack((signalWin, signalWin, signalWin))
+#    edgeCanny = feature.canny(sigStack,
+#                              sigma=sigma)
+#    # getting the actual edge points instead of just the indices
+#    edgeTimes = timeWin[edgeCanny[1,:]]
+#    edges = signalWin[edgeCanny[1,:]]
+#    
+#    if plot:
+#        # plotting detected edges
+#        plt.plot(time, signal, label='signal')
+#        plt.scatter(edgeTimes, edges, color='red', label='edges')
+#        plt.legend()
+#        plt.title(title)
+#        plt.show()
+#    
+#    return edgeTimes, edges
+#
+#
+#
+#
+#
+#def highestTwo(signal, peakIdxs):
+#    """
+#    Get the two highest peaks given signal, and peak
+#    indices found using a peaking finder such as find_peaks_cwt().
+#    """
+#    # search for first tallest peak
+#    peakSignals1 = signal[peakIdxs]
+#    peak1Idx = np.argmax(peakSignals1)
+#    # convert back into index for signals
+#    peak1SigIdx = peakIdxs[peak1Idx]
+#    # remove first tallest peak
+#    peakIdxs2 = np.delete(peakIdxs, peak1Idx)
+#    peakSignals2 = signal[peakIdxs2]
+#    # search for second tallest peak
+#    peak2Idx = np.argmax(peakSignals2)
+#    # convert back into index for signals
+#    peak2SigIdx = peakIdxs2[peak2Idx]
+#    # return positions of two tallest peaks in order of which is first
+#    if peak1SigIdx < peak2SigIdx:
+#        return peak1SigIdx, peak2SigIdx
+#    else:
+#        return peak2SigIdx, peak1SigIdx
\ No newline at end of file
diff --git a/fiducia/response.py b/fiducia/response.py
new file mode 100644
index 0000000..e541e90
--- /dev/null
+++ b/fiducia/response.py
@@ -0,0 +1,94 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:23:02 2019
+
+Utilities for working with DANTE response functions (e.g. plotting, locating
+edges).
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["knotFind"]
+
+
+def knotFind(channels,
+             responseFrame,
+             forceKnot=np.array([]),
+             knotBoundary=0,
+             boundary='y0'):
+    r"""Find knot points.
+    
+    Find knot points for cubic splines based on positions of K-edges
+    of each DANTE channel filter.
+    
+    Parameters
+    ----------
+    channels: list, numpy.ndarray
+        List or array of relevant channels
+    
+    responseFrame: pandas.core.frame.DataFrame
+        Pandas dataFrame containing response functions for each DANTE
+        channel. See loadResponses().
+        
+    forceKnot: numpy.ndarray
+        Numpy array where first column is channelNumber and second column
+        is the corresponding photonEnergy we wish to force. Use this for
+        channels that do not have a distinct K-edge.
+        
+    knotBoundary: float
+        Photon energy value for either y_0 or y_{n+1} boundary condition.
+        This value gets appended to the array of photon energies otherwise
+        found by knotFind().
+        
+    Returns 
+    -------
+    knotsAppend: numpy.ndarray
+        An array of knot points, with each element corresponding to a channel
+        or boundary condition.
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    knots = np.zeros(len(channels))
+    for idx, channel in enumerate(channels):
+        # if forceKnot isn't an empty list, then we go about forcing the user
+        # provided values.
+        if forceKnot.size != 0:
+            if channel in forceKnot[:,0]:
+                # if the channel is forced by the user then put the user provided
+                # photon energy into the knots array
+                forceIdx = np.where(forceKnot[:,0] == channel)
+                knots[idx] = forceKnot[forceIdx, 1]
+            else:
+                # finding largest negative gradient, which should correpond to the
+                # K-edge of the DANTE channel filter.
+                grad = -np.gradient(responseFrame[channel])
+                maxIndex = np.argmax(grad)
+                knots[idx] = responseFrame['Energy(eV)'][maxIndex]
+        else:
+            # finding largest negative gradient, which should correpond to the
+            # K-edge of the DANTE channel filter.
+            grad = -np.gradient(responseFrame[channel])
+            maxIndex = np.argmax(grad)
+            knots[idx] = responseFrame['Energy(eV)'][maxIndex]
+        
+    if boundary == 'y0':
+        knotsAppend = np.append([knotBoundary], knots)
+    elif boundary == 'yn+1':
+        knotsAppend = np.append(knots, [knotBoundary])
+    else:
+        raise Exception(f"No method found for boundary {boundary}.")
+    return knotsAppend
\ No newline at end of file
diff --git a/fiducia/stats.py b/fiducia/stats.py
new file mode 100644
index 0000000..a0f61b8
--- /dev/null
+++ b/fiducia/stats.py
@@ -0,0 +1,400 @@
+ # -*- coding: utf-8 -*-
+"""
+Created on Tue Jul 21 12:11:12 2020
+
+Common statistic operations
+
+@author: Myles T. Brophy
+"""
+# python modules
+import numpy as np
+ 
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["simpsVariance",
+           "trapzVariance",
+           "gradientVariance",
+           "dotVariance",
+           "interpVariance",
+           ]
+
+
+def simpsVariance(yUnc, x=None, dx=1.0):
+    r"""
+    Propagates variance (:math:`\sigma^2`) through Simpson's rule
+    numerical integration.
+    
+    NOTE: THIS FUNCTION IS INCOMPLETE AND HAS NOT BEEN VERIFIED
+    2020-08-21 PMK.
+    
+    Parameters
+    ----------
+    yUnc: list, numpy.ndarray
+        Uncertainties in the vertical axis.
+        
+    x: list, numpy.ndarray, optional
+        Horizontal coordinates corresponding to yUnc. Default is None,
+        which generates a uniformly spaced linear array of horizontal
+        coordinates based on the length of yUnc and the value of dx.
+        
+    dx: float, optional
+        Uniform spacing between horizontal coordinates corresponding to
+        yUnc. Default is 1.0.
+        
+    Returns
+    -------
+    variance: float
+        Variance (:math:`\sigma^2`) on value of integral from Simpson's
+        rule numerical integration.
+        
+    Notes
+    -----
+    Based on a modified verison of:
+        https://en.wikipedia.org/wiki/Simpson's_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    
+    """
+    N = len(yUnc) - 1
+    yUnc = np.asarray(yUnc)
+    if x is None:
+        d = np.full(N, dx)
+    else:
+        d = np.diff(np.asarray(x))
+
+    variance = 0.0
+    #TODO validate with hand calculations
+    for i in range(1, N, 2):
+        hph = d[i] + d[i - 1]
+        variance += yUnc[i]**2 * ( ( d[i]**3 + d[i - 1]**3
+                           + 3. * d[i] * d[i - 1] * hph )\
+                     / ( 6 * d[i] * d[i - 1] ) )**2
+        variance += yUnc[i - 1]**2 * ( ( 2. * d[i - 1]**3 - d[i]**3
+                              + 3. * d[i] * d[i - 1]**2)\
+                     / ( 6 * d[i - 1] * hph) )**2
+        variance += yUnc[i + 1]**2 * ( ( 2. * d[i]**3 - d[i - 1]**3
+                              + 3. * d[i - 1] * d[i]**2)\
+                     / ( 6 * d[i] * hph ) )**2
+
+    if (N + 1) % 2 == 0:
+        variance += yUnc[N]**2 * ( ( 2 * d[N - 1]**2
+                          + 3. * d[N - 2] * d[N - 1])\
+                     / ( 6 * ( d[N - 2] + d[N - 1] )) )**2
+        variance += yUnc[N - 1]**2 * ( ( d[N - 1]**2
+                           + 3*d[N - 1]* d[N - 2] )\
+                     / ( 6 * d[N - 2] ) )**2
+        variance -= yUnc[N - 2]**2 * (d[N - 1]**3\
+                     / ( 6 * d[N - 2] * ( d[N - 2] + d[N - 1] )))**2
+    return variance
+
+
+def trapzVariance(yUnc, x=None, dx=1.0):
+    r"""
+    Error propogation for Trapezoidal rule integration using uniform or non-uniform grids.
+        
+    Parameters
+    ----------
+    yUnc : list, numpy.ndarray
+        The list of uncertainities, referenced as :math:`\sigma_i`.
+        
+    x : list, numpy.ndarray, optional
+        The sampling points for which the uncertainites ''y'' were found. Must be the
+        same length as ''y''. If none are provided, then the step size will be uniform 
+        and set with ''dx''. The default is None.
+    dx : int, float, optional
+        Step size. Only applies if sampling points aren't specified
+        with ''x''. The default is 1.0.
+
+    Returns
+    -------
+    variance : float
+        The total variance (:math:`\sigma^2`) found by propagating ''y''.
+        
+    Notes
+    -----
+    Trap rule integration with non uniform spacing takes the form
+    
+    .. math::
+        \sum_{k=1}^N \frac{\Delta x_i}{2} \left(f(x)_{i-1} + f(y)_i \right)
+    
+    Propogating the uncertainties through this integration results in
+    
+    .. math::
+        \sigma^2 = \frac{1}{4} \left(\sum_{k=1}^N \Delta x_i \sigma_{i-1}^2 + \sigma_i^2 + 2\sum_{k=1}^{N-1} \Delta x_i \Delta x_i+1 \sigma_i^2 \right)
+    
+    The equation is generalized and applies to uniform and non-uniform step sizes.
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    N = len(yUnc)
+    yUnc = np.asarray(yUnc)
+    if x is None:
+        d = np.full(N-1, dx)
+    else:
+        d = np.diff(np.asarray(x))
+    variance = 0.0
+    #add variance of indiviual trapezoid variances
+    for i in range(1, N):
+        variance += d[i-1] ** 2 * 0.25 * (yUnc[i-1] ** 2 + yUnc[i] ** 2) 
+    #add covariant terms because of edge overlap
+    for i in range(1, N-1):
+        variance += 0.5 * d[i-1] * d[i] * yUnc[i] ** 2
+    return variance
+
+
+def gradientVariance(yUnc, x=None, dx=1.0):
+    r"""
+    Propogates uncertainty for the gradient operator of an array of a given step size.
+
+    Parameters
+    ----------
+    yUnc : list, numpy.ndarray
+        The list of uncertainities, referenced as :math:'\sigma_i'.
+        
+    x : list, numpy.ndarray, optional
+        The sampling points for which the uncertainites ''y'' were found. Must be the
+        same length as ''y''. If none are provided, then the step size will be uniform 
+        and set with ''dx''. The default is None.
+        
+    dx : list, numpy.ndarray, optional
+        Step size. Only applies if sampling points aren't specified
+        with ''x''. The default is 1.0.
+
+    Returns
+    -------
+    variance : float
+        The total variance (:math:`\sigma^2`) found by propagating ''y''.
+
+    Notes
+    -----
+    
+    .. math::
+        \operatorname{Var}(\nabla y_i) = \frac{h_{i-1}^2 \sigma_{i+1}^2 + (h_i^2 + h_{i-1}^2)^2 \sigma_i^2 - h_i^4 \sigma_{i-1}^2}{(h_i h_{i-1}(h_i + h_{i-1}))^2}
+
+    At the boundaries
+    
+    .. math::
+        \operatorName{Var}(\nabla y_0) = \frac{\sigma_1^2 - \sigma_0^2}{h_0^2}, \operatorName{Var}(\nabla y_{N-1}) = \frac{\sigma_{N-1}^2 - \sigma_{N-2}^2}{h_{N-2}^2} 
+   
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    N = len(yUnc)
+    np.asarray(yUnc)
+    if x is None:
+        h = np.full(N-1, dx)
+    else:
+        h = np.diff(np.asarray(x))
+    
+    grad = np.zeros(N)
+    #square the uncertainites, cleaner to do it here than in the loop
+    yUnc = yUnc**2
+    grad[0] == (yUnc[1]-yUnc[0])/(h[0]**2)
+    grad[N-1] == (yUnc[N-1]-yUnc[N-2])/(h[N-2]**2)
+    for i in range(1, N-1):
+        firstTerm = h[i-1]**4*yUnc[i+1]
+        secondTerm = (h[i]**2 - h[i-1]**2)**2 * yUnc[i]
+        thirdTerm = h[i]**4*yUnc[i-1]
+        denom = (h[i]*h[i-1]*(h[i]+h[i-1]))**2
+        
+        grad[i] = (firstTerm + secondTerm - thirdTerm) / denom  
+
+    return grad
+
+
+def dotVariance(a, b, aUncertainty=None, bUncertainty=None):
+    r"""
+    Propogate uncertainty for the dot product of matrix a and 1D vector b.
+    
+    Propogate uncertainty for the dot product of a matrix and a 1D vector.
+    Assumes no covariance between `a` and `b`. Methodology is similar to :func:`numpy.dot` where:
+        
+    - If both `a` and `b` are 1D, the uncertainty of the inner product of vectors
+      is returned.
+    - If 'a' is N dimensional (Where :math: `N>=2`) and `b` is 1D, the uncertainty of the sum
+      product of the last axis of a with b is returned.
+    
+    0-D (scalar) arrays are not supported.
+    `b` arrays that have more than one axis are not supported.
+    `a` and `b` must have the same shape as `aUncertainty` and `bUncertainty`,
+    respectively.
+    
+    Parameters
+    ----------
+    a : numpy.ndarray, list
+        Matrix or vector to dot with 'b'.
+  
+    b : numpy.ndarray, list
+        Vector that 'a' will be dotted with. Must be the same size as 
+        the last axis of a.
+   
+    aUncertainty : numpy.ndarray, optional
+        Uncertainty of each element in 'a'. The default is None.
+   
+    bUncertainty : numpy.ndarray, optional
+        Uncertainty of each element in 'b'. The default is None.
+
+    Returns
+    -------
+    variance : float
+    
+    Notes
+    -----
+    
+    .. math::
+        \operatorname{Var}(A \cdot B) = \sum_{i=1}^N \operatorname{Var}(a_i b_i)
+    
+    Assuming covariance between independent variables
+    
+    .. math:: 
+        \sum_{i=1}^N (a_i\sigma_{b_i})^2 + (b_i\sigma_{a_i})^2
+    
+    See also
+    --------
+    
+    Examples
+    --------
+
+    """
+    a = np.asarray(a)
+    b = np.asarray(b)
+    
+    if aUncertainty is None:
+        aUncertainty = np.zeros(a.shape)
+    else:
+        aUncertainty = np.asarray(aUncertainty)
+        
+    if bUncertainty is None:
+        bUncertainty = np.zeros(b.shape)
+    else:
+        bUncertainty = np.asarray(bUncertainty)
+    
+    #first see if the shapes are compatible 
+    #also checks that 
+    try:
+        if b.ndim != 1:
+            raise NotImplementedError("Propogating uncertainty not implemented for b arguments that aren't 1D.")
+        if b.shape != bUncertainty.shape or a.shape != aUncertainty.shape:
+            raise ValueError("Array does not have the same shape as it's uncertainty array.")
+        np.dot(a, b)
+        np.dot(a, bUncertainty)
+        np.dot(aUncertainty, b)
+    except ValueError:
+        raise
+    else:
+        #uncertainty propagation for dotting vectors
+        if a.ndim == 1:
+            variance = np.dot(b ** 2, aUncertainty ** 2) + np.dot(a ** 2, bUncertainty ** 2)
+            return variance
+        #uncertainty propagation for dotting matrix with vector
+        elif a.ndim >=2:
+            variance = np.zeros(a.shape[:-1])
+            for i in range(a.shape[0]):
+                variance[i] = dotVariance(a[i], b, aUncertainty[i], bUncertainty)
+            return variance
+    
+
+def interpVariance(x, xp, fpUnc, leftVar=None, rightVar=None, period=None):
+    r"""
+    Propagate uncertainty for linear interpolation.
+    
+
+    Parameters
+    ----------
+    x : numpy.ndarray, list
+        The x-coordinates at which to evaluate the interpolated values.
+        
+    xp : numpy.ndarray, list
+        The 1D x-coordinates of the data points, must be increasing order
+        
+    fpUnc : numpy.ndarray, list
+        The uncertainty in the y-coordinates of the data points, same length as `xp`.
+        
+    leftVar : float, optional
+        Variance to return for `x < xp[0]`. If not given, the first `yUnc` element 
+        will be used. Default is `None`
+        
+    rightVar : float, optional
+        Variance to return for `x > xp[-1]`. If not given, the last `yUnc` element 
+        will be used. Default is `None`
+
+
+    Returns
+    -------
+    yVar : numpy.ndarray
+        1D array containing the variance for each interpolated `x`.
+
+    Notes
+    -----
+    Variance of interpolated point, assuming no uncertainty in `x` and `xp`,
+    and no covariance between y-coordinates, is given by
+    
+    .. math::
+        \operatorname{Var}(y) = \frac{1}{(x_1 - x_0)^2} ( (x_1 - x)^2 \sigma_{y_0}^2 + (x-x_0)^2 \sigma_{y_1}^2 )
+        
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    #conver to ndarray
+    x = np.asarray(x)
+    xp = np.asarray(xp)
+    fpUnc = np.asarray(fpUnc)
+    
+    #make space for final answer
+    yVar = np.zeros(x.shape)
+    
+    #find the indexs where x is between xp_i and xp_{i+1}
+    indices = np.searchsorted(xp, x, side='left')
+    
+    for i, index in enumerate(indices):
+        #if index is 0 and is not equal to the first data point
+        if index == 0 and x[i] != xp[0]: 
+            #set to default value  
+            if leftVar:
+                variance = leftVar
+            else:
+                variance = fpUnc[0]**2
+        #if the index is after the last datapoint
+        elif index >= len(xp): 
+            #set to default value
+            if rightVar:
+                variance = rightVar
+            else:
+                variance = fpUnc[-1]**2
+        else:
+           rightIndex = index
+           leftIndex = index-1
+           #if we already have this value
+           if xp[rightIndex] == x[i]:
+               
+               variance = fpUnc[rightIndex] ** 2
+           elif xp[leftIndex] == x[i]:
+               variance = fpUnc[leftIndex] ** 2
+           else:
+               #new value, actually need to do interpolation propagation
+               term1 = 1 / (xp[rightIndex] - xp[leftIndex]) ** 2
+               term2 = (xp[rightIndex] - x[i]) ** 2 * fpUnc[leftIndex] ** 2 + (x[i] - xp[leftIndex]) ** 2 *fpUnc[rightIndex] ** 2 
+               variance = term1 * term2
+
+        yVar[i] = variance
+        
+    return yVar
+        
+        
+    
\ No newline at end of file
diff --git a/fiducia/visualization.py b/fiducia/visualization.py
new file mode 100644
index 0000000..487216b
--- /dev/null
+++ b/fiducia/visualization.py
@@ -0,0 +1,352 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 10:49:17 2019
+
+Utilities for visualizing DANTE data.
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import matplotlib.pyplot as plt
+
+# custom modules
+import fiducia.pltDefaults
+
+
+# listing all functions declared in this file so that sphinx-automodapi
+# correctly documents them and doesn't document imported functions.
+__all__ = ["plotResponse",
+           "plotTraces",
+           "plotStreak",
+           "signalImg"]
+
+
+def plotResponse(channels,
+                 responseFrame,
+                 knots,
+                 solid=True,
+                 title='Dante Response Functions'):
+    r"""
+    Plots response function curves with knot locations identified as
+    vertical dashed lines.
+    
+    channels: list, numpy.ndarray
+        List or array of relevant channels
+    
+    responseFrame: pandas.core.frame.DataFrame
+        Pandas dataFrame containing response functions for each DANTE
+        channel. See loadResponses().
+        
+    knots: list, numpy.ndarray
+        List or array of knot point photon energy value. See knotFind().
+        
+    solid: Bool
+        Includes solid angle in response function value if true.
+        Necessary for plotting responses with correct units.
+        
+    Parameters
+    ----------
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+        
+    """
+    # getting sensible figure bounds
+    yMax = np.max(np.max(responseFrame[channels]))
+    yMin = yMax / 1e11
+    xMin = 1e1
+    xMax = 1e5
+    # plotting figure
+    fig = plt.figure()
+    for idx, channel in enumerate(channels):
+        plt.loglog(responseFrame['Energy(eV)'],
+                   responseFrame[channel],
+                   label=str(channel))
+    for idx, _ in enumerate(knots):
+        plt.loglog((knots[idx], knots[idx]), (yMin, yMax), '--', color='grey')
+    plt.xlim((xMin, xMax))
+    plt.ylim((yMin, yMax))
+    plt.legend(frameon=False,
+               labelspacing=0.001,
+               borderaxespad=0.1)
+    plt.xlabel('Energy (eV)')
+    if solid:
+        plt.ylabel('Response (V/GW/sr)')
+    else:
+        plt.ylabel('Response (V/GW)')        
+    plt.title(title)
+    plt.show()
+    return
+
+
+def plotTraces(channels, measurementFrame, scale='regular'):
+    r"""
+    Given a dataframe of Dante channel data, plot all the signal
+    traces onto a single plot.
+    
+    channels: list, numpy.ndarray
+        List or array of relevant channels
+        
+    measurementFrame: pandas.core.frame.DataFrame
+        Pandas dataframe containing DANTE measurement data. See
+        readDanteData() and readDanProcessed().
+        
+    Parameters
+    ----------
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    if scale == 'regular':
+        for channel in channels:
+            plt.plot(measurementFrame['Time' + str(channel)],
+                     measurementFrame['Signal' + str(channel)],
+                     label=str(channel))
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.xlabel('Time (ns)')
+        plt.ylabel('Signal (V)')
+        plt.title('Dante Measurement Data')
+        plt.show()
+    elif scale == 'log':
+        for channel in channels:
+            plt.semilogy(measurementFrame['Time' + str(channel)],
+                         measurementFrame['Signal' + str(channel)],
+                         label=str(channel))
+        plt.legend(frameon=False,
+                   labelspacing=0.001,
+                   borderaxespad=0.1)
+        plt.xlabel('Time (ns)')
+        plt.ylabel('Signal (V)')
+        plt.title('Dante Measurement Data')
+        plt.show()
+    else:
+        raise Exception(f"No such method for scale: {scale}")
+    return
+
+
+def uniformStreak(times, energies, spectra):
+    r"""
+    Convert streak spectra from an irregularly spaced energy sampling
+    to uniformly spaced energies via linear interpolation. This is
+    primarily useful for visualizing the streaked spectrum in plotStreak().
+    
+    times : numpy.ndarray
+        1D array of times corresponding to streaked spectrum
+        
+    energies : numpy.ndarray
+        2D array of energies which are nonuniformly spaced. Each column
+        corresponds to the photon energy axis for a particularly time.
+        Typically all the columns are identical.
+        
+    spectra : numpy.ndarray
+        2D array of streaked spectra corresponding to energies array.
+    """
+    # generate a uniform energy axis
+    energyMin = int(np.min(energies)) # eV
+    energyMax = int(np.max(energies)) # eV
+    energyStep = 1 # eV
+    energyNew = np.arange(start=energyMin, stop=energyMax, step=energyStep)
+    newShape = (len(energyNew), len(times))
+    energiesNew = np.zeros(newShape)
+    spectraNew = np.zeros(newShape)
+    for idx, time in enumerate(times):
+        energyOld = energies[:, idx]
+        spectrumOld = spectra[:, idx]
+        spectrumNew = np.interp(x=energyNew, xp=energyOld, fp=spectrumOld)
+        # storing values
+        energiesNew[:, idx] = energyNew
+        spectraNew[:, idx] = spectrumNew
+    return energiesNew, spectraNew
+
+
+def plotStreak(times, energies, spectra):
+    r"""
+    Plot streak of unfolded Dante spectra. See analyzeStreak().
+    
+    Parameters
+    ----------
+    times: numpy.ndarray
+        Array of times for which the unfold was analyzed.
+        
+    energies: numpy.ndarray
+        Array of photon energies corresponding to the unfolded spectra.
+        
+    spectra: numpy.ndarray
+        The unfolded spectral intensities as a 2D array. See analyzeStreak().
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    energiesInterp, spectraInterp = uniformStreak(times, energies, spectra)
+    extent=[energiesInterp[0,0],
+            energiesInterp[-1,0],
+            times[0],
+            times[-1]]
+    
+    plt.imshow(spectraInterp.transpose(),
+               extent=extent,
+               origin='lower',
+               aspect='auto')
+    plt.ylabel('Time (ns)')
+    plt.xlabel('Photon Energy (eV)')
+    cbar = plt.colorbar()
+    cbar.set_label("Spectrum (GW/sr/eV)")
+    plt.show()
+    return
+
+
+def signalImg(signalsArr):
+    r"""
+    Visualize dante signals as an image.
+    
+    Parameters
+    ----------
+    signalsArr: numpy.ndarray
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    plt.imshow(signalsArr.transpose(), aspect='auto')
+    plt.xlabel('Time')
+    plt.ylabel('channel')
+    cbar = plt.colorbar()
+    cbar.set_label('Signal')
+    plt.show()
+    return
+
+
+def individualChPlot(channels, timesFrame, signalsFrame, cut=None):
+    r"""
+    Plot individual channels on separate plots.
+    
+    Parameters
+    ----------
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    if cut:
+        timesFrame = timesFrame.iloc[cut:-cut]
+        signalsFrame = signalsFrame.iloc[cut:-cut]
+    for chNum in channels:
+        plt.plot(timesFrame[chNum] * 1e9, signalsFrame[chNum])
+        plt.title(f"Ch {chNum}")
+        plt.xlabel("Time  (ns)")
+        plt.ylabel("Signal (V)")
+        plt.show()
+    return
+
+
+def traceGrid(channels,
+              timesFrame,
+              signalsFrame,
+              chStatus=None,
+              saveName="",
+              cut=None,
+              shotNum=""):
+    r"""
+    Plot individual channels on separate plots arranged into a grid.
+    
+    Parameters
+    ----------
+    
+    Returns
+    -------
+        
+    Notes
+    -----
+    
+    See also
+    --------
+    
+    Examples
+    --------
+    """
+    # cutting leading and trailing edges of signal where noise exists
+    if cut:
+        timesFrame = timesFrame.iloc[cut:-cut]
+        signalsFrame = signalsFrame.iloc[cut:-cut]
+        
+    # setting up subplots for all 18 channels
+    fig, ax = plt.subplots(6, 3, sharex=True, sharey=True, figsize=(15,10))
+    fig.subplots_adjust(hspace=0, wspace=0)
+    
+    # creating array of all possible dante channels
+    allCh = np.arange(18) + 1
+    # annotating each subplot window with the channel number
+    for chNum in allCh:
+        col = int((chNum - 1)/6)
+        row = (chNum - 1) % 6
+        if not chStatus:
+            ax[row][col].annotate(f'Ch {chNum}', xy=(0, 0.8), xytext=(0, 0.8))
+        else:
+            ax[row][col].annotate(f'Ch {chNum} - {chStatus[chNum]}', xy=(0, 0.8), xytext=(0, 0.8))
+    # plotting traces for only user specified channels
+    for chNum in channels:
+        col = int((chNum - 1)/6)
+        row = (chNum - 1) % 6
+        ax[row][col].plot(timesFrame[chNum] * 1e9,
+                          signalsFrame[chNum] / np.max(signalsFrame[chNum]))
+        ax[row][col].set_ylim((-0.1, 1.1))
+#        ax[row][col].annotate(f'Ch {chNum}', xy=(0, 0.8), xytext=(0, 0.8))
+    fig.suptitle(f"Shot {shotNum}")
+    fig.text(0.5, 0.04, 'Time (ns, arbitrary offset)', ha='center')
+    fig.text(0.04, 0.5, 'Signal (normalized to ch max)', va='center', rotation='vertical')
+    
+    if saveName:
+        # if there is a non-empty name, then try to save
+        plt.savefig(saveName)
+    else:
+        # no name, then just show
+        plt.show()
+    return
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..51ef4ac
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,6 @@
+numpy (>= 1.15.0)
+matplotlib (>= 3.0.0)
+pandas (>= 0.23.0)
+astropy (>= 3.1)
+scipy (>= 1.1.0)
+xarray (>=0.15.1)
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..8796ee3
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Feb 26 16:11:09 2019
+
+Fiducia setup script.
+
+@author: Pawel M. Kozlowski
+"""
+
+from setuptools import setup
+
+setup(name='fiducia',
+      version='0.2',
+      description='DANTE data unfolder using cubic-spline algorithm',
+      url='https://github.com/lanl/fiducia',
+      author='Pawel M. Kozlowski, et al.',
+      author_email='pkozlowski@lanl.gov',
+      license='BSD',
+      packages=['fiducia'],
+      zip_safe=True)
\ No newline at end of file


Test patch:
diff --git a/fiducia/tests/__init__.py b/fiducia/tests/__init__.py
new file mode 100644
index 0000000..f97aad1
--- /dev/null
+++ b/fiducia/tests/__init__.py
@@ -0,0 +1,8 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 25 13:14:52 2019
+
+@author: Pawel M. Kozlowski
+"""
+
diff --git a/fiducia/tests/conftest.py b/fiducia/tests/conftest.py
new file mode 100644
index 0000000..645faa2
--- /dev/null
+++ b/fiducia/tests/conftest.py
@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 21 18:10:57 2020
+
+Place for setting up fixtures to be used across pytests.
+Loading example dante datasets should be done here so that they
+don't have to be repeated across the various test modules.
+
+@author: Pawel M. Kozlowski
+"""
+
+import pytest
+import numpy as np
+
+
+@pytest.fixture(scope="session")
+def example_info(tmpdir_factory):
+    r"""
+    Loads example information on Dante data (calibration files, and raw data)
+    from Omega-60, which can be used across multiple pytests.
+    """
+    # response function curves for Dante channels.
+    responseFile = 'example_response.csv'
+    # attenuator information loopkup by ID number.
+    attenuatorsFile = 'TableAttenuators.xls'
+    # Dante oscilloscope voltage axis offsets.
+    offsetsFile = 'Offset.xls'
+    # Raw Dante data file from Omega-60.
+    dataFile = 'dante_dante86455.dat'
+    # emitting area of blackbody in mm^2.
+    area = np.pi * 0.6 ** 2
+    # dante viewing angle relative to normal of emitting surface, in degrees.
+    angle = 69.09
+    # collecting info into a dict which can be passed to various pytests.
+    info = {"Response File" : responseFile,
+            "Attenuators File" : attenuatorsFile,
+            "Offsets File" : offsetsFile,
+            "Data File" : dataFile,
+            "Area" : area,
+            "Angle" : angle}
+    return info
+
+
+
+@pytest.fixture(scope="session")
+def load_example_data(example_info):
+    r"""
+    Loads the example Dante data from example_info() so that it only
+    has to be done once across the session.
+    """
+    # load raw dante data into corrected/aligned format
+    from fiducia.rawProcess import loadCorrected
+    loaded = loadCorrected(danteFile=example_info["Data File"],
+                           attenuatorsFile=example_info["Attenuators File"],
+                           offsetsFile=example_info["Offsets File"],
+                           plot=False)
+    timesFrame, dfAtten, onChList, hf, dfVolt = loaded
+    dataDict = {"Times Frame" : timesFrame,
+                "Attenuators Frame" : dfAtten,
+                "On Channels" : onChList,
+                "Header Frame" : hf,
+                "Signal Frame" : dfVolt}
+    return dataDict
\ No newline at end of file
diff --git a/fiducia/tests/test_conftest.py b/fiducia/tests/test_conftest.py
new file mode 100644
index 0000000..5a5bc2e
--- /dev/null
+++ b/fiducia/tests/test_conftest.py
@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 21 19:17:24 2020
+
+A couple of tests to ensure that fixtures from conftest.py are loading
+correctly for the pytest sessions, and that the returned data are
+behaving as expected.
+
+@author: Pawel M. Kozlowski
+"""
+
+import pytest
+
+def test_info_fixture(example_info):
+    """Testing whether conftest.py fixture is working as expected."""
+    assert example_info["Angle"] == 69.09
+    
+
+def test_data_fixture(load_example_data):
+    """Testing whether conftest.py fixture is working as expected."""
+    expectedChOn = [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 17]
+    testTrue = load_example_data["On Channels"] == expectedChOn
+    errStr = (f"Header loaded incorrectly. Expected on channels are "
+              f"{expectedChOn}, but got {load_example_data['On Channels']} "
+              "instead.")
+    assert testTrue, errStr
\ No newline at end of file
diff --git a/fiducia/tests/test_cspline.py b/fiducia/tests/test_cspline.py
new file mode 100644
index 0000000..20a7702
--- /dev/null
+++ b/fiducia/tests/test_cspline.py
@@ -0,0 +1,270 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:42:24 2019
+
+Tests for cspline.py
+
+@author: Pawel M. Kozlowski, Myles T. Brophy
+"""
+
+# python modules
+import numpy as np
+import pytest
+import os
+import inspect
+
+# custom modules
+import fiducia
+import fiducia.cspline as cspline
+import fiducia.loader as loader
+import scipy.sparse as sparse
+
+class Test_splineCoords(object):
+    """
+    Testing class for the splineCoords and splineCorrdsInv functions
+    """
+    def setup_class(self):
+        """Initialize parameters for tests."""
+        self.energyMin = 10
+        self.energyMax = 20
+        self.values1 = np.arange(9) + 1
+        self.energy1 = self.values1 + self.energyMin
+        self.known1 = self.values1 / 10
+            
+    def test_splineCoords(self):
+        """
+        Testing for expected return value.
+        """
+        methodVals = cspline.splineCoords(self.energy1, self.energyMin, self.energyMax)
+        testTrue = np.allclose(methodVals, self.known1, rtol=1e-16, atol=0.0)
+        errStr = (f"Spline coordinate value test gives {methodVals} but "
+                  f"should be {self.known1}")
+        assert testTrue, errStr
+        
+    def test_splineCoordsInv(self):
+        """Testing for expected return value."""
+        methodVals = cspline.splineCoordsInv(self.known1, self.energyMin, self.energyMax)
+        testTrue = np.allclose(methodVals, self.energy1,
+                               rtol=1e-16, atol=0.0)
+        errStr = (f"Spline coordinate inverted value test gives {methodVals} but "
+                  f"should be {self.energy1}")
+        assert testTrue, errStr
+
+class Test_coeffArr(object):
+    """
+    Testing class for the yCoeffArr and dCoeffArr functions
+    """
+    def setup_class(self):
+        """Initialize parameters for tests."""
+        self.energyNorm = 0.5;
+        self.chLen = 9
+        self.trueYCoeffArr = sparse.csc_matrix([
+            [0.875,0.125,0.,0.,0.,0.,0.,0.,0.,0.],
+            [0.,0.5,0.5,0.,0.,0.,0.,0.,0.,0.],
+            [0.,0.,0.5,0.5,0.,0.,0.,0.,0.,0.],
+            [0.,0.,0.,0.5,0.5,0.,0.,0.,0.,0.],
+            [0.,0.,0.,0.,0.5,0.5,0.,0.,0.,0.],
+            [0.,0.,0.,0.,0.,0.5,0.5,0.,0.,0.],
+            [0.,0.,0.,0.,0.,0.,0.5,0.5,0.,0.],
+            [0.,0.,0.,0.,0.,0.,0.,0.5,0.5,0.],
+            [0.,0.,0.,0.,0.,0.,0.,0.,0.875,0.125]])
+        self.trueDCoeffArr = sparse.csc_matrix([
+            [0.375,-0.125,0.,0.,0.,0.,0.,0.,0.,0.],
+            [0.,0.125,-0.125,0.,0.,0.,0.,0.,0.,0.],
+            [0.,0.,0.125,-0.125,0.,0.,0.,0.,0.,0.],
+            [0.,0.,0.,0.125,-0.125,0.,0.,0.,0.,0.],
+            [0.,0.,0.,0.,0.125,-0.125,0.,0.,0.,0.],
+            [0.,0.,0.,0.,0.,0.125,-0.125,0.,0.,0.],
+            [0.,0.,0.,0.,0.,0.,0.125,-0.125,0.,0.],
+            [0.,0.,0.,0.,0.,0.,0.,0.125,-0.125,0.],
+            [0.,0.,0.,0.,0.,0.,0.,0.,0.25,0.125]])
+        
+        
+    def test_yCoeffArr(self):
+        """Testing for expected value"""
+        methodVals = cspline.yCoeffArr(self.energyNorm, self.chLen)
+        testTrue = np.allclose(methodVals.toarray(), 
+                               self.trueYCoeffArr.toarray(),
+                               rtol=1e-16, atol=0.0)
+        errStr = (f"Y coefficient array test gives {methodVals.toarray()} but "
+                  f"should be {self.trueYCoeffArr.toarray()}")
+        assert testTrue, errStr
+        
+    def test_dCoeffArr(self):
+        """Testing for expected value"""
+        methodVals = cspline.dCoeffArr(self.energyNorm, self.chLen)
+        testTrue = np.allclose(methodVals.toarray(),
+                               self.trueDCoeffArr.toarray(),
+                               rtol=1e-16, atol=0.0)
+        errStr = (f"D coefficient array test gives {methodVals.toarray()} but "
+                  f"should be {self.trueDCoeffArr.toarray()}")
+        assert testTrue, errStr
+
+
+class Test_chArr(object):
+    """Testing class for the __chi1Arr__ and __chi3Arr__ functions"""
+    def setup_class(self):
+        """Initialize parameters for test"""
+        self.chLen = 6
+        self.trueChi1 = sparse.csc_matrix([
+            [2.,1.,0.,0.,0.,0.,0.],
+            [1.,4.,1.,0.,0.,0.,0.],
+            [0.,1.,4.,1.,0.,0.,0.],
+            [0.,0.,1.,4.,1.,0.,0.],
+            [0.,0.,0.,1.,4.,1.,0.],
+            [0.,0.,0.,0.,1.,4.,1.],
+            [0.,0.,0.,0.,0.,1.,2.]])
+        self.trueChi3 = sparse.csc_matrix([
+            [-1.,1.,0.,0.,0.,0.,0.],
+            [-1.,0.,1.,0.,0.,0.,0.],
+            [0.,-1.,0.,1.,0.,0.,0.],
+            [0.,0.,-1.,0.,1.,0.,0.],
+            [0.,0.,0.,-1.,0.,1.,0.],
+            [0.,0.,0.,0.,-1.,0.,1.],
+            [0.,0.,0.,0.,0.,-1.,1.]])
+    
+    def test___chi1Arr__(self):
+        """Testing for expected return value."""
+        methodVals = cspline.__chi1Arr__(self.chLen)
+        testTrue = np.allclose(methodVals.toarray(), self.trueChi1.toarray(),
+                               rtol=1e-16, atol=0.0)
+        errStr = (f"Chi1 array test gives {methodVals.toarray()} but "
+                  f"should be {self.trueChi1.toarray()}")
+        assert testTrue, errStr
+        
+    def test___chi3Arr__(self):
+        """Testing for expected return value."""
+        methodVals = cspline.__chi3Arr__(self.chLen)
+        testTrue = np.allclose(methodVals.toarray(), self.trueChi3.toarray(),
+                               rtol=1e-16, atol=0.0)
+        errStr = (f"Chi3 array test gives {methodVals.toarray()} but "
+                  f"should be {self.trueChi3.toarray()}")
+        assert testTrue, errStr
+
+class Test_dToyArr(object):
+    """Testing class for the dToyArr function"""
+    def setup_class(self):
+        self.chLen = 4
+        self.trueDToyArr = [[-1.26785714,  1.60714286, -0.42857143,  0.10714286, -0.01785714],
+       [-0.46428571, -0.21428571,  0.85714286, -0.21428571,  0.03571429],
+       [ 0.125     , -0.75      ,  0.        ,  0.75      , -0.125     ],
+       [-0.03571429,  0.21428571, -0.85714286,  0.21428571,  0.46428571],
+       [ 0.01785714, -0.10714286,  0.42857143, -1.60714286,  1.26785714]]
+        
+    def test_dToyArr(self):
+        """Testing for expected return value."""
+        methodVals = cspline.dToyArr(self.chLen)
+        testTrue = np.allclose(methodVals, self.trueDToyArr, 
+                               rtol=1e-6, atol=0.0)
+        errStr = (f"D to Y test gives {methodVals} but "
+                  f"should be {self.trueDToyArr}")
+        assert testTrue, errStr
+        
+class Test_responseInterp(object):
+    """
+    Testing class for the responseInterp function.
+    """
+    def setup_class(self):
+        "Initialize parameters for testing"
+        self.energyNorm = 0.2
+        self.energyMin = 10
+        self.energyMax = 80
+        #Have to set path this way in order to work locally on pytest and
+        #on the GitLab pipeline. Not pretty.
+        self.fiduciaPath = os.path.abspath(os.path.dirname(inspect.getfile(fiducia)))
+        self.responseFilePath = os.path.join(self.fiduciaPath, "tests", "example_response.csv")
+        
+        self.channels= [2, 3, 4, 5, 6, 7, 8, 9, 12, 13, 14]
+        self.responseFrame = loader.loadResponses(self.channels, self.responseFilePath)
+        self.trueResponsesInterpolated = [[8.19720e-39,0.00000e+00,0.00000e+00,0.00000e+00,2.82800e-38,
+                                           7.03465e-25,2.07858e-09,6.20039e-08,0.00000e+00,0.00000e+00,
+                                           0.00000e+00]]
+        
+    def test_responseInterp(self):
+        """Testing for expected return value."""
+        methodVals = cspline.responseInterp(self.energyNorm, self.energyMin, 
+                                            self.energyMax, self.responseFrame,
+                                            self.channels)
+        testTrue = np.allclose(methodVals, self.trueResponsesInterpolated, 
+                               rtol=1e-13, atol=0.0)
+        errStr = (f"Response interpolation test gives {methodVals} but "
+                  f"should be {self.trueResponsesInterpolated}")
+        assert testTrue, errStr
+    
+class Test_yChiCoeff(object):
+    """
+    Testing class for the yChiCoeffArr and yChiCoeffArrEnergies functions
+    """
+    def setup_class(self):
+        "Initialize parameters for testing"
+        self.chLen = 4
+        self.energyNorm = 0.2
+        self.energyNorms = [0.2,0.3,0.5, 0.9]
+        self.dToYArr = [[-1.26785714,1.60714286,-0.42857143,0.10714286,-0.01785714],
+                        [-0.46428571,-0.21428571,0.85714286,-0.21428571,0.03571429],
+                        [0.125,-0.75,0.,0.75,-0.125],
+                        [-0.03571429,0.21428571,-0.85714286,0.21428571,0.46428571],
+                        [0.01785714,-0.10714286,0.42857143,-1.60714286,1.26785714]]
+    
+        self.trueyChiCoeffArr = [[0.76342857,0.32342857,-0.10971429,0.02742857,-0.00457143],
+                                 [-0.06342857,0.89257143,0.21371429,-0.05142857,0.00857143],
+                                 [0.01714286,-0.10285714,0.92342857,0.19314286,-0.03085714],
+                                 [-0.00514286,0.03085714,-0.12342857,0.97485714,0.12285714]]
+        self.trueyChiCoeffArrEnergies = [[[0.76342857,0.32342857,-0.10971429,0.02742857,
+                                           -0.00457143],
+                                          [-0.06342857,0.89257143,0.21371429,-0.05142857,
+                                           0.00857143],
+                                          [0.01714286,-0.10285714,0.92342857,0.19314286,
+                                           -0.03085714],
+                                          [-0.00514286,0.03085714,-0.12342857,0.97485714,
+                                           0.12285714]],
+
+                                         [[0.656125,0.47925,-0.171,0.04275,
+                                           -0.007125],
+                                          [-0.076125,0.79975,0.342,-0.07875,
+                                           0.013125],
+                                          [0.020625,-0.12375,0.838,0.31275,
+                                           -0.047625],
+                                          [-0.006375,0.03825,-0.153,0.91675,
+                                           0.204375]],
+
+                                         [[0.45758929,0.75446429,-0.26785714,0.06696429,
+                                           -0.01116071],
+                                          [-0.07366071,0.56696429,0.60714286,-0.12053571,
+                                           0.02008929],
+                                          [0.02008929,-0.12053571,0.60714286,0.56696429,
+                                           -0.07366071],
+                                          [-0.00669643,0.04017857,-0.16071429,0.72767857,
+                                           0.39955357]],
+
+                                        [[0.09180357,1.02117857,-0.14271429,0.03567857,
+                                          -0.00594643],
+                                         [-0.01430357,0.08682143,0.97971429,-0.06267857,
+                                          0.01044643],
+                                         [0.00401786,-0.02410714,0.09742857,0.96139286,
+                                          -0.03873214],
+                                         [-0.00176786,0.01060714,-0.04242857,0.16010714,
+                                          0.87348214]]]
+        
+    def test_yChiCoeffArr(self):
+        """Testing for expected return value."""
+        methodVals = cspline.yChiCoeffArr(self.energyNorm, self.chLen, self.dToYArr)
+        testTrue = np.allclose(methodVals, self.trueyChiCoeffArr, 
+                               rtol=1e-6, atol=0.0)
+        errStr = (f"y Chi Coefficient Array test gives {methodVals} but "
+                  f"should be {self.trueyChiCoeffArr}")
+        assert testTrue, errStr
+    
+    def test_yChiCoeffArrEnergies(self):
+        """Testing for expected return value."""
+        methodVals = cspline.yChiCoeffArrEnergies(self.energyNorms, self.chLen, self.dToYArr)
+        testTrue = np.allclose(methodVals, self.trueyChiCoeffArrEnergies, 
+                               rtol=1e-5, atol=0.0)
+        errStr = (f"y Chi Coefficient Energies Array test gives {methodVals} but "
+                  f"should be {self.trueyChiCoeffArrEnergies}")
+        assert testTrue, errStr
+        
+        
+        
+        
\ No newline at end of file
diff --git a/fiducia/tests/test_loader.py b/fiducia/tests/test_loader.py
new file mode 100644
index 0000000..5ff42f9
--- /dev/null
+++ b/fiducia/tests/test_loader.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:21:18 2019
+
+Tests for loader.py
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import pytest
+
+# custom modules
+import fiducia.loader as loader
\ No newline at end of file
diff --git a/fiducia/tests/test_main.py b/fiducia/tests/test_main.py
new file mode 100644
index 0000000..7f110d6
--- /dev/null
+++ b/fiducia/tests/test_main.py
@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 25 12:22:22 2019
+
+Tests for main.py
+
+Tests should be split and marked as:
+    -type check tests, which check function execution and error handling
+    against inputs of different types.
+    - numerical tests, which check whether given input values produce
+    expected output values.
+    - property tests, which check mathematical properties of functions,
+    such as associativity, commutativity, identity, etc.
+    - physical tests, which check physical intution, such as higher power
+    should produce larger radiation temperature, or larger uncertainties
+    on the inputs should produce larger uncertainties on the outputs.
+
+
+@author: Pawel M. Kozlowski, Myles T. Brophy
+"""
+
+# python modules
+import numpy as np
+import pytest
+from hypothesis import given
+import hypothesis.strategies as s
+import hypothesis.extra.numpy as hnp
+import hypothesis.extra.pandas as hpd
+
+# custom modules
+import fiducia.main as main
+
+
+# categorize this test as part of the group of tests which check
+# expected values of functions.
+@pytest.mark.value
+# parametrize the test to run multiple times over different inputs.
+@pytest.mark.parametrize("power,area,angle,expectedTRad",
+                         [(0, 0.6, 30, 0),
+                          (1, 0.6, 30, 36.98703164858843),
+                          (10, 0.6, 30, 65.7732768191428),
+                          (100, 0.6, 30, 116.963263898272)])
+def test_inferRadTemp_val(power, area, angle, expectedTRad):
+    """Testing for expected return value."""
+    methodVals, methodValsVariance = main.inferRadTemp(power, 
+                                                       area,
+                                                       angle)
+    testTrue = np.allclose(methodVals,
+                           expectedTRad, 
+                           rtol=1e-7,
+                           atol=0.0)
+    errStr = (f"Infer radiation temperture test gives {methodVals} but "
+              f"should be {expectedTRad}")
+    assert testTrue, errStr
+        
+
+@given(power=s.floats(min_value=0),
+       powerUncertainty=s.floats(min_value=0))
+def test_inferRadTemp_prop(power, powerUncertainty):
+    """
+    Property-based testing of inferRadTemp() using hypothesis.
+    Given a physical value of power (>=0),
+    expect radiation temperature to be phyiscal (>=0).
+    """
+    area = 0.6
+    angle = 30
+    methodVals, methodValsVariance = main.inferRadTemp(power=power, 
+                                                       area=area,
+                                                       angle=angle,
+                                                       powerUncertainty=powerUncertainty)
+    assert methodVals >= 0 and methodValsVariance >= 0
\ No newline at end of file
diff --git a/fiducia/tests/test_misc.py b/fiducia/tests/test_misc.py
new file mode 100644
index 0000000..9265b07
--- /dev/null
+++ b/fiducia/tests/test_misc.py
@@ -0,0 +1,25 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:25:23 2019
+
+Tests for misc.py
+
+@author: Pawel M. Kozlowski, Myles T. Brophy
+"""
+
+# python modules
+import numpy as np
+import pytest
+
+# custom modules
+import fiducia.misc as misc
+
+
+def test_find_nearest():
+    """
+    Testing for expected return value.
+    """
+    ex = [1,10,100,500,1000];
+    assert misc.find_nearest(ex, 40) == (1,10)
+    assert misc.find_nearest(ex, 5000) == (4,1000)
diff --git a/fiducia/tests/test_response.py b/fiducia/tests/test_response.py
new file mode 100644
index 0000000..82d7e53
--- /dev/null
+++ b/fiducia/tests/test_response.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 09:23:30 2019
+
+Tests for response.py
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import pytest
+
+# custom modules
+import fiducia.response as response
\ No newline at end of file
diff --git a/fiducia/tests/test_stats.py b/fiducia/tests/test_stats.py
new file mode 100644
index 0000000..5fcd10b
--- /dev/null
+++ b/fiducia/tests/test_stats.py
@@ -0,0 +1,43 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jul 21 12:16:31 2020
+
+Testing for stats.py
+
+@author: Myles T. Brophy
+"""
+
+# python modules
+import numpy as np
+import pytest
+
+# custom modules
+import fiducia.stats as stats
+
+class Test_integrationError(object):
+    def setup_class(self):
+        """Initialize parameters for tests."""
+        self.rectangle = np.full(2,5)
+        self.rectangleHeightUncertainty = np.full(2, 1) 
+        
+        self.diagonal = np.arange(5)*3 + 1
+    
+    def test_trapz(self):
+        #check that the trap integration of a rectangle equals width*height
+        assert np.trapz(self.rectangle) == (len(self.rectangle) -1)*self.rectangle[0]
+        #check that the trap integration of a diagonal equals the area of a right trapezoid.
+        assert np.trapz(self.diagonal) == 0.5 * (len(self.diagonal)-1) * (self.diagonal[0] + self.diagonal[-1])
+    
+    def test_trapzUncertainty(self):
+        """Testing for expected return value."""
+        assert stats.trapzVariance([1,2,3]) == 6.5
+        assert stats.trapzVariance([1,2,3], x=[5,7,10]) == 46.25
+        assert stats.trapzVariance([1,2], x=[1,3]) == 5.0 
+        assert stats.trapzVariance([1], x=[2]) == 0.0
+        #assert misc.trapzVariance(self.rectangleHeightUncertainty) == self.rectangleHeightUncertainty[0]*(len(self.rectangleHeightUncertainty)-1)
+    
+    def test_interpVariance(self):
+        """Testing for expected return value."""
+        #checked by hand
+        assert np.array_equal(stats.interpVariance([1,4,6], [1,2,3,4,5,6,7,8], [2,4,3,5,2,5,1,2]),[4.0, 25.0,25.0])
+        assert np.array_equal(stats.interpVariance([0,4,4.5,6,9], [1,2,3,4,5,6,7,8], [2,4,3,5,2,5,1,2]),[4.0, 25.0, 7.25,25.0,4.0])
\ No newline at end of file
diff --git a/fiducia/tests/test_visualization.py b/fiducia/tests/test_visualization.py
new file mode 100644
index 0000000..96856ce
--- /dev/null
+++ b/fiducia/tests/test_visualization.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  8 10:49:45 2019
+
+Tests for visualization.py
+
+@author: Pawel M. Kozlowski
+"""
+
+# python modules
+import numpy as np
+import pytest
+
+# custom modules
+import fiducia.visualization as vis
\ No newline at end of file