Commit 73516134 authored by Christopher Wipf's avatar Christopher Wipf

Suspension refactoring

This commit keeps the restructuring work from this series, but
all changes that affect noise curves should now have been reverted.
Corrections that should be reinstated later are flagged with FIXME comments.
susptherm has been re-renamed to suspension_thermal.
parent 7c9ed00a
Pipeline #118442 failed with stage
in 1 minute and 37 seconds
This diff is collapsed.
......@@ -278,7 +278,7 @@ class SuspensionThermal(nb.Noise):
def calc(self):
precomp_suspension(self.freq, self.ifo)
n = noise.suspensionthermal.susptherm(self.freq, self.ifo.Suspension)
n = noise.suspensionthermal.suspension_thermal(self.freq, self.ifo.Suspension)
return n * 4
......
......@@ -9,7 +9,7 @@ from ..const import kB
from ..suspension import getJointParams
def susptherm(f, sus):
def suspension_thermal(f, sus):
"""Suspension thermal noise.
:f: frequency array in Hz
......
from __future__ import division
from numpy import pi, sqrt, sin, cos, tan, real, imag, zeros
from numpy import pi, sqrt, sin, cos, tan, real, imag
import numpy as np
import logging
......@@ -101,6 +101,7 @@ def getJointParams(sus, n):
TempUpper = stages[n-1].Temp
else:
TempUpper = sus.Temp
TempUpper = Temp # FIXME: reproduces the old calculation
##############################
# material parameters
......@@ -114,6 +115,10 @@ def getJointParams(sus, n):
elif 'WireMaterial' in stage:
wireMatUpper = stage.WireMaterial
wireMatLower = stage.WireMaterial
# FIXME: reproduces the old calculation
elif n == last_stage and sus.Type == 'BQuad':
wireMatUpper = 'Silicon'
wireMatLower = 'Silicon'
elif n == last_stage:
wireMatUpper = 'Silica'
wireMatLower = 'Silica'
......@@ -128,6 +133,9 @@ def getJointParams(sus, n):
# support blade (upper joint only)
if 'BladeMaterial' in stage:
bladeMat = stage.BladeMaterial
# FIXME: reproduces the old calculation
elif n == last_stage and sus.Type == 'BQuad':
bladeMat = 'Silicon'
else:
bladeMat = 'MaragingSteel'
......@@ -148,7 +156,7 @@ def wireGeometry(r, N, RibbonThickness=None, TaperedEndRadius=None, **kwargs):
TaperedEndRadius when tapered fibers are used.
Returns the cross-sectional area and moment of inertia,
and the surface to volume ratios (vertical and horizontal).
and the modified surface to volume ratios (vertical and horizontal).
"""
## Usual case: round wire/fiber
......@@ -160,6 +168,7 @@ def wireGeometry(r, N, RibbonThickness=None, TaperedEndRadius=None, **kwargs):
# surface to volume ratio, vertical
# Surface = 2*pi*r*h
# Volume = pi*r^2*h
# Geometrical factor \mu is unity (corresponding to uniform strain energy)
mu_v = 2 / r
# surface to volume ratio, horizontal
# Ref: Gretarsson and Harry, Rev. Sci. Inst. 70 (1999) 4081-4087
......@@ -174,8 +183,9 @@ def wireGeometry(r, N, RibbonThickness=None, TaperedEndRadius=None, **kwargs):
if RibbonThickness is not None:
W = r
t = RibbonThickness
xsect = W * t # cross-sectional area
xII = (W * t**3)/12 # x-sectional moment of inertia
xsect = W * t # cross-sectional area
xsectEnd = xsect # only differs for tapered fiber
xII = (W * t**3)/12 # x-sectional moment of inertia
# surface to volume ratio, vertical
# Surface = 2*(W+t)*h
# Volume = W*t*h
......@@ -192,13 +202,13 @@ def wireGeometry(r, N, RibbonThickness=None, TaperedEndRadius=None, **kwargs):
xsectEnd = pi * r_end**2 # cross-sectional area (for delta_h)
xII = pi * r_end**4 / 4 # x-sectional moment of inertia
mu_h = 4 / r_end # surface to volume ratio, horizontal
# mu_v is unchanged (assumes the middle part of the fiber dominates)
return (xsect, xsectEnd, xII, mu_v, mu_h)
def wireTELoss(w, tension, xsect, xsectEnd, xII, Temp,
alpha, beta, rho, C, K, Y,
RibbonThickness=None, TaperedEndRadius=None, **kwargs):
def wireTELoss(w, tension, xsectEnd, xII, Temp, alpha, beta, rho, C, K, Y, xsect,
RibbonThickness=None, **kwargs):
"""Thermoelastic calculations for wires
Repeated for upper and lower joint of each stage.
......@@ -212,12 +222,15 @@ def wireTELoss(w, tension, xsect, xsectEnd, xII, Temp,
"""
# horizontal TE time constant, wires
# The constant 7.37e-2 is 1/(4*q0^2) from eq 12, C. Zener 10.1103/PhysRev.53.90 (1938)
tau = 7.37e-2 * 4 * (rho * C * xsectEnd) / (pi * K)
# FIXME: reproduces the old calculation
# then xsect can be dropped from the argument list
tau = 7.37e-2 * 4 * (rho * C * xsect) / (pi * K)
# deal with ribbon geometry
if RibbonThickness is not None:
t = RibbonThickness
tau_h = (rho * C * t**2) / (K * pi**2)
tau = (rho * C * t**2) / (K * pi**2)
# delta: TE factor
# The first term expresses the cancellation of the CTE and dY/dT
......@@ -226,12 +239,7 @@ def wireTELoss(w, tension, xsect, xsectEnd, xII, Temp,
# Cagnoli G and Willems P 2002 Phys. Rev. B 65 174111
# Cumming A V et al 2009 Class. Quantum Grav. 26 215012
# horizontal delta, wires
delta = (alpha - tension * beta / (xsect * Y))**2 * Y * Temp / (rho * C)
# deal with tapered geometry
if TaperedEndRadius is not None:
# use end xsect for thermo-elastic noise
delta = (alpha - tension * beta / (xsectEnd * Y))**2 * Y * Temp / (rho * C)
delta = (alpha - tension * beta / (xsectEnd * Y))**2 * Y * Temp / (rho * C)
phi_TE = delta * tau * w / (1 + w**2 * tau**2)
return phi_TE
......@@ -299,16 +307,16 @@ def continuumWireKh(w, N, length, tension, xsect, xII, rho, Y, phi):
# = T k (cos(k L) + k delta sin(k L))
# for w -> 0, this reduces to N_w * T * k
khnum = N * tension * k * (coskl + dk * sinkl)
if True: # FIXME: use previous calculation so tests pass
khnum = N * tension * k * (1 + dk**2) * (coskl + dk * sinkl)
# FIXME: reproduces the old calculation
khnum = N * tension * k * (1 + dk**2) * (coskl + dk * sinkl)
# denominator, horiz spring constant
# D after equation 8 in GG
# D = sin(k L) - 2 k delta cos(k L)
# for w -> 0, this reduces to k (L - 2 delta)
khden = sinkl - 2 * dk * coskl
if True: # FIXME: use previous calculation so tests pass
khden = (1 - dk**2) * sinkl - 2 * dk * coskl
# FIXME: reproduces the old calculation
khden = (1 - dk**2) * sinkl - 2 * dk * coskl
return khnum/khden
......@@ -324,11 +332,14 @@ def continuumWireKv(w, N, length, xsect, xsectEnd, rho, Y, phi,
kv = N * xsect * Y * k / tan(k * length)
# deal with tapered geometry
if TaperedEndLength is not None:
l_end = 2 * TaperedEndLength
l_mid = length - l_end
l_end = TaperedEndLength
l_mid = length - 2*l_end
kv_mid = N * xsect * Y * k / tan(k * l_mid)
kv_end = N * xsectEnd * Y * k / tan(k * l_end)
kv = kv_mid * kv_end / (kv_mid + kv_end)
kv = 1/(2/kv_end + 1/kv_mid)
# FIXME: reproduces the old calculation
kv_end = N * xsectEnd * Y * k / tan(k * 2*l_end)
kv = 1/(1/kv_end + 1/kv_mid)
return kv
......@@ -385,12 +396,12 @@ def suspQuad(f, sus):
##############################
# Complex spring constants
khr = zeros([len(stages), len(w)])
khi = zeros([len(stages), 2, len(w)])
kh = zeros([len(stages), len(w)], dtype=np.complex_)
kvr = zeros([len(stages), len(w)])
kvi = zeros([len(stages), 2, len(w)])
kv = zeros([len(stages), len(w)], dtype=np.complex_)
khr = np.zeros([len(stages), len(w)])
khi = np.zeros([len(stages), 2, len(w)])
kh = np.zeros([len(stages), len(w)], dtype=np.complex_)
kvr = np.zeros([len(stages), len(w)])
kvi = np.zeros([len(stages), 2, len(w)])
kv = np.zeros([len(stages), len(w)], dtype=np.complex_)
for n, stage in enumerate(stages):
......@@ -481,8 +492,8 @@ def suspQuad(f, sus):
# (mu_h*ds_w): surface loss (brownian motion)
# nominally this term becomes zero for steel wires as ds_w is zero
# The last term is for TE loss
phih_TE = wireTELoss(w, tension, xsect, xsectEnd, xII,
Temp, alpha_w, beta_w, rho_w, C_w, K_w, Y_w,
phih_TE = wireTELoss(w, tension, xsectEnd, xII, Temp,
alpha_w, beta_w, rho_w, C_w, K_w, Y_w, xsect,
**wireShape)
phih = phi_w * (1 + mu_h * ds_w) + phih_TE
......@@ -493,6 +504,7 @@ def suspQuad(f, sus):
if not isLower:
# vertical loss factor, blades
# The bulk loss term and the TE loss term
# wire loss and stiffness are neglected here
phiv_TE = bladeTELoss(w, t_b, Temp,
alpha_b, beta_b, rho_b, C_b, K_b, Y_b)
phiv_blade = phi_b + phiv_TE
......@@ -504,6 +516,7 @@ def suspQuad(f, sus):
if n == last_stage:
# loss factor, vertical (no blades)
# wire vertical thermoelastic loss is negligible (Saulson 1990)
phiv = phi_w * (1 + mu_v * ds_w)
kv4 = continuumWireKv(w, N_w, length, xsect, xsectEnd,
......@@ -531,13 +544,13 @@ def suspQuad(f, sus):
# Calculate horizontal/vertical TFs turning on the loss of each
# joint one by one, for the thermal noise calculation
hForce = zeros([2*len(stages), len(w)], dtype=complex)
vForce = zeros([2*len(stages), len(w)], dtype=complex)
hForce = np.zeros([2*len(stages), len(w)], dtype=complex)
vForce = np.zeros([2*len(stages), len(w)], dtype=complex)
for m in range(2*len(stages)):
# turn on only the loss of the current joint
lossy_region_lower = zeros((len(stages), 1))
lossy_region_upper = zeros((len(stages), 1))
lossy_region_lower = np.zeros((len(stages), 1))
lossy_region_upper = np.zeros((len(stages), 1))
n = int(m/2) # stage number
if m % 2:
lossy_region_upper[n] = 1
......
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