diff --git a/gwinc/gwinc.py b/gwinc/gwinc.py
index dab83d96976ed7c4807d103780b872759c4354cb..6cd0128a31a10960aef82c1fc82383fd9af7e44d 100644
--- a/gwinc/gwinc.py
+++ b/gwinc/gwinc.py
@@ -49,7 +49,7 @@ def noise_calc(ifo, f):
     noises['Coating Brownian']  = noise.coatingthermal.coatbrownian(f, ifo)
     noises['Substrate Thermo-Elastic']  = noise.substratethermal.subtherm(f, ifo)
     noises['Newtonian Gravity']  = noise.newtonian.gravg(f, ifo)
-    noises['Seismic'] = noise.seismic.seismic(f, ifo)
+    noises['Seismic'] = noise.seismic.seismic(f, ifo)[0]
     noises['Coating Thermo-Optic'] = noise.coatingthermal.thermooptic(f, ifo)
 
     # calc semiconductor noise sources
diff --git a/gwinc/noise/seismic.py b/gwinc/noise/seismic.py
index c49cd436e2000f7c255bdf384d0c414f89fc3e39..cb66bf493480dcef7ffb51ccd8af9900752e283e 100644
--- a/gwinc/noise/seismic.py
+++ b/gwinc/noise/seismic.py
@@ -5,28 +5,16 @@ from scipy.interpolate import interp1d
 
 
 def seismic(f, ifo):
-    """seismic noise psd at frequencies f for given ifo.
-      n = seismic(f,ifo)
-      [nh,nv] = seismic(f,ifo)
-      [n,nh,nv] = seismic(f,ifo)
-    
-    Modified to include realistic SEI + SUS models (Rana, 11/2005)"""
+    """Seismic noise.
 
-    # Interpolate the log10 onto the ifo frequency array
-    #   n = interp1(ifo.Seismic.darmseis_f, ...
-    #     log10(ifo.Seismic.darmseis_x), f, 'cubic', -30);
+    Return (noise, noise_vertical, noise_horizontal)
 
-    ##########################################################
-    # Suspension TFs
-    ##########################################################
+    """
     hTable = ifo.Suspension.hTable
     vTable = ifo.Suspension.vTable
 
-    # and vertical to beamline coupling angle
     theta = ifo.Suspension.VHCoupling.theta
 
-    ##########################################################
-
     # noise input, horizontal and vertical
     nx, np_ = seisBSC(f)
 
@@ -44,15 +32,15 @@ def seismic(f, ifo):
     nv *= 4 / ifo.Infrastructure.Length**2
     n  *= 4 / ifo.Infrastructure.Length**2
 
-    return n
+    return n, nh, nv
 
 
 def seisBSC(f):
-    """get rough ISI noise source spectra
-    nx - ISI translational DOFs
-    np - ISI rotational DOFs"""
+    """Rough ISI noise source spectra.
 
+    Returns (ISI translational DOFs, ISI rotational DOFs)
 
+    """
     # translational DOFs (from Rana's bsc_seismic.m)
     SEI_F = np.array([0.01, 0.03, 0.1, 0.2, 0.5, 1, 10, 30, 300])
     SEI_X = np.array([3e-6, 1e-6, 2e-7, 2e-7, 8e-10, 1e-11, 3e-13, 3e-14, 3e-14])