diff --git a/pycbc/waveform/waveform_modes.py b/pycbc/waveform/waveform_modes.py
index 3601a2034c1..ca1f70ebee0 100644
--- a/pycbc/waveform/waveform_modes.py
+++ b/pycbc/waveform/waveform_modes.py
@@ -232,6 +232,66 @@ def get_nrhybsur_modes(**params):
 get_nrsur_modes.__doc__ = _formatdocstr(get_nrsur_modes.__doc__)
 get_nrhybsur_modes.__doc__ = _formatdocstr(get_nrhybsur_modes.__doc__)
 
+def get_lalsimulation_modes(**params):
+    """Generates approximant waveform mode-by-mode.
+
+    All waveform parameters should be provided as keyword arguments.
+    Recognized parameters are listed below. Unrecognized arguments are ignored.
+
+    Parameters
+    ----------
+    template: object
+        An object that has attached properties. This can be used to substitute
+        for keyword arguments. A common example would be a row in an xml table.
+    approximant : str
+        The approximant to generate. Must be available in ``lalsimulation``.
+    {delta_t}
+    {mass1}
+    {mass2}
+    {spin1x}
+    {spin1y}
+    {spin1z}
+    {spin2x}
+    {spin2y}
+    {spin2z}
+    {f_lower}
+    {f_ref}
+    {distance}
+    {mode_array}
+    {ell_max}
+
+    Returns
+    -------
+    dict :
+        Dictionary of ``(l, m)`` -> ``(h_+, -h_x)`` ``TimeSeries``.
+    """
+    ell_max = 5
+    if 'ell_max' in params:
+        ell_max = params['ell_max']
+    laldict = _check_lal_pars(params)
+    ret = lalsimulation.SimInspiralChooseTDModes(
+        params['coa_phase'],
+        params['delta_t'],
+        params['mass1']*lal.MSUN_SI,
+        params['mass2']*lal.MSUN_SI,
+        params['spin1x'],
+        params['spin1y'],
+        params['spin1z'],
+        params['spin2x'],
+        params['spin2y'],
+        params['spin2z'],
+        params['f_lower'], params['f_ref'],
+        params['distance']*1e6*lal.PC_SI, laldict,
+        ell_max,
+        params['approx']
+    )
+    hlms = {}
+    while ret:
+        hlm = TimeSeries(ret.mode.data.data, delta_t=ret.mode.deltaT,
+                         epoch=ret.mode.epoch)
+        hlms[ret.l, ret.m] = (hlm.real(), hlm.imag())
+        ret = ret.next
+    return hlms
 
 def get_imrphenomxh_modes(**params):
     """Generates ``IMRPhenomXHM`` waveforms mode-by-mode. """
@@ -271,6 +331,8 @@ def get_imrphenomxh_modes(**params):
 
 _mode_waveform_td = {'NRSur7dq4': get_nrsur_modes,
                      'NRHybSur3dq8': get_nrhybsur_modes,
+                     'SEOBNRv4P': get_lalsimulation_modes,
+                     'SEOBNRv4PHM': get_lalsimulation_modes,
                      }
 _mode_waveform_fd = {'IMRPhenomXHM': get_imrphenomxh_modes,
                      }