Residual gas sub-budgets and updates

Adds residual gas damping and additional molecular species to the residual gas calculations. The gas scattering and gas damping are now tracked separately for each species in an excess gas sub-budget. The damping calculation is done in the infinite volume limit by default but adds corrections for squeezed film damping if specified in the yaml file. The aLIGO excess gas sub-budget is shown below as an example.

The new yaml file fields are explained in this comment. The yaml files in this merge request are simpler and none of the interferometers include squeezed film damping.

For aLIGO, Aplus, and Voyager, H2, N2, and H2O are considered. The numbers for aLIGO are taken from the link in this comment. The numbers for Aplus and Voyager are copied from aLIGO and should be checked. The Cosmic Explorer interferometers additionally consider CO2 and Ar.

Closes #63 (closed).

gwinc/noise/residualgas.py

@@ -5,18 +5,21 @@ from __future__ import division
 from numpy import sqrt, log, pi
 
 from .. import const
+from .. import Struct
 
 
-def residual_gas_cavity(f, ifo):
-    """Residual gas noise strain spectrum in an evacuated cavity.
+def residual_gas_scattering_arm(f, ifo, cavity, species):
+    """Residual gas noise strain spectrum due to scattering from one arm
 
     Noise caused by the passage of residual gas molecules through the
-    laser beams in a cavity.
+    laser beams in one arm cavity due to scattering.
 
     :f: frequency array in Hz
     :ifo: gwinc IFO structure
+    :cavity: arm cavity structure
+    :species: molecular species structure
 
-    :returns: arm displacement noise power spectrum at :f:
+    :returns: arm strain noise power spectrum at :f:
 
     The method used here is presented by Rainer Weiss, Micheal
     E. Zucker, and Stanley E. Whitcomb in their paper Optical
@@ -28,26 +31,19 @@ def residual_gas_cavity(f, ifo):
     expansion of exp, to speed it up.
 
     """
-    Lambda = ifo.Laser.Wavelength
     L = ifo.Infrastructure.Length
     kT = ifo.Infrastructure.Temp * const.kB
-    P = ifo.Infrastructure.ResidualGas.pressure
-    M = ifo.Infrastructure.ResidualGas.mass
-    alpha = ifo.Infrastructure.ResidualGas.polarizability
-    R1 = ifo.Optics.Curvature.ITM
-    R2 = ifo.Optics.Curvature.ETM
-
-    rho = P /(kT)                                   # number density of Gas
-    v0 = sqrt(2*kT / M)                             # mean speed of Gas
-
-    g1 = 1 - L/R1                                   # first resonator g-parameter of the ITM
-    g2 = 1 - L/R2                                   # second resonator g-parameter of the ETM
-    waist = L * Lambda / pi
-    waist = waist * sqrt(((g1*g2)*(1-g1*g2))/((g1+g2-2*g1*g2)**2))
-    waist = sqrt(waist)                             # Gaussian beam waist size
-    zr = pi*waist**2 / Lambda                       # Rayleigh range
-    z1 = -((g2*(1-g1))/(g1+g2-2*g1*g2))*L           # location of ITM relative to the waist
-    z2 =  ((g1*(1-g2))/(g1+g2-2*g1*g2))*L           # location of ETM relative to the waist
+    P = species.BeamtubePressure
+    M = species.mass
+    alpha = species.polarizability
+
+    rho = P / (kT)                   # number density of Gas
+    v0 = sqrt(2*kT / M)              # mean speed of Gas
+
+    waist = cavity.w0                # Gaussian beam waist size
+    zr = cavity.zr                   # Rayleigh range
+    z1 = -cavity.zBeam_ITM           # location of ITM relative to the waist
+    z2 = cavity.zBeam_ETM            # location of ETM relative to the waist
 
     # The exponential of Eq. 1 of P940008 is expanded to first order; this
     # can be integrated analytically
@@ -60,3 +56,78 @@ def residual_gas_cavity(f, ifo):
     zint[zint < 0] = 0
 
     return zint
+
+
+def residual_gas_damping_test_mass(f, ifo, species, sustf, squeezed_film):
+    """Noise due to residual gas damping for one test mass
+
+    :f: frequency array in Hz
+    :ifo: gwinc IFO structure
+    :species: molecular species structure
+    :sustf: suspension transfer function structure
+    :squeezed_film: squeezed film damping structure
+
+    :returns: displacement noise
+    """
+    sus = ifo.Suspension
+    if 'Temp' in sus.Stage[0]:
+        kT = sus.Stage[0].Temp * const.kB
+    else:
+        kT = sus.Temp * const.kB
+
+    mass = species.mass
+    radius = ifo.Materials.MassRadius
+    thickness = ifo.Materials.MassThickness
+    thermal_vel = sqrt(kT/mass)  # thermal velocity
+
+    # pressure in the test mass chambers; possibly different from the pressure
+    # in the arms due to outgassing near the test mass
+    pressure = species.ChamberPressure
+
+    # infinite volume viscous damping coefficient for a cylinder
+    # table 1 of https://doi.org/10.1016/j.physleta.2010.06.041
+    beta_inf = pi * radius**2 * pressure/thermal_vel
+    beta_inf *= sqrt(8/pi) * (1 + thickness/(2*radius) + pi/4)
+
+    force_noise = 4 * kT * beta_inf
+
+    # add squeezed film damping if necessary as parametrized by
+    # Eq (5) of http://dx.doi.org/10.1103/PhysRevD.84.063007
+    if squeezed_film.keys():
+        # the excess force noise and diffusion time are specified directly
+        if 'ExcessDamping' in squeezed_film:
+            # ExcessDamping is the ratio of the total gas damping noise at DC
+            # to damping in the infinite volume limit (in amplitude)
+            DeltaS0 = (squeezed_film.ExcessDamping**2 - 1) * force_noise
+            if DeltaS0 < 0:
+                raise ValueError('ExcessDamping must be > 1')
+
+            try:
+                diffusion_time = squeezed_film.DiffusionTime
+            except AttributeError:
+                msg = 'If excess residual gas damping is given a diffusion ' \
+                    + 'time must be specified as well'
+                raise ValueError(msg)
+
+        # if a gap between the test mass and another object is specified
+        # use the approximate model of section IIIA and B
+        elif 'gap' in squeezed_film:
+            gap = squeezed_film.gap
+
+            # Eq (14)
+            diffusion_time = sqrt(pi/2) * radius**2 / (gap * thermal_vel)
+            diffusion_time /= log(1 + (radius/gap)**2)
+
+            # Eq (11) factoring out the low pass cutoff as in (5)
+            DeltaS0 = 4 * kT * pi*radius**2 * pressure * diffusion_time / gap
+
+        else:
+            raise ValueError('Must specify either excess damping or a gap')
+
+        # Eq (5)
+        force_noise += DeltaS0 / (1 + (2*pi*f*diffusion_time)**2)
+
+    # convert force to displacement noise using the suspension susceptibility
+    noise = force_noise * abs(sustf.tst_suscept)**2
+
+    return noise