Koji's restructured quad suspension code

These changes improve the modeling of suspensions that have a
temperature gradient along the fibers.  Each stage is split into an
upper and a lower joint, which can have different material properties
assigned in the ifo structure.  The temperature of a joint is that of
the mass to which it is attached.  (However, temperature gradients are
not yet supported within the top or bottom stages.)  The average of
the upper and lower joint dissipation is used to compute the noise.

gwinc/ifo/noises.py

@@ -141,21 +141,9 @@ def precomp_suspension(f, ifo):
         ifo.Suspension.VHCoupling = Struct()
         ifo.Suspension.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
 
-    if ifo.Suspension.Type == 'BQuad':
-        material = 'Silicon'
-    else:
-        material = 'Silica'
-    hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension, material)
-    try:
-        # full TF (conventional)
-        ifo.Suspension.hForce = hForce.fullylossy
-        ifo.Suspension.vForce = vForce.fullylossy
-        # TFs with each stage lossy
-        ifo.Suspension.hForce_singlylossy = hForce.singlylossy
-        ifo.Suspension.vForce_singlylossy = vForce.singlylossy
-    except:
-        ifo.Suspension.hForce = hForce
-        ifo.Suspension.vForce = vForce
+    hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension)
+    ifo.Suspension.hForce = hForce
+    ifo.Suspension.vForce = vForce
     ifo.Suspension.hTable = hTable
     ifo.Suspension.vTable = vTable
 

gwinc/noise/suspensionthermal.py

@@ -5,11 +5,12 @@ from __future__ import division
 import numpy as np
 from numpy import pi, imag
 
-from .. import const
+from ..const import kB
+from ..suspension import getJointParams
 
 
 def suspension_thermal(f, sus):
-    """Suspention thermal noise for a single suspended test mass
+    """Suspension thermal noise.
 
     :f: frequency array in Hz
     :sus: gwinc suspension structure
@@ -17,67 +18,45 @@ def suspension_thermal(f, sus):
     :returns: displacement noise power spectrum at :f:
 
     :Temp: must either be set for each stage individually, or globally
-    in :sus.Temp:.  The latter will be preferred if specified, so if
-    you wish to use per-stage tempurature you must remove :sus.Temp:.
+    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.
 
     """
-    # Assign Physical Constants
-    kB = const.kB
 
     # and vertical to beamline coupling angle
     theta = sus.VHCoupling.theta
 
-    noise = np.zeros((1, f.size))
+    noise = np.zeros_like(f)
 
-    # if the temperature is uniform along the suspension
-    if 'Temp' in sus:
-        ##########################################################
-        # Suspension TFs
-        ##########################################################
+    ##########################################################
+    # Suspension TFs
+    ##########################################################
 
-        hForce = sus.hForce
-        vForce = sus.vForce
+    hForce = sus.hForce
+    vForce = sus.vForce
 
-        ##########################################################
-        # Thermal Noise Calculation
-        ##########################################################
-
-        # 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 = hForce + theta**2 * vForce
-
-        # thermal noise (m^2/Hz) for one suspension
-        w = 2*pi*f
-        noise = 4 * kB * sus.Temp * abs(imag(dxdF)) / w
-
-    # if the temperature is set for each suspension stage
-    else:
-        ##########################################################
-        # Suspension TFs
-        ##########################################################
-
-        hForce = sus.hForce_singlylossy[:, :]
-        vForce = sus.vForce_singlylossy[:, :]
-
-        ##########################################################
-        # Thermal Noise Calculation
-        ##########################################################
-
-        dxdF = np.zeros(hForce.shape, dtype=complex)
-        for n, stage in enumerate(reversed(sus.Stage)):
-            # add up the contribution from each stage
+    ##########################################################
+    # Thermal Noise Calculation
+    ##########################################################
 
+    # 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
-            dxdF[n, :] = hForce[n, :] + theta**2 * vForce[n, :]
-
-            # thermal noise (m^2/Hz) for one suspension
-            w = 2*pi*f
-            noise += 4 * kB * stage.Temp * abs(imag(dxdF[n, :])) / w
+            m = 2*n + isLower
+            dxdF[m, :] = hForce[m, :] + theta**2 * vForce[m, :]
+            noise += 4 * kB * Temp * abs(imag(dxdF[m, :])) / w
 
-    return np.squeeze(noise)
+    # thermal noise (m^2/Hz) for one suspension
+    return noise