diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py
index 74e4598678a855fc76083862964a4807cb5e7975..37a7c2f8e6ffbd8ce670939740a5facecfe23fa6 100644
--- a/gwinc/ifo/noises.py
+++ b/gwinc/ifo/noises.py
@@ -24,7 +24,18 @@ class Seismic(nb.Noise):
     )
 
     def calc(self):
-        return noise.seismic.seismic(self.freq, self.ifo)
+        if 'PlatformMotion' in self.ifo.Seismic:
+            if self.ifo.Seismic.PlatformMotion == 'BSC':
+                nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
+            elif self.ifo.Seismic.PlatformMotion == '6D':
+                nt, nr = noise.seismic.seismic_BSC_ISI_6D(self.freq)
+            else:
+                nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
+        else:
+            nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
+        n, nh, nv = noise.seismic.seismic_suspension_fitered(
+            self.ifo.Suspension, nt)
+        return n * 4 * self.ifo.gwinc.dhdl_sqr
 
 
 class Newtonian(nb.Noise):
diff --git a/gwinc/noise/newtonian.py b/gwinc/noise/newtonian.py
index 0fc6819254ec5bd4d83e5bd0c7d3f9e49e02a401..3ab420304d7c9eab5da6fd1c32ebf02d4ba7a7a4 100644
--- a/gwinc/noise/newtonian.py
+++ b/gwinc/noise/newtonian.py
@@ -4,7 +4,7 @@ import numpy as np
 import scipy.integrate as scint
 
 from numpy import pi, sqrt, exp
-from .seismic import seisNLNM
+from .seismic import seismic_ground_NLNM
 from .. import const
 
 
@@ -131,7 +131,7 @@ def gravg_pwave(f, ifo):
     kP = (2 * pi * f) / cP
 
     rho_ground = ifo.Seismic.Rho
-    psd_ground_pwave = (levelP * seisNLNM(f))**2
+    psd_ground_pwave = (levelP * seismic_ground_NLNM(f))**2
 
     tmheight = ifo.Seismic.TestMassHeight 
     xP = np.abs(kP * tmheight)
@@ -165,7 +165,7 @@ def gravg_swave(f, ifo):
     kS = (2 * pi * f) / cS
 
     rho_ground = ifo.Seismic.Rho
-    psd_ground_swave = (levelS * seisNLNM(f))**2
+    psd_ground_swave = (levelS * seismic_ground_NLNM(f))**2
 
     tmheight = ifo.Seismic.TestMassHeight 
     xS = np.abs(kS * tmheight)
diff --git a/gwinc/noise/seismic.py b/gwinc/noise/seismic.py
index fbee12da06100d4cecf54220cb6865b3d104079f..ac856f713b1f183a1a83e1f9082443779e61f913 100644
--- a/gwinc/noise/seismic.py
+++ b/gwinc/noise/seismic.py
@@ -1,90 +1,83 @@
+'''Functions to calculate seismic noise in suspended optics.
+
+'''
 from __future__ import division
 import numpy as np
-from numpy import log10
 from scipy.interpolate import PchipInterpolator as interp1d
 
 
-def seismic(f, ifo):
-    """Seismic noise.
-    """
-    return seismicAll(f, ifo)[0]
-
+def seismic_suspension_fitered(sus, in_trans):
+    """Seismic displacement noise for single suspended test mass.
 
-def seismicAll(f, ifo):
-    """Seismic noise.
+    :sus: gwinc suspension structure
+    :in_trans: input translational displacement spectrum
 
-    Return (noise, noise_vertical, noise_horizontal)
+    :returns: tuple of displacement noise power spectrum at :f:, and
+    horizontal and vertical components.
 
     """
-    hTable = ifo.Suspension.hTable
-    vTable = ifo.Suspension.vTable
-
-    theta = ifo.Suspension.VHCoupling.theta
-
-    # noise input, horizontal and vertical
-    if 'PlatformMotion' in ifo.Seismic:
-        if ifo.Seismic.PlatformMotion == 'BSC':
-            nt, nr = seisBSC(f)
-        elif ifo.Seismic.PlatformMotion == '6D':
-            nt, nr = seis6D(f)
-        else:
-            nt, nr = seisBSC(f)
-    else:
-        nt, nr = seisBSC(f)
+    hTable = sus.hTable
+    vTable = sus.vTable
+
+    theta = sus.VHCoupling.theta
 
     # horizontal noise total
-    nh = (abs(hTable)**2) * nt**2
+    nh = (abs(hTable)**2) * in_trans**2
 
     # vertical noise total
-    nv = (abs(theta * vTable)**2) * nt**2
+    nv = (abs(theta * vTable)**2) * in_trans**2
 
     # new total noise
     n = nv + nh
 
-    # Convert into Strain PSD (4 TMs)
-    nh *= 4 * ifo.gwinc.dhdl_sqr
-    nv *= 4 * ifo.gwinc.dhdl_sqr
-    n *= 4 * ifo.gwinc.dhdl_sqr
-
     return n, nh, nv
 
 
-def seisBSC(f):
-    """Rough ISI noise source spectra.
+def seismic_BSC_ISI(f):
+    """Rough seismic noise spectra on aLIGO BSC ISI table.
+
+    :f: frequency array in Hz
 
-    Returns ISI (translational, rotational) DOFs
+    :returns: tuple of displacement noise power spectrum at :f: for
+    translational and rotational DOFs.
 
     """
     SEI_F = np.array([0.01, 0.03, 0.1, 0.2, 0.5, 1, 10, 30, 300])
 
     # translational DOFs
     SEI_T = np.array([3e-6, 1e-6, 2e-7, 2e-7, 8e-10, 1e-11, 3e-13, 3e-14, 3e-14])
-    nt = 10**(interp1d(SEI_F, log10(SEI_T))(f))
+    nt = 10**(interp1d(SEI_F, np.log10(SEI_T))(f))
 
     # rotational DOFs
     SEI_R = np.array([1e-8, 3e-8, 2e-8, 1e-8, 4e-10, 1e-11, 3e-13, 3e-14, 3e-14])
-    nr = 10**(interp1d(SEI_F, log10(SEI_R))(f))
+    nr = 10**(interp1d(SEI_F, np.log10(SEI_R))(f))
 
     return nt, nr
 
 
-def seis6D(f):
-    """ISI noise source spectra with a 6D seismometer.
+def seismic_BSC_ISI_6D(f):
+    """Rough seismic noise spectra on aLIGO BSC ISI table with a 6D seismometer.
 
     This largely follows Mow-Lowry and Martynov, arXiv:1801.01468.
 
+    :f: frequency array in Hz
+
+    :returns: tuple of displacement noise power spectrum at :f: for
+    translational and rotational DOFs.
+
     """
+    # FIXME: merge this with above, using flag
 
     SEI_F = np.array([0.01, 0.03, 0.1, 0.2, 0.5, 1, 10, 100, 300])
 
     SEI_T_self = np.array([1e-7, 1e-9, 3e-11, 6e-12, 3e-13, 1e-13, 3e-14, 1e-14, 1e-14])
-    nt_self = 10**(interp1d(SEI_F, log10(SEI_T_self))(f))
-    nt_gnd = 10*seisNLNM(f)
+    nt_self = 10**(interp1d(SEI_F, np.log10(SEI_T_self))(f))
+    nt_gnd = 10*seismic_ground_NLNM(f)
     blend_t = np.abs(100/(1+1j*f/0.01)**4)
     nt = np.sqrt(nt_self**2 + (blend_t * nt_gnd)**2)
 
     SEI_R_self = np.array([2e-11, 5e-12, 1e-12, 6e-13, 3e-13, 2e-13, 6e-14, 2e-14, 2e-14])
-    nr_self = 10**(interp1d(SEI_F, log10(SEI_R_self))(f))
+    nr_self = 10**(interp1d(SEI_F, np.log10(SEI_R_self))(f))
     nr_gnd = np.abs(1e-7/(1+1j*f/0.001))
     blend_r = np.abs(100/(1+1j*f/0.01)**4)
     nr = np.sqrt(nr_self**2 + (blend_r * nr_gnd)**2)
@@ -92,35 +85,39 @@ def seis6D(f):
     return nt, nr
 
 
-def seisNLNM(f):
-    """The Peterson New Low-Noise Model.
+def seismic_ground_NLNM(f):
+    """The Peterson new generic ground motion low noise model.
+
+    :f: frequency array in Hz
 
-    Returns a displacement ASD.
+    :returns: displacement noise amplitude spectrum at :f:
 
     """
     Pl = np.array([
-       1.00e-02, 1.00e-01, 1.70e-01, 4.00e-01, 8.00e-01, 1.24e+00,
-       2.40e+00, 4.30e+00, 5.00e+00, 6.00e+00, 1.00e+01, 1.20e+01,
-       1.56e+01, 2.19e+01, 3.16e+01, 4.50e+01, 7.00e+01, 1.01e+02,
-       1.54e+02, 3.28e+02, 6.00e+02, 1.00e+04])
+        1.00e-02, 1.00e-01, 1.70e-01, 4.00e-01, 8.00e-01, 1.24e+00,
+        2.40e+00, 4.30e+00, 5.00e+00, 6.00e+00, 1.00e+01, 1.20e+01,
+        1.56e+01, 2.19e+01, 3.16e+01, 4.50e+01, 7.00e+01, 1.01e+02,
+        1.54e+02, 3.28e+02, 6.00e+02, 1.00e+04])
     Al = np.array([
-       -156.72, -162.36, -166.7 , -170.  , -166.4 , -168.6 , -159.98,
-       -141.1 ,  -71.36,  -97.26, -132.18, -205.27,  -37.65, -114.37,
-       -160.58, -187.5 , -216.47, -185.  , -168.34, -217.43, -258.28,
-       -346.88])
+        -156.72, -162.36, -166.7, -170.0, -166.4, -168.6, -159.98,
+        -141.1, -71.36, -97.26, -132.18, -205.27, -37.65, -114.37,
+        -160.58, -187.5, -216.47, -185.0, -168.34, -217.43, -258.28,
+        -346.88])
     Bl = np.array([
-          5.64,    5.64,    0.  ,   -8.3 ,   28.9 ,   52.48,   29.81,
-          0.  ,  -99.77,  -66.49,  -31.57,   36.16, -104.33,  -47.1 ,
-        -16.28,    0.  ,   15.7 ,    0.  ,   -7.61,   11.9 ,   26.6 ,
-         48.75])
+        5.64, 5.64, 0.0, -8.3, 28.9, 52.48, 29.81,
+        0.0, -99.77, -66.49, -31.57, 36.16, -104.33, -47.1,
+        -16.28, 0.0, 15.7, 0.0, -7.61, 11.9, 26.6,
+        48.75])
     nlnm = 10**(np.interp(1/f, Pl, Al+Bl*np.log10(Pl))/20) / (2 * np.pi * f)**2
     return nlnm
 
 
-def seisNHNM(f):
-    """The Peterson New High-Noise Model.
+def seismic_ground_NHNM(f):
+    """The Peterson new generic ground motion high noise model.
 
-    Returns a displacement ASD.
+    :f: frequency array in Hz
+
+    :returns: displacement noise amplitude spectrum at :f:
 
     """
 
@@ -128,16 +125,16 @@ def seisNHNM(f):
         1.00e-01, 2.20e-01, 3.20e-01, 8.00e-01, 3.80e+00,
         4.60e+00, 6.30e+00, 7.90e+00, 1.54e+01, 2.00e+01,
         3.54e+02,
-        ])  
+        ])
     Al = np.array([
         -108.73, -150.34, -122.31, -116.85, -108.48,
-         -74.66,    0.66,  -93.37,   73.54, -151.52,
+        -74.66, 0.66, -93.37, 73.54, -151.52,
         -206.66,
         ])
     Bl = np.array([
-         -17.23,  -80.50,  -23.87,   32.51,   18.08,
-         -32.95, -127.18,  -22.42, -162.98,   10.01,
-          31.63,
+        -17.23, -80.50, -23.87, 32.51, 18.08,
+        -32.95, -127.18, -22.42, -162.98, 10.01,
+        31.63,
         ])
     nhnm = 10**(np.interp(1/f, Pl, Al+Bl*np.log10(Pl))/20) / (2 * np.pi * f)**2
     return nhnm