Commit ff20cbfd authored by Christopher Wipf's avatar Christopher Wipf Committed by Christopher Wipf

Additional suspension code restructuring and bug fixes

Responding to some MR !61 comments:
- docstrings were improved
- most bottom-stage special cases were rolled into new helper functions (`wireGeometry`, `wireTELoss`, `continuumWireK{h,v}`)
- suspensionthermal.py uses `zeros_like` where appropriate
- comment arithmetic was corrected in the Voyager parameters
- `TempLower` and `TempUpper` are returned correctly by `getJointParams`
- `material` kw argument to suspQuad was removed (again)
- `fullylossy`/`singlylossy` outputs was removed (again)
- dilution factor was checked (I think it's correct: see ref GG sec. 2.5)
- vertical spring constant from ref GG was checked (I think it's correct: GG assumed a 2-wire suspension)

Support was added for top stage temperature gradients.  sus.Temp, if defined, is assumed to be the temperature of the topmost joint.

Thermal noise is now computed by summing over the noise of each individual joint (no more "temperature regions").

Support for the stage.FiberRadius parameter has been removed.  None of the canonical IFOs use this.  If anyone ever used it, they can switch to stage.WireRadius instead.

The "questionable approximation" removed in the previous commit was reverted.  It should now agree with what is in master, so that tests pass whenever the suspension temperature is uniform.  This was done to keep the MR focused on restructuring and temperature-gradient related issues.  This approximation (and any other separable issues) can be addressed in another MR.

This commit does not yet deal with Kevin's comment on the splitting of the bottom stage noise.  I would like to find a way to address that while also allowing for a temperature gradient in the bottom stage.

I suggest opening new bugs for Kevin's comments on the vertical TE geometric factors and tapered fiber approximations.
parent 20306fd9
Pipeline #109806 failed with stage
in 1 minute and 12 seconds
......@@ -102,6 +102,8 @@ Suspension:
# This corresponds to 1.14e4 N/m. There are 4 blades.
# Silicon Blade for max 300MPa stress has 9.6mm sag under 50kg load.
# This corresponds to 1.14e4 N/m. There are 4 blades.
# Silicon Blade (40cmx8cm) for max 100MPa stress has 7.5mm sag under 50kg load.
# This corresponds to 6.5e4 N/m. There are 4 blades.
#
# Wire radii scaled from the aLIGO values
# according to the suspended mass by each stage
......@@ -118,7 +120,7 @@ Suspension:
Length: 0.7824 # m; sus_param(1)
Temp: 123.0
Dilution: .nan
K: 2.6e5 # N/m; 1.14e4*4
K: 2.6e5 # N/m; 6.5e4*4
WireRadius: .nan
Blade: 12e-3 # blade thickness
NWires: 4
......
......@@ -142,21 +142,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
......
......@@ -3,9 +3,10 @@
'''
from __future__ import division
import numpy as np
from numpy import pi, imag, zeros
from numpy import pi, imag, zeros_like
from .. import const
from ..const import kB
from ..suspension import getJointParams
def susptherm(f, sus):
......@@ -17,20 +18,18 @@ def susptherm(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 = zeros((1, f.size))
noise = zeros_like(f)
##########################################################
# Suspension TFs
......@@ -43,21 +42,21 @@ def susptherm(f, sus):
# Thermal Noise Calculation
##########################################################
dxdF = zeros(hForce.shape, dtype=complex)
for n, stage in enumerate(reversed(sus.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[n, :] = hForce[n, :] + theta**2 * vForce[n, :]
# thermal noise (m^2/Hz) for one suspension
if 'Temp' in stage:
T = stage.Temp
else:
T = sus.Temp
w = 2*pi*f
noise += 4 * kB * T * abs(imag(dxdF[n, :])) / w
return np.squeeze(noise)
# 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 = 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
# thermal noise (m^2/Hz) for one suspension
return noise
This diff is collapsed.
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment