diff --git a/enterprise_extensions/blocks.py b/enterprise_extensions/blocks.py
index 8b38ae71..69283a2b 100644
--- a/enterprise_extensions/blocks.py
+++ b/enterprise_extensions/blocks.py
@@ -11,13 +11,18 @@
 from enterprise.signals import utils
 from enterprise.signals import gp_bases as gpb
 from enterprise.signals import gp_priors as gpp
+from enterprise.signals import deterministic_signals
 from enterprise import constants as const
 from . import gp_kernels as gpk
 from . import chromatic as chrom
 from . import model_orfs
 
+from enterprise_extensions import deterministic as ee_deterministic
+
 __all__ = ['white_noise_block',
            'red_noise_block',
+           'bwm_block',
+           'bwm_sglpsr_block'
            'dm_noise_block',
            'scattering_noise_block',
            'chromatic_noise_block',
@@ -101,7 +106,7 @@ def white_noise_block(vary=False, inc_ecorr=False, gp_ecorr=False,
 def red_noise_block(psd='powerlaw', prior='log-uniform', Tspan=None,
                     components=30, gamma_val=None, coefficients=False,
                     select=None, modes=None, wgts=None,
-                    break_flat=False, break_flat_fq=None):
+                    break_flat=False, break_flat_fq=None, logmin=None, logmax=None):
     """
     Returns red noise model:
         1. Red noise modeled as a power-law with 30 sampling frequencies
@@ -124,15 +129,21 @@ def red_noise_block(psd='powerlaw', prior='log-uniform', Tspan=None,
     if psd in ['powerlaw', 'powerlaw_genmodes', 'turnover',
                'tprocess', 'tprocess_adapt', 'infinitepower']:
         # parameters shared by PSD functions
-        if prior == 'uniform':
-            log10_A = parameter.LinearExp(-20, -11)
-        elif prior == 'log-uniform' and gamma_val is not None:
-            if np.abs(gamma_val - 4.33) < 0.1:
-                log10_A = parameter.Uniform(-20, -11)
+        if logmin is not None and logmax is not None:
+            if prior == 'uniform':
+                log10_A = parameter.LinearExp(logmin, logmax)
+            elif prior == 'log-uniform':
+                log10_A = parameter.Uniform(logmin, logmax)
+        else:
+            if prior == 'uniform':
+                log10_A = parameter.LinearExp(-20, -11)
+            elif prior == 'log-uniform' and gamma_val is not None:
+                if np.abs(gamma_val - 4.33) < 0.1:
+                    log10_A = parameter.Uniform(-20, -11)
+                else:
+                    log10_A = parameter.Uniform(-20, -11)
             else:
                 log10_A = parameter.Uniform(-20, -11)
-        else:
-            log10_A = parameter.Uniform(-20, -11)
 
         if gamma_val is not None:
             gamma = parameter.Constant(gamma_val)
@@ -214,6 +225,83 @@ def red_noise_block(psd='powerlaw', prior='log-uniform', Tspan=None,
 
     return rn
 
+def bwm_block(Tmin, Tmax, amp_prior='log-uniform',
+              skyloc=None, logmin=-18, logmax=-11,
+              name='bwm'):
+    """
+    Returns deterministic GW burst with memory model:
+        1. Burst event parameterized by time, sky location,
+        polarization angle, and amplitude
+    :param Tmin:
+        Min time to search, probably first TOA (MJD).
+    :param Tmax:
+        Max time to search, probably last TOA (MJD).
+    :param amp_prior:
+        Prior on log10_A. Default if "log-uniform". Use "uniform" for
+        upper limits.
+    :param skyloc:
+        Fixed sky location of BWM signal search as [cos(theta), phi].
+        Search over sky location if ``None`` given.
+    :param logmin:
+        log of minimum BWM amplitude for prior (log10)
+    :param logmax:
+        log of maximum BWM amplitude for prior (log10)
+    :param name:
+        Name of BWM signal.
+    """
+
+    # BWM parameters
+    amp_name = '{}_log10_A'.format(name)
+    if amp_prior == 'uniform':
+        log10_A_bwm = parameter.LinearExp(logmin, logmax)(amp_name)
+    elif amp_prior == 'log-uniform':
+        log10_A_bwm = parameter.Uniform(logmin, logmax)(amp_name)
+
+    pol_name = '{}_pol'.format(name)
+    pol = parameter.Uniform(0, np.pi)(pol_name)
+
+    t0_name = '{}_t0'.format(name)
+    t0 = parameter.Uniform(Tmin, Tmax)(t0_name)
+
+    costh_name = '{}_costheta'.format(name)
+    phi_name = '{}_phi'.format(name)
+    if skyloc is None:
+        costh = parameter.Uniform(-1, 1)(costh_name)
+        phi = parameter.Uniform(0, 2*np.pi)(phi_name)
+    else:
+        costh = parameter.Constant(skyloc[0])(costh_name)
+        phi = parameter.Constant(skyloc[1])(phi_name)
+
+
+    # BWM signal
+    bwm_wf = ee_deterministic.bwm_delay(log10_h=log10_A_bwm, t0=t0,
+                            cos_gwtheta=costh, gwphi=phi, gwpol=pol)
+    bwm = deterministic_signals.Deterministic(bwm_wf, name=name)
+
+    return bwm
+
+def bwm_sglpsr_block(Tmin, Tmax, amp_prior='log-uniform',
+               logmin=-17, logmax=-12, name='ramp', fixed_sign=None):
+
+    if fixed_sign is None:
+        sign = parameter.Uniform(-1, 1)("sign")
+    else:
+        sign = np.sign(fixed_sign)
+
+    amp_name = '{}_log10_A'.format(name)
+    if amp_prior == 'uniform':
+        log10_A_ramp = parameter.LinearExp(logmin, logmax)(amp_name)
+    elif amp_prior == 'log-uniform':
+        log10_A_ramp = parameter.Uniform(logmin, logmax)(amp_name)
+
+    t0_name = '{}_t0'.format(name)
+    t0 = parameter.Uniform(Tmin, Tmax)(t0_name)
+
+
+    ramp_wf = ee_deterministic.bwm_sglpsr_delay(log10_A=log10_A_ramp, t0=t0, sign = sign)
+    ramp = deterministic_signals.Deterministic(ramp_wf, name=name)
+
+    return ramp
 
 def dm_noise_block(gp_kernel='diag', psd='powerlaw', nondiag_kernel='periodic',
                    prior='log-uniform', dt=15, df=200,
@@ -493,6 +581,7 @@ def chromatic_noise_block(gp_kernel='nondiag', psd='powerlaw',
 def common_red_noise_block(psd='powerlaw', prior='log-uniform',
                            Tspan=None, components=30,
                            log10_A_val = None, gamma_val=None, delta_val=None,
+                           logmin = None, logmax = None,
                            orf=None, orf_ifreq=0, leg_lmax=5,
                            name='gw', coefficients=False,
                            pshift=False, pseed=None):
@@ -519,7 +608,11 @@ def common_red_noise_block(psd='powerlaw', prior='log-uniform',
         models. By default spectral index is varied of range [0,7]
     :param delta_val:
         Value of spectral index for high frequencies in broken power-law
-        and turnover models. By default spectral index is varied in range [0,7].
+        and turnover models. By default spectral index is varied in range [0,7].\
+    :param logmin:
+        Specify the lower bound of the prior on the amplitude.
+    :param logmax:
+        Specify the upper bound of the prior on the amplitude
     :param orf:
         String representing which overlap reduction function to use.
         By default we do not use any spatial correlations. Permitted
@@ -539,6 +632,7 @@ def common_red_noise_block(psd='powerlaw', prior='log-uniform',
 
     """
 
+
     orfs = {'crn': None, 'hd': model_orfs.hd_orf(),
             'gw_monopole': model_orfs.gw_monopole_orf(),
             'gw_dipole': model_orfs.gw_dipole_orf(),
@@ -565,6 +659,18 @@ def common_red_noise_block(psd='powerlaw', prior='log-uniform',
         amp_name = '{}_log10_A'.format(name)
         if log10_A_val is not None:
             log10_Agw = parameter.Constant(log10_A_val)(amp_name)
+
+        if logmin is not None and logmax is not None:
+            if prior == 'uniform':
+                log10_Agw = parameter.LinearExp(logmin, logmax)(amp_name)
+            elif prior == 'log-uniform' and gamma_val is not None:
+                if np.abs(gamma_val - 4.33) < 0.1:
+                    log10_Agw = parameter.Uniform(logmin, logmax)(amp_name)
+                else:
+                    log10_Agw = parameter.Uniform(logmin, logmax)(amp_name)
+            else:
+                log10_Agw = parameter.Uniform(logmin, logmax)(amp_name)
+
         else:
             if prior == 'uniform':
                 log10_Agw = parameter.LinearExp(-18, -11)(amp_name)
diff --git a/enterprise_extensions/deterministic.py b/enterprise_extensions/deterministic.py
index b84af05a..5aacba7b 100644
--- a/enterprise_extensions/deterministic.py
+++ b/enterprise_extensions/deterministic.py
@@ -10,60 +10,6 @@
 from enterprise import constants as const
 
 
-def bwm_block(Tmin, Tmax, amp_prior='log-uniform',
-              skyloc=None, logmin=-18, logmax=-11,
-              name='bwm'):
-    """
-    Returns deterministic GW burst with memory model:
-        1. Burst event parameterized by time, sky location,
-        polarization angle, and amplitude
-    :param Tmin:
-        Min time to search, probably first TOA (MJD).
-    :param Tmax:
-        Max time to search, probably last TOA (MJD).
-    :param amp_prior:
-        Prior on log10_A. Default if "log-uniform". Use "uniform" for
-        upper limits.
-    :param skyloc:
-        Fixed sky location of BWM signal search as [cos(theta), phi].
-        Search over sky location if ``None`` given.
-    :param logmin:
-        log of minimum BWM amplitude for prior (log10)
-    :param logmax:
-        log of maximum BWM amplitude for prior (log10)
-    :param name:
-        Name of BWM signal.
-    """
-
-    # BWM parameters
-    amp_name = '{}_log10_A'.format(name)
-    if amp_prior == 'uniform':
-        log10_A_bwm = parameter.LinearExp(logmin, logmax)(amp_name)
-    elif amp_prior == 'log-uniform':
-        log10_A_bwm = parameter.Uniform(logmin, logmax)(amp_name)
-
-    pol_name = '{}_pol'.format(name)
-    pol = parameter.Uniform(0, np.pi)(pol_name)
-
-    t0_name = '{}_t0'.format(name)
-    t0 = parameter.Uniform(Tmin, Tmax)(t0_name)
-
-    costh_name = '{}_costheta'.format(name)
-    phi_name = '{}_phi'.format(name)
-    if skyloc is None:
-        costh = parameter.Uniform(-1, 1)(costh_name)
-        phi = parameter.Uniform(0, 2*np.pi)(phi_name)
-    else:
-        costh = parameter.Constant(skyloc[0])(costh_name)
-        phi = parameter.Constant(skyloc[1])(phi_name)
-
-
-    # BWM signal
-    bwm_wf = utils.bwm_delay(log10_h=log10_A_bwm, t0=t0,
-                            cos_gwtheta=costh, gwphi=phi, gwpol=pol)
-    bwm = deterministic_signals.Deterministic(bwm_wf, name=name)
-
-    return bwm
 
 
 def fdm_block(Tmin, Tmax, amp_prior='log-uniform', name='fdm',
@@ -461,6 +407,74 @@ def cw_delay(toas, pos, pdist,
 
     return res
 
+@signal_base.function
+def bwm_delay(toas, pos, log10_h=-14.0, cos_gwtheta=0.0, gwphi=0.0, gwpol=0.0, t0=55000,
+                antenna_pattern_fn=None):
+    """
+    Function that calculates the earth-term gravitational-wave
+    burst-with-memory signal, as described in:
+    Seto et al, van haasteren and Levin, phsirkov et al, Cordes and Jenet.
+    This version uses the F+/Fx polarization modes, as verified with the
+    Continuous Wave and Anisotropy papers.
+
+    :param toas: Time-of-arrival measurements [s]
+    :param pos: Unit vector from Earth to pulsar
+    :param log10_h: log10 of GW strain
+    :param cos_gwtheta: Cosine of GW polar angle
+    :param gwphi: GW azimuthal polar angle [rad]
+    :param gwpol: GW polarization angle
+    :param t0: Burst central time [day]
+    :param antenna_pattern_fn:
+        User defined function that takes `pos`, `gwtheta`, `gwphi` as
+        arguments and returns (fplus, fcross)
+
+    :return: the waveform as induced timing residuals (seconds)
+    """
+
+    # convert
+    h = 10 ** log10_h
+    gwtheta = np.arccos(cos_gwtheta)
+    t0 *= const.day
+
+    # antenna patterns
+    if antenna_pattern_fn is None:
+        apc = create_gw_antenna_pattern(pos, gwtheta, gwphi)
+    else:
+        apc = antenna_pattern_fn(pos, gwtheta, gwphi)
+
+    # grab fplus, fcross
+    fp, fc = apc[0], apc[1]
+
+    # combined polarization
+    pol = np.cos(2 * gwpol) * fp + np.sin(2 * gwpol) * fc
+
+    # Return the time-series for the pulsar
+    return pol * h * np.heaviside(toas - t0, 0.5) * (toas - t0)
+
+@signal_base.function
+def bwm_sglpsr_delay(toas, sign, log10_A=-15, t0=55000):
+
+    """
+    Function that calculates the earth-term gravitational-wave
+    burst-with-memory signal for an optimally oriented source in a single pulsar
+
+    :param toas: Time-of-arrival measurements [s]
+    :param log10_A: log10 of the amplitude of the ramp (delta_f/f)
+    :param t0: Burst central time [day]
+
+    :return: the waveform as induced timing residuals (seconds)
+    """
+
+
+    A = 10 ** log10_A
+    t0 *= const.day
+    # Return the time-series for the pulsar
+    heaviside = lambda x: 0.5 * (np.sign(x) + 1)
+
+    #return 0 #Fix the return to 0 in order to test what the heck is wrong with red noise detection in bwm
+    return A * np.sign(sign) * heaviside(toas - t0) * (toas - t0)
+
+
 
 @signal_base.function
 def compute_eccentric_residuals(toas, theta, phi, cos_gwtheta, gwphi,
diff --git a/enterprise_extensions/models.py b/enterprise_extensions/models.py
index 33179d12..860b575f 100644
--- a/enterprise_extensions/models.py
+++ b/enterprise_extensions/models.py
@@ -19,10 +19,14 @@
 from enterprise_extensions.blocks import (white_noise_block, red_noise_block,
                                           dm_noise_block,
                                           chromatic_noise_block,
-                                          common_red_noise_block)
+                                          common_red_noise_block,
+                                          bwm_block,
+                                          bwm_sglpsr_block)
 from enterprise_extensions.chromatic.solar_wind import solar_wind_block
 from enterprise_extensions import chromatic as chrom
 from enterprise_extensions import dropout as do
+from enterprise.signals.signal_base import LogLikelihood
+#from enterprise.signals.signal_base import LookupLikelihood
 
 
 def model_singlepsr_noise(psr, tm_var=False, tm_linear=False,
@@ -572,9 +576,9 @@ def model_2a(psrs, psd='powerlaw', noisedict=None, components=30,
 def model_general(psrs, tm_var=False, tm_linear=False, tmparam_list=None,
                   tm_svd=False, tm_norm=True, noisedict=None, white_vary=False,
                   Tspan=None, modes=None, wgts=None, logfreq=False, nmodes_log=10,
-                  common_psd='powerlaw', common_components=30, 
+                  common_psd='powerlaw', common_components=30,
                   log10_A_common=None, gamma_common=None,
-                  orf='crn', orf_names=None, orf_ifreq=0, leg_lmax=5, 
+                  orf='crn', orf_names=None, orf_ifreq=0, leg_lmax=5,
                   upper_limit_common=None, upper_limit=False,
                   red_var=True, red_psd='powerlaw', red_components=30, upper_limit_red=None,
                   red_select=None, red_breakflat=False, red_breakflat_fq=None,
@@ -619,7 +623,7 @@ def model_general(psrs, tm_var=False, tm_linear=False, tmparam_list=None,
         [default = 'powerlaw']
     :param common_components: number of frequencies starting at 1/T for common process.
         [default = 30]
-    :param log10_A_common: value of fixed log10_A_common parameter for 
+    :param log10_A_common: value of fixed log10_A_common parameter for
         fixed amplitude analyses.
         [default = None]
     :param gamma_common: fixed common red process spectral index value. By default we
@@ -627,18 +631,18 @@ def model_general(psrs, tm_var=False, tm_linear=False, tmparam_list=None,
         [default = None]
     :param orf: comma de-limited string of multiple common processes with different orfs.
         [default = crn]
-    :param orf_names: comma de-limited string of process names for different orfs. Manual 
+    :param orf_names: comma de-limited string of process names for different orfs. Manual
         control of these names is useful for embedding model_general within a hypermodel
-        analysis for a process with and without hd correlations where we want to avoid 
+        analysis for a process with and without hd correlations where we want to avoid
         parameter duplication.
         [default = None]
     :param orf_ifreq:
-        Frequency bin at which to start the Hellings & Downs function with 
+        Frequency bin at which to start the Hellings & Downs function with
         numbering beginning at 0. Currently only works with freq_hd orf.
         [default = 0]
     :param leg_lmax:
-        Maximum multipole of a Legendre polynomial series representation 
-        of the overlap reduction function. 
+        Maximum multipole of a Legendre polynomial series representation
+        of the overlap reduction function.
         [default = 5]
     :param upper_limit_common: perform upper limit on common red noise amplitude. Note
         that when perfoming upper limits it is recommended that the spectral index also
@@ -791,12 +795,12 @@ def model_general(psrs, tm_var=False, tm_linear=False, tmparam_list=None,
         else:
             log10_A_val = None
         crn.append(common_red_noise_block(psd=common_psd, prior=amp_prior_common, Tspan=Tspan,
-                                          components=common_components, 
+                                          components=common_components,
                                           log10_A_val=log10_A_val, gamma_val=gamma_common,
                                           delta_val=None, orf=elem, name='gw_{}'.format(elem_name),
                                           orf_ifreq=orf_ifreq, leg_lmax=leg_lmax,
                                           coefficients=coefficients, pshift=pshift, pseed=None))
-                                          # orf_ifreq only affects freq_hd model. 
+                                          # orf_ifreq only affects freq_hd model.
                                           # leg_lmax only affects (zero_diag_)legendre_orf model.
     crn = functools.reduce((lambda x,y:x+y), crn)
     s += crn
@@ -2093,12 +2097,11 @@ def model_chromatic(psrs, psd='powerlaw', noisedict=None, white_vary=False,
     return pta
 
 
-def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
-              Tmin_bwm=None, Tmax_bwm=None, skyloc=None,
-              red_psd='powerlaw', components=30,
+def model_bwm(psrs, likelihood=LogLikelihood,lookupdir=None, noisedict=None, tm_svd=False,
+              Tmin_bwm=None, Tmax_bwm=None, skyloc=None, logmin=None, logmax=None,
+              burst_logmin=-17, burst_logmax=-12, red_psd='powerlaw', components=30,
               dm_var=False, dm_psd='powerlaw', dm_annual=False,
-              upper_limit=False, bayesephem=False, is_wideband=False,
-              use_dmdata=False):
+              upper_limit=False, bayesephem=False, wideband=False):
     """
     Reads in list of enterprise Pulsar instance and returns a PTA
     instantiated with BWM model:
@@ -2119,8 +2122,6 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
     :param noisedict:
         Dictionary of pulsar noise properties for fixed white noise.
         Can provide manually, or the code will attempt to find it.
-    :param white_vary:
-        boolean for varying white noise or keeping fixed.
     :param tm_svd:
         boolean for svd-stabilised timing model design matrix
     :param Tmin_bwm:
@@ -2130,6 +2131,14 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
     :param skyloc:
         Fixed sky location of BWM signal search as [cos(theta), phi].
         Search over sky location if ``None`` given.
+    :param logmin:
+        Lower bound on log10_A of the red noise process in each pulsar`
+    :param logmax:
+        Upper bound on log10_A of the red noise process in each pulsar
+    :param burst_logmin:
+        Lower bound on the log10_A of the burst amplitude in each pulsar
+    :param burst_logmax:
+        Upper boudn on the log10_A of the burst amplitude in each pulsar
     :param red_psd:
         PSD to use for per pulsar red noise. Available options
         are ['powerlaw', 'turnover', tprocess, 'spectrum'].
@@ -2147,11 +2156,6 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
         set to False for a 'detection' run.
     :param bayesephem:
         Include BayesEphem model.
-    :param is_wideband:
-        Whether input TOAs are wideband TOAs; will exclude ecorr from the white
-        noise model.
-    :param use_dmdata: whether to use DM data (WidebandTimingModel) if
-        is_wideband.
     :return: instantiated enterprise.PTA object
     """
 
@@ -2168,7 +2172,7 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
         Tmax_bwm = tmax/const.day
 
     # red noise
-    s = red_noise_block(prior=amp_prior, psd=red_psd, Tspan=Tspan, components=components)
+    s = red_noise_block(prior=amp_prior, psd=red_psd, Tspan=Tspan, components=components, logmin=logmin, logmax=logmax)
 
     # DM variations
     if dm_var:
@@ -2181,7 +2185,7 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
         dmexp = chrom.dm_exponential_dip(tmin=54500, tmax=54900)
 
     # GW BWM signal block
-    s += deterministic.bwm_block(Tmin_bwm, Tmax_bwm,
+    s += bwm_block(Tmin_bwm, Tmax_bwm, logmin=burst_logmin, logmax=burst_logmax,
                                          amp_prior=amp_prior,
                                          skyloc=skyloc, name='bwm')
 
@@ -2190,35 +2194,18 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
         s += deterministic_signals.PhysicalEphemerisSignal(use_epoch_toas=True)
 
     # timing model
-    if (is_wideband and use_dmdata):
-        dmjump = parameter.Constant()
-        if white_vary:
-            dmefac = parameter.Uniform(pmin=0.1, pmax=10.0)
-            log10_dmequad = parameter.Uniform(pmin=-7.0, pmax=0.0)
-            #dmjump = parameter.Uniform(pmin=-0.005, pmax=0.005)
-        else:
-            dmefac = parameter.Constant()
-            log10_dmequad = parameter.Constant()
-            #dmjump = parameter.Constant()
-        s += gp_signals.WidebandTimingModel(dmefac=dmefac,
-                log10_dmequad=log10_dmequad, dmjump=dmjump,
-                dmefac_selection=selections.Selection(selections.by_backend),
-                log10_dmequad_selection=selections.Selection(
-                    selections.by_backend),
-                dmjump_selection=selections.Selection(selections.by_frontend))
-    else:
-        s += gp_signals.TimingModel(use_svd=tm_svd)
+    s += gp_signals.TimingModel(use_svd=tm_svd)
 
     # adding white-noise, and acting on psr objects
     models = []
     for p in psrs:
-        if 'NANOGrav' in p.flags['pta'] and not is_wideband:
-            s2 = s + white_noise_block(vary=white_vary, inc_ecorr=True)
+        if 'NANOGrav' in p.flags['pta'] and not wideband:
+            s2 = s + white_noise_block(vary=False, inc_ecorr=True)
             if dm_var and 'J1713+0747' == p.name:
                 s2 += dmexp
             models.append(s2(p))
         else:
-            s3 = s + white_noise_block(vary=white_vary, inc_ecorr=False)
+            s3 = s + white_noise_block(vary=False, inc_ecorr=False)
             if dm_var and 'J1713+0747' == p.name:
                 s3 += dmexp
             models.append(s3(p))
@@ -2227,12 +2214,129 @@ def model_bwm(psrs, noisedict=None, white_vary=False, tm_svd=False,
     pta = signal_base.PTA(models)
 
     # set white noise parameters
-    if not white_vary or (is_wideband and use_dmdata):
-        if noisedict is None:
-            print('No noise dictionary provided!...')
-        else:
-            noisedict = noisedict
-            pta.set_default_params(noisedict)
+    if noisedict is None:
+        print('No noise dictionary provided!...')
+    else:
+        noisedict = noisedict
+        pta.set_default_params(noisedict)
+
+    return pta
+
+def model_bwm_sglpsr (psr, likelihood=LogLikelihood, lookupdir=None, noisedict=None, tm_svd=False,
+              Tmin_bwm=None, Tmax_bwm=None,  burst_logmin=-17, burst_logmax=-12, fixed_sign=None,
+              red_psd='powerlaw', logmin=None, logmax=None, components=30,
+              dm_var=False, dm_psd='powerlaw', dm_annual=False,
+              upper_limit=False, bayesephem=False, wideband=False):
+    """
+    Burst-With-Memory model for single pulsar runs
+    Because all of the geometric parameters (pulsar_position, source_position, gw_pol) are all degenerate with each other in a single pulsar BWM search,
+    this model can only search over burst epoch and residual-space ramp amplitude (t0, ramp_amplitude)
+
+    Reads in list of enterprise Pulsar instance and returns a PTA
+    instantiated with single-pulsar BWM model (called a ramp):
+
+    per pulsar:
+        1. fixed EFAC per backend/receiver system
+        2. fixed EQUAD per backend/receiver system
+        3. fixed ECORR per backend/receiver system (if NG channelized)
+        4. Red noise modeled by a specified psd
+        5. Linear timing model.
+        6. Optional DM-variation modeling
+        7. Deterministic GW burst with memory signal for this pulsar
+
+    :param psr:
+        enterprise.Pulsar objects for PTA. This model is only for one pulsar at a time.
+    :param likelihood:
+        The likelihood function to use. The options are [enterprise.signals.signal_base.LogLikelihood, enterprise.signals.signal_base.LookupLikelihood]
+    :param noisedict:
+        Dictionary of pulsar noise properties for fixed white noise.
+        Can provide manually, or the code will attempt to find it.
+    :param tm_svd:
+        boolean for svd-stabilised timing model design matrix
+    :param Tmin_bwm:
+        Min time to search for BWM (MJD). If omitted, uses first TOA.
+    :param Tmax_bwm:
+        Max time to search for BWM (MJD). If omitted, uses last TOA.
+    :param red_psd:
+        PSD to use for per pulsar red noise. Available options
+        are ['powerlaw', 'turnover', tprocess, 'spectrum'].
+    :param components:
+        number of modes in Fourier domain processes (red noise, DM
+        variations, etc)
+    :param dm_var:
+        include gaussian process DM variations
+    :param dm_psd:
+        power-spectral density for gp DM variations
+    :param dm_annual:
+        include a yearly period DM variation
+    :param upper_limit:
+        Perform upper limit on BWM amplitude. By default this is
+        set to False for a 'detection' run.
+    :param bayesephem:
+        Include BayesEphem model.
+    :return: instantiated enterprise.PTA object
+
+
+    """
+    amp_prior = 'uniform' if upper_limit else 'log-uniform'
+
+    # find the maximum time span to set frequency sampling
+    tmin = psr.toas.min()
+    tmax = psr.toas.max()
+    Tspan = tmax - tmin
+
+    if Tmin_bwm is None:
+        Tmin_bwm = tmin/const.day
+    if Tmax_bwm is None:
+        Tmax_bwm = tmax/const.day
+
+    # red noise
+    s = red_noise_block(prior=amp_prior, psd=red_psd, Tspan=Tspan, components=components, logmin=logmin, logmax=logmax)
+
+    # DM variations
+    if dm_var:
+        s += dm_noise_block(psd=dm_psd, prior=amp_prior, components=components,
+                            gamma_val=None)
+        if dm_annual:
+            s += chrom.dm_annual_signal()
+
+        # DM exponential dip for J1713's DM event
+        dmexp = chrom.dm_exponential_dip(tmin=54500, tmax=54900)
+
+    # GW BWM signal block
+    s += bwm_sglpsr_block(Tmin_bwm, Tmax_bwm,amp_prior=amp_prior, name='ramp', logmin=burst_logmin, logmax=burst_logmax, fixed_sign=fixed_sign)
+
+    # ephemeris model
+    if bayesephem:
+        s += deterministic_signals.PhysicalEphemerisSignal(use_epoch_toas=True)
+
+    # timing model
+    s += gp_signals.TimingModel(use_svd=tm_svd)
+
+    # adding white-noise, and acting on psr objects
+    models = []
+
+    if 'NANOGrav' in psr.flags['pta'] and not wideband:
+        s2 = s + white_noise_block(vary=False, inc_ecorr=True)
+        if dm_var and 'J1713+0747' == p.name:
+            s2 += dmexp
+        models.append(s2(psr))
+    else:
+        s3 = s + white_noise_block(vary=False, inc_ecorr=False)
+        if dm_var and 'J1713+0747' == p.name:
+            s3 += dmexp
+        models.append(s3(psr))
+
+    # set up PTA
+    #TODO: decide on a way to handle likelihood
+    pta = signal_base.PTA(models)
+
+    # set white noise parameters
+    if noisedict is None:
+        print('No noise dictionary provided!...')
+    else:
+        noisedict = noisedict
+        pta.set_default_params(noisedict)
 
     return pta
 
diff --git a/tests/test_models.py b/tests/test_models.py
index 0c610cd7..035d302d 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -299,6 +299,37 @@ def test_model1(dmx_psrs,caplog):
     m1=models.model_1(dmx_psrs,noisedict=noise_dict)
     assert hasattr(m1,'get_lnlikelihood')
 
+def test_modelbwmsglpsr(nodmx_psrs, caplog):
+    tmax = max([p.toas.max() for p in nodmx_psrs])
+    tmin = min([p.toas.min() for p in nodmx_psrs])
+    Tspan = tmax - tmin
+
+    nodmx_psr=nodmx_psrs[0]
+
+    m=models.model_bwm_sglpsr(nodmx_psr) # should I be testing the Log and Lookup Likelihoods?
+                                    # If this test belongs in enterprise/tests instead, do
+                                    # I need to include the lookup table in tests/data?
+    assert hasattr(m, 'get_lnlikelihood')
+    assert "ramp_log10_A" in m.param_names
+    assert "ramp_t0" in m.param_names
+
+def test_modelbwm(nodmx_psrs, caplog):
+    tmax = max([p.toas.max() for p in nodmx_psrs])
+    tmin = min([p.toas.min() for p in nodmx_psrs])
+    Tspan = tmax - tmin
+
+
+    m=models.model_bwm(nodmx_psrs) # should I be testing the Log and Lookup Likelihoods?
+                                    # If this test belongs in enterprise/tests instead, do
+                                    # I need to include the lookup table in tests/data?
+    assert hasattr(m, 'get_lnlikelihood')
+    assert "bwm_log10_A" in m.param_names
+    assert "bwm_t0" in m.param_names
+    assert "bwm_phi" in m.param_names
+    assert "bwm_pol" in m.param_names
+    assert "bwm_costheta" in m.param_names
+
+
 @pytest.mark.filterwarnings('ignore::DeprecationWarning')
 def test_model1(dmx_psrs,caplog):
     # caplog.set_level(logging.CRITICAL)