Suspension thermal noise sub-budgets

gwinc/ifo/noises.py

@@ -461,20 +461,168 @@ class NewtonianInfrasound(nb.Noise):
 # suspension thermal
 #########################
 
-class SuspensionThermal(nb.Noise):
+class SuspensionThermalHorizTop(nb.Noise):
+    """Horizontal suspension thermal around the top mass
+
+    """
+    style = dict(
+        label='Horiz. Top',
+        color='xkcd:orangeish',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 0, 'horiz')
+        return abs(n) * 4
+
+
+class SuspensionThermalHorizAPM(nb.Noise):
+    """Horizontal suspension thermal around the upper intermediate mass
+
+    """
+    style = dict(
+        label='Horiz. APM',
+        color='xkcd:mustard',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 1, 'horiz')
+        return abs(n) * 4
+
+
+class SuspensionThermalHorizPUM(nb.Noise):
+    """Horizontal suspension thermal around the penultimate mass
+
+    """
+    style = dict(
+        label='Horiz. PUM',
+        color='xkcd:turquoise',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 2, 'horiz')
+        return abs(n) * 4
+
+
+class SuspensionThermalHorizTM(nb.Noise):
+    """Horizontal suspension thermal around the test
+
+    """
+    style = dict(
+        label='Horiz. Test mass',
+        color='xkcd:bright purple',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        precomp_suspension(self.freq, self.ifo)
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 3, 'horiz')
+        return abs(n) * 4
+
+
+class SuspensionThermalVertTop(nb.Noise):
+    """Vertical suspension thermal around the top mass
+
+    """
+    style = dict(
+        label='Vert. Top',
+        color='xkcd:orangeish',
+        linestyle='--',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 0, 'vert')
+        return abs(n) * 4
+
+
+class SuspensionThermalVertAPM(nb.Noise):
+    """Vertical suspension thermal around the upper intermediate mass
+
+    """
+    style = dict(
+        label='Vert. APM',
+        color='xkcd:mustard',
+        linestyle='--',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 1, 'vert')
+        return abs(n) * 4
+
+
+class SuspensionThermalVertPUM(nb.Noise):
+    """Vertical suspension thermal around the penultimate mass
+
+    """
+    style = dict(
+        label='Vert. PUM',
+        color='xkcd:turquoise',
+        linestyle='--',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 2, 'vert')
+        return abs(n) * 4
+
+
+class SuspensionThermalVertTM(nb.Noise):
+    """Vertical suspension thermal around the test
+
+    """
+    style = dict(
+        label='Vert. Test mass',
+        color='xkcd:bright purple',
+        linestyle='--',
+        alpha=0.7,
+    )
+
+    @nb.precomp(sustf=precomp_suspension)
+    def calc(self, sustf):
+        n = noise.suspensionthermal.susptherm_stage(
+            self.freq, self.ifo.Suspension, sustf, 3, 'vert')
+        return abs(n) * 4
+
+
+class SuspensionThermal(nb.Budget):
     """Suspension Thermal
 
     """
+
+    name = 'SuspensionThermal'
+
     style = dict(
         label='Suspension Thermal',
         color='#0d75f8',
     )
 
-    @nb.precomp(sustf=precomp_suspension)
-    def calc(self, sustf):
-        n = noise.suspensionthermal.suspension_thermal(
-            self.freq, self.ifo.Suspension, sustf)
-        return n * 4
+    noises = [
+        SuspensionThermalHorizTop,
+        SuspensionThermalHorizAPM,
+        SuspensionThermalHorizPUM,
+        SuspensionThermalHorizTM,
+        SuspensionThermalVertTop,
+        SuspensionThermalVertAPM,
+        SuspensionThermalVertPUM,
+    ]
 
 
 #########################

gwinc/noise/suspensionthermal.py

@@ -9,13 +9,46 @@ from ..const import kB
 from ..suspension import getJointParams
 
 
-def suspension_thermal(f, sus, sustf):
-    """Suspension thermal noise.
+def getJointTempSusceptibility(sus, sustf, stage_num, isUpper, direction):
+    """Get the product of the temperature and the susceptibility of a joint
 
-    :f: frequency array in Hz
-    :sus: gwinc suspension structure
-    :sustf: sus transfer function Struct
+    :sus: the suspension struct
+    :sustf: suspension transfer function struct
+    :stage_num: stage number with 0 at the top
+    :isUpper: whether to get the upper or lower joint
+    :direction: either 'horiz' or 'vert'
+    """
+    jointParams = getJointParams(sus, stage_num)
+    for (isLower, Temp, wireMat, bladeMat) in jointParams:
+        if isUpper == isLower:
+            continue
+
+        ind = 2*stage_num + isLower
+        # Compute the force on the TM along the beamline
+        if direction == 'horiz':
+            dxdF = sustf.hForce[ind]
+
+        elif direction == 'vert':
+            # vertical to beamline coupling angle
+            theta = sustf.VHCoupling.theta
+
+            # convert to beam line motion. theta is squared because
+            # we rotate by theta into the suspension basis, and by
+            # theta to rotate back to the beam line basis
+            dxdF = theta**2 * sustf.vForce[ind]
+
+        return Temp * abs(imag(dxdF))
 
+
+def susptherm_stage(f, sus, sustf, stage_num, direction):
+    """Suspension thermal noise of a single stage
+
+    :f: frequency array in Hz
+    :sus: gwinc suspension struct
+    :sustf: suspension transfer function struct
+    :stage_num: stage number with 0 at the top
+    :direction: either 'horiz' for horizontal or 'vert' for vertical
+    :comp: component of the loss 0) bulk, 1) surface, 2) TE
     :returns: displacement noise power spectrum at :f:
 
     :Temp: must either be set for each stage individually, or globally
@@ -26,29 +59,48 @@ def suspension_thermal(f, sus, sustf):
     pre-calculated and populated into the relevant `sus` fields.
 
     """
-    # suspension TFs
-    hForce = sustf.hForce
-    vForce = sustf.vForce
-    # and vertical to beamline coupling angle
-    theta = sustf.VHCoupling.theta
-
-    noise = np.zeros_like(f)
 
-    # Here the suspension stage list is reversed so that the last stage
-    # in the list is the test mass
-    stages = sus.Stage[::-1]
     w = 2*pi*f
-    dxdF = np.zeros_like(hForce)  # complex valued
-    for n, stage in enumerate(stages):
-        # add up the contribution from each joint
-        jointParams = getJointParams(sus, n)
-        for (isLower, Temp, wireMat, bladeMat) in jointParams:
-            # convert to beam line motion.  theta is squared because
-            # we rotate by theta into the suspension basis, and by
-            # theta to rotate back to the beam line basis
-            m = 2*n + isLower
-            dxdF[m, :] = hForce[m, :] + theta**2 * vForce[m, :]
-            noise += 4 * kB * Temp * abs(imag(dxdF[m, :])) / w
+    last_stage = len(sus.Stage) - 1
+
+    # get the noise from the lower joint
+    noise = getJointTempSusceptibility(sus, sustf, stage_num, False, direction)
+
+    # if this is the top mass, also get the noise from the upper joint
+    if stage_num == 0:
+        noise += getJointTempSusceptibility(sus, sustf, 0, True, direction)
+
+    # if this isn't the test mass, also get the noise from the upper joint of
+    # the mass below this one
+    if stage_num < last_stage:
+        noise += getJointTempSusceptibility(
+            sus, sustf, stage_num + 1, True, direction)
 
     # thermal noise (m^2/Hz) for one suspension
+    noise = 4 * kB * noise / w
+
     return noise
+
+
+def suspension_thermal(f, sus, sustf):
+    """Total suspension thermal noise of one suspension
+
+    :f: frequency array in Hz
+    :sus: gwinc suspension structure
+    :sustf: suspension transfer function struct
+    :returns: displacement noise power spectrum at :f:
+
+    :Temp: must either be set for each stage individually, or globally
+    in :sus.Temp:.  If both are specified, :sus.Temp: is interpreted as
+    the temperature of the upper joint of the top stage.
+
+    Assumes suspension transfer functions and V-H coupling have been
+    pre-calculated and populated into the relevant `sus` fields.
+
+    """
+    nstage = len(sus.Stage)
+    noise = np.zeros_like(f)
+    for stage_num in range(nstage):
+        noise += susptherm_stage(f, sus, sustf, stage_num, 'horiz')
+        noise += susptherm_stage(f, sus, sustf, stage_num, 'vert')
+    return abs(noise)

gwinc/suspension.py

@@ -553,8 +553,8 @@ def suspQuad(f, sus):
             kvi[n, 0, :] = imag(kv4)
             kvi[n, 1, :] = 0
             khr[n, :] = real(kh4)
-            khi[n, 0, :] = imag(kh4)
-            khi[n, 1, :] = 0
+            khi[n, 0, :] = imag(kh4) / 2
+            khi[n, 1, :] = imag(kh4) / 2
 
 
     ###############################################################