Merge branch 'seismic-sub-nb' into 'master'

Seismic noise sub-budgets

See merge request !110

gwinc/ifo/noises.py

@@ -14,7 +14,6 @@ from .. import suspension
 # helper functions
 ############################################################
 
-
 def mirror_struct(ifo, tm):
     """Create a "mirror" Struct for a LIGO core optic
 
@@ -126,7 +125,6 @@ def ifo_power(ifo, PRfixed=True):
 
 ##################################################
 
-
 def precomp_suspension(f, ifo):
     pc = Struct()
     pc.VHCoupling = Struct()
@@ -142,7 +140,6 @@ def precomp_suspension(f, ifo):
     return pc
 
 
-
 def precomp_quantum(f, ifo):
     pc = Struct()
     mirror_mass = mirror_struct(ifo, 'ETM').MirrorMass
@@ -345,35 +342,58 @@ class StandardQuantumLimit(nb.Noise):
         ETM = mirror_struct(self.ifo, 'ETM')
         return 8 * const.hbar / (ETM.MirrorMass * (2 * np.pi * self.freq) ** 2)
 
+
 #########################
 # seismic
 #########################
 
-class Seismic(nb.Noise):
-    """Seismic
+class SeismicHorizontal(nb.Noise):
+    """Horizontal seismic noise
 
     """
     style = dict(
-        label='Seismic',
-        color='#855700',
+        label='Horizontal',
+        color='xkcd:muted blue',
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        nt, nr = noise.seismic.platform_motion(self.freq, self.ifo)
+        n = noise.seismic.seismic_suspension_filtered(sustf, nt, 'horiz')
+        return n * 4
+
+
+class SeismicVertical(nb.Noise):
+    """Vertical seismic noise
+
+    """
+    style = dict(
+        label='Vertical',
+        color='xkcd:brick red',
     )
 
     @nb.precomp(sustf=precomp_suspension)
     def calc(self, sustf):
-        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, sustf, nt)
+        nt, nr = noise.seismic.platform_motion(self.freq, self.ifo)
+        n = noise.seismic.seismic_suspension_filtered(sustf, nt, 'vert')
         return n * 4
 
 
+class Seismic(nb.Budget):
+    """Seismic
+
+    """
+    style = dict(
+        label='Seismic',
+        color='#855700',
+    )
+
+    noises = [
+        SeismicHorizontal,
+        SeismicVertical,
+    ]
+
+
 #########################
 # Newtonian
 #########################
@@ -572,7 +592,6 @@ class SubstrateThermoElastic(nb.Noise):
 # residual gas
 #########################
 
-
 class ExcessGas(nb.Noise):
     """Excess Gas
 

gwinc/noise/seismic.py

@@ -6,32 +6,50 @@ import numpy as np
 from scipy.interpolate import PchipInterpolator as interp1d
 
 
-def seismic_suspension_fitered(sus, sustf, in_trans):
+def seismic_suspension_filtered(sustf, in_trans, direction):
     """Seismic displacement noise for single suspended test mass.
 
-    :sus: suspension Struct
     :sustf: sus transfer function Struct
     :in_trans: input translational displacement spectrum
+    :direction: 'horiz' for horizontal or 'vert' for vertical
 
-    :returns: tuple of displacement noise power spectrum at :f:, and
-    horizontal and vertical components.
+    :returns: tuple of displacement noise power spectrum at :f:
 
     """
-    hTable = sustf.hTable
-    vTable = sustf.vTable
+    if direction == 'horiz':
+        # horizontal noise total
+        n = (abs(sustf.hTable)**2) * in_trans**2
 
-    theta = sustf.VHCoupling.theta
+    elif direction == 'vert':
+        # vertical to horizontal coupling
+        theta = sustf.VHCoupling.theta
 
-    # horizontal noise total
-    nh = (abs(hTable)**2) * in_trans**2
+        # vertical noise total
+        n = (abs(theta * sustf.vTable)**2) * in_trans**2
 
-    # vertical noise total
-    nv = (abs(theta * vTable)**2) * in_trans**2
+    return n
 
-    # new total noise
-    n = nv + nh
 
-    return n, nh, nv
+def platform_motion(f, ifo):
+    """Compute the platform motion
+
+    :f: frequency array in Hz
+    :ifo: the IFO struct
+
+    :returns: tuple of displacement noise power spectrum at :f: for
+    translational and rotational DOFs.
+    """
+    if 'PlatformMotion' in ifo.Seismic:
+        if ifo.Seismic.PlatformMotion == 'BSC':
+            nt, nr = seismic_BSC_ISI(f)
+        elif ifo.Seismic.PlatformMotion == '6D':
+            nt, nr = seismic_BSC_ISI_6D(f)
+        else:
+            nt, nr = seismic_BSC_ISI(f)
+    else:
+        nt, nr = seismic_BSC_ISI(f)
+
+    return nt, nr
 
 
 def seismic_BSC_ISI(f):