add suspension thermal sub-budgets

Noise from each stage is tracked separately, but the total calculation is unchanged.

gwinc/ifo/CE1/__init__.py

@@ -21,6 +21,29 @@ class Newtonian(nb.Budget):
     ]
 
 
+class SuspensionThermal(nb.Budget):
+    """Suspension Thermal
+
+    """
+
+    name = 'SuspensionThermal'
+
+    style = dict(
+        label='Suspension Thermal',
+        color='#0d75f8',
+    )
+
+    noises = [
+        SuspensionThermalHorizTop,
+        SuspensionThermalHorizAPM,
+        SuspensionThermalHorizPUM,
+        SuspensionThermalHorizTM,
+        SuspensionThermalVertTop,
+        SuspensionThermalVertAPM,
+        SuspensionThermalVertPUM,
+    ]
+
+
 class Coating(nb.Budget):
     """Coating Thermal
 

gwinc/ifo/CE2/__init__.py

@@ -21,6 +21,29 @@ class Newtonian(nb.Budget):
     ]
 
 
+class SuspensionThermal(nb.Budget):
+    """Suspension Thermal
+
+    """
+
+    name = 'SuspensionThermal'
+
+    style = dict(
+        label='Suspension Thermal',
+        color='#0d75f8',
+    )
+
+    noises = [
+        SuspensionThermalHorizTop,
+        SuspensionThermalHorizAPM,
+        SuspensionThermalHorizPUM,
+        SuspensionThermalHorizTM,
+        SuspensionThermalVertTop,
+        SuspensionThermalVertAPM,
+        SuspensionThermalVertPUM,
+    ]
+
+
 class Coating(nb.Budget):
     """Coating Thermal
 

gwinc/ifo/noises.py

@@ -9,8 +9,10 @@ from .. import nb
 from .. import noise
 from .. import suspension
 
-##################################################
 
+############################################################
+# helper functions
+############################################################
 
 def coating_thickness(materials, optic):
     if 'CoatLayerOpticalThickness' in optic:
@@ -111,32 +113,9 @@ def ifo_power(ifo, PRfixed=True):
     return pbs, parm, finesse, prfactor, Tpr
 
 
-def dhdl(f, armlen):
-    """Strain to length conversion for noise power spetra
-
-    This takes into account the GW wavelength and is only important
-    when this is comparable to the detector arm length.
-
-    From R. Schilling, CQG 14 (1997) 1513-1519, equation 5,
-    with n = 1, nu = 0.05, ignoring overall phase and cos(nu)^2.
-    A small value of nu is used instead of zero to avoid infinities.
-
-    Returns the square of the dh/dL function, and the same divided by
-    the arm length squared.
-
-    """
-    c = const.c
-    nu_small = 15*pi/180
-    omega_arm = pi * f * armlen / c
-    omega_arm_f = (1 - sin(nu_small)) * pi * f * armlen / c
-    omega_arm_b = (1 + sin(nu_small)) * pi * f * armlen / c
-    sinc_sqr = 4 / abs(sin(omega_arm_f) * exp(-1j * omega_arm) / omega_arm_f
-                       + sin(omega_arm_b) * exp(1j * omega_arm) / omega_arm_b)**2
-    dhdl_sqr = sinc_sqr / armlen**2
-    return dhdl_sqr, sinc_sqr
-
-##################################################
-
+############################################################
+# precomp functions
+############################################################
 
 def precomp_suspension(f, ifo):
     pc = Struct()
@@ -179,14 +158,49 @@ def precomp_cavity(f, ifo):
     pc.wBeam_ETM = wBeam_ETM
     return pc
 
-##################################################
+
+############################################################
+# calibration
+############################################################
+
+def dhdl(f, armlen):
+    """Strain to length conversion for noise power spetra
+
+    This takes into account the GW wavelength and is only important
+    when this is comparable to the detector arm length.
+
+    From R. Schilling, CQG 14 (1997) 1513-1519, equation 5,
+    with n = 1, nu = 0.05, ignoring overall phase and cos(nu)^2.
+    A small value of nu is used instead of zero to avoid infinities.
+
+    Returns the square of the dh/dL function, and the same divided by
+    the arm length squared.
+
+    """
+    c = const.c
+    nu_small = 15*pi/180
+    omega_arm = pi * f * armlen / c
+    omega_arm_f = (1 - sin(nu_small)) * pi * f * armlen / c
+    omega_arm_b = (1 + sin(nu_small)) * pi * f * armlen / c
+    sinc_sqr = 4 / abs(sin(omega_arm_f) * exp(-1j * omega_arm) / omega_arm_f
+                       + sin(omega_arm_b) * exp(1j * omega_arm) / omega_arm_b)**2
+    dhdl_sqr = sinc_sqr / armlen**2
+    return dhdl_sqr, sinc_sqr
+
 
 class Strain(nb.Calibration):
     def calc(self):
         dhdl_sqr, sinc_sqr = dhdl(self.freq, self.ifo.Infrastructure.Length)
         return dhdl_sqr
 
-##################################################
+
+############################################################
+# noise sources
+############################################################
+
+#########################
+# quantum
+#########################
 
 class QuantumVacuum(nb.Noise):
     """Quantum Vacuum
@@ -217,6 +231,9 @@ class StandardQuantumLimit(nb.Noise):
     def calc(self, mirror):
         return 8 * const.hbar / (mirror.ETM.MirrorMass * (2 * np.pi * self.freq) ** 2)
 
+#########################
+# seismic
+#########################
 
 class Seismic(nb.Noise):
     """Seismic
@@ -243,6 +260,10 @@ class Seismic(nb.Noise):
         return n * 4
 
 
+#########################
+# Newtonian
+#########################
+
 class Newtonian(nb.Noise):
     """Newtonian Gravity
 
@@ -300,6 +321,10 @@ class NewtonianInfrasound(nb.Noise):
         return n * 2
 
 
+#########################
+# suspension thermal
+#########################
+
 class SuspensionThermal(nb.Noise):
     """Suspension Thermal
 
@@ -316,6 +341,151 @@ class SuspensionThermal(nb.Noise):
         return n * 4
 
 
+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
+
+
+#########################
+# coating thermal
+#########################
+
 class CoatingBrownian(nb.Noise):
     """Coating Brownian
 
@@ -363,6 +533,10 @@ class CoatingThermoOptic(nb.Noise):
         return (nITM + nETM) * 2
 
 
+#########################
+# substrate thermal
+#########################
+
 class ITMThermoRefractive(nb.Noise):
     """ITM Thermo-Refractive
 
@@ -419,6 +593,11 @@ class SubstrateThermoElastic(nb.Noise):
         return (nITM + nETM) * 2
 
 
+#########################
+# residual gas
+#########################
+
+
 class ExcessGas(nb.Noise):
     """Excess Gas
 

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)