suspension: update suspQuad function from matgwinc supporting stage iteration

This is a rote copy of the new suspQuad code in matgwinc that iterates over
the stages, and supports tapered fibers:

https://git.ligo.org/gwinc/matgwinc/commit/a59f01192a1ffc1644a45be536e31bc0eb78bd7c#9cb6e594e14c1d61dbb3c91549166a2d48db194f

This allows for setting per-stage temperature and wire and blade properties.

This version differs slightly in some generalizations about the number of
stages, although it doesn't quite yet support an arbitrary numer.

gwinc/ifo/Voyager.yaml

@@ -78,12 +78,6 @@ Seismic:
 
 Suspension:
   Type: 'BQuad'
-  # Suspension fiber temperatures [TOP UIM PUM TST]
-  Temp:
-    - 300.0
-    - 300.0
-    - 300.0
-    - 123.0
   VHCoupling:
     theta: 1e-3 # vertical-horizontal x-coupling
   FiberType: 'Ribbon'
@@ -101,6 +95,7 @@ Suspension:
     #susmat = loadmat('CryogenicLIGO/QuadModel/quad_optimized_masses_for_PUM_with_springs.mat')
     - Mass: 200.0 # kg; susmat['testmass_mass'][0,0]
       Length: 0.4105 # m
+      Temp: 123.0
       Dilution: .nan
       K: .nan
       WireRadius: .nan
@@ -108,6 +103,7 @@ Suspension:
       NWires: 4
     - Mass: 65.9 # kg; susmat['PUMmass'][0,0]
       Length: 0.4105 # m
+      Temp: 300.0
       Dilution: 106.0
       K: 5200.0 # N/m; vertical spring constant
       WireRadius: 310e-6
@@ -115,6 +111,7 @@ Suspension:
       NWires: 4
     - Mass: 87.6 # kg; susmat['UIMmass'][0,0]
       Length: 0.4105 # m
+      Temp: 300.0
       Dilution: 80.0
       K: 3900.0 # N/m; vertical spring constant
       WireRadius: 350e-6
@@ -122,6 +119,7 @@ Suspension:
       NWires: 4
     - Mass: 116.5 # kg; susmat['topmass_mass'][0,0]
       Length: 0.4105 # m
+      Temp: 300.0
       Dilution: 87.0
       K: 3400.0 # N/m; vertical spring constant
       WireRadius: 520e-6

gwinc/noise/suspensionthermal.py

@@ -7,6 +7,11 @@ import scipy.constants
 def susptherm(f, ifo):
     """Suspention thermal noise.
 
+    'Temp' must either be set for each stage individually, or globally
+    in ifo.Suspension.Temp.  The latter will be preferred if
+    specified, so if you wish to use per-stage tempurature you must
+    remove Suspension.Temp.
+
     Assumes suspension transfer functions and V-H coupling have been
     pre-calculated and populated into the relevant `ifo` struct
     fields.
@@ -14,7 +19,6 @@ def susptherm(f, ifo):
     """
     # Assign Physical Constants
     kB   = scipy.constants.k
-    Temp = ifo.Suspension.Temp
 
     # and vertical to beamline coupling angle
     theta = ifo.Suspension.VHCoupling.theta
@@ -22,10 +26,11 @@ def susptherm(f, ifo):
     noise = zeros((1, f.size))
 
     # if the temperature is uniform along the suspension
-    if np.isscalar(Temp) or len(Temp) == 1:
+    if 'Temp' in ifo.Suspension:
         ##########################################################
         # Suspension TFs
         ##########################################################
+
         hForce = ifo.Suspension.hForce
         vForce = ifo.Suspension.vForce
 
@@ -40,32 +45,33 @@ def susptherm(f, ifo):
 
         # thermal noise (m^2/Hz) for one suspension
         w = 2*pi*f
-        noise = 4 * kB * Temp * abs(imag(dxdF)) / w
+        noise = 4 * kB * ifo.Suspension.Temp * abs(imag(dxdF)) / w
 
     # if the temperature is set for each suspension stage
     else:
         ##########################################################
         # Suspension TFs
         ##########################################################
-        hForce = ifo.Suspension.hForce_singlylossy[:,:]
-        vForce = ifo.Suspension.vForce_singlylossy[:,:]
+
+        hForce = ifo.Suspension.hForce_singlylossy[:, :]
+        vForce = ifo.Suspension.vForce_singlylossy[:, :]
 
         ##########################################################
         # Thermal Noise Calculation
         ##########################################################
 
         dxdF = zeros(hForce.shape, dtype=complex)
-        for ii in range(len(Temp)):
+        for n, stage in enumerate(ifo.Suspension.Stage):
             # add up the contribution from each stage
 
-            # 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[ii,:] = hForce[ii,:] + theta**2 * vForce[ii,:]
+            # 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 * Temp[ii] * abs(imag(dxdF[ii,:])) / w
+            noise += 4 * kB * stage.Temp * abs(imag(dxdF[n, :])) / w
 
     # 4 masses, turn into gravitational wave strain
     noise *= 4 * ifo.gwinc.dhdl_sqr