Suspension thermal noise sub-budgets

Adds suspension thermal noise sub-budgets. This addresses 1/3 of #51. Suspension thermal noise can now be broken up into horizontal and vertical noise of each of the suspension stages. With one exception, which does not change the total noise, the calculation is the same as the current calculation but the contributions of each stage and direction are tracked separately.

In the current calculation the horizontal and vertical spring constants of each stage are broken up into an upper and lower joint and are stored separately. The total suspension thermal noise is then calculated by adding the contributions of each of these joints using the appropriate temperature of each joint. Instead of adding the contributions of each joint all at once, this MR makes it possible to track each stage with separate `Noise` classes and add them all to a sub-budget.

The noise for each stage is the sum of the noise from that stage's lower joint and the noise from the upper joint of the stage below. In the case of the top mass, the contribution from the upper joint is also included. In the current calculation, for the upper stages, the imaginary part of the horizontal spring constant is equally distributed between the upper and lower joints while it is only attributed to the upper joint for the test mass. The one change this MR makes to the calculation itself is to distribute this equally between the top and bottom joints as is done for the upper stages.

The imaginary part of the vertical spring constant is only attributed to the upper joint of each stage. The blade springs are the only contribution to the vertical spring constant for the upper stages while the lower stage also includes the continuum fibers or ribbons. It's unclear how to break this up into an upper and lower stage as is done in figure 14 of https://arxiv.org/abs/2001.11173, so this MR does not have vertical noise of the test mass. (The noise is included, it's just attributed entirely to the PUM.)

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
 
 
     ###############################################################