README.md

@@ -113,13 +113,22 @@ For custom plotting, parameter optimization, etc. all functionality can be
 accessed directly through the `gwinc` library interface:
 ```python
 >>> import gwinc
->>> import numpy as np
->>> freq = np.logspace(1, 3, 1000)
->>> budget = gwinc.load_budget('aLIGO', freq)
+>>> budget = gwinc.load_budget('aLIGO')
 >>> trace = budget.run()
 >>> fig = gwinc.plot_budget(trace)
 >>> fig.show()
 ```
+A default frequency array is used, but alternative frequencies can be
+provided to `load_budget()` either in the form of a numpy array:
+```python
+>>> import numpy as np
+>>> freq = np.logspace(1, 3, 1000)
+>>> budget = gwinc.load_budget('aLIGO', freq=freq)
+```
+or frequency specification string ('FLO:[NPOINTS:]FHI'):
+```
+>>> budget = gwinc.load_budget('aLIGO', freq='10:1000:1000')
+```
 
 The `load_budget()` function takes most of the same inputs as the
 command line interface (e.g. IFO names, budget module paths, YAML
@@ -134,6 +143,15 @@ properties.  The budget sub-traces are available through a dictionary
 (`trace['QuantumVacuum']`) interface and via attributes
 (`trace.QuantumVacumm`).
 
+The budget `freq` and `ifo` attributes can be updated at run time by
+passing them as keyword arguments to the `run()` method:
+```python
+>>> budget = load_budget('aLIGO')
+>>> freq = np.logspace(1, 3, 1000)
+>>> ifo = Struct.from_file('/path/to/ifo_alt.yaml')
+>>> trace = budget.run(freq=freq, ifo=ifo)
+```
+
 
 ## noise functions
 
@@ -340,12 +358,6 @@ If a budget module defined as a package includes an `ifo.yaml`
 [parameter file](#parameter-files) in the package directory, the
 `load_budget()` function will automatically load the YAML data into an
 `ifo` `gwinc.Struct` and assign it to the `budget.ifo` attribute.
-Alternate ifos can be specified at run time:
-```python
-budget = load_budget('/path/to/MyBudget', freq)
-ifo = Struct.from_file('/path/to/MyBudget.ifo')
-trace = budget.run(ifo=ifo)
-```
 
 The IFOs included in `gwinc.ifo` provide examples of the use of the
 budget interface (e.g. [gwinc.ifo.aLIGO](gwinc/ifo/aLIGO)).

gwinc/__init__.py

@@ -14,6 +14,31 @@ from .plot import plot_noise
 logger = logging.getLogger('gwinc')
 
 
+DEFAULT_FREQ = '5:3000:6000'
+
+
+def freq_from_spec(spec=None):
+    """logarithmicly spaced frequency array, based on specification string
+
+    Specification string should be of form 'START:[NPOINTS:]STOP'.  If
+    `spec` is an array, then the array is returned as-is, and if it's
+    None the DEFAULT_FREQ specification is used.
+
+    """
+    if isinstance(spec, np.ndarray):
+        return spec
+    elif spec is None:
+        spec = DEFAULT_FREQ
+    fspec = spec.split(':')
+    if len(fspec) == 2:
+        fspec = fspec[0], DEFAULT_FREQ.split(':')[1], fspec[1]
+    return np.logspace(
+        np.log10(float(fspec[0])),
+        np.log10(float(fspec[2])),
+        int(fspec[1]),
+    )
+
+
 def load_module(name_or_path):
     """Load module from name or path.
 
@@ -43,19 +68,24 @@ def load_budget(name_or_path, freq=None):
 
     Accepts either the name of a built-in canonical budget (see
     gwinc.IFOS), the path to a budget package (directory) or module
-    (ending in .py), or the path to an IFO struct definition file (see
+    (ending in .py), or the path to an IFO Struct definition file (see
     gwinc.Struct).
 
     If the budget is a package directory which includes an 'ifo.yaml'
     file the ifo Struct will be loaded from that file and assigned to
-    the budget.ifo attribute.  If a struct definition file is provided
-    the base aLIGO budget definition and ifo will be assumed.
+    the budget.ifo attribute.  If a Struct definition file is provided
+    the base aLIGO budget definition will be assumed.
+
+    Returns an instantiated Budget object.  If a frequency array or
+    frequency specification string (see `freq_from_spec()`) is
+    provided, the budget will be instantiated with the provided array.
+    If a frequency array is not provided and the Budget class
+    definition includes a `freq` attribute defining either an array or
+    specification string, then that array will be used.  Otherwise a
+    default array will be provided (see DEFAULT_FREQ).
 
-    Returns a Budget object instantiated with the provided frequency
-    array (if specified), and with any included ifo assigned as an
-    object attribute.  If a frequency array is not provided and the
-    budget class does not define it's own, the frequency array must be
-    provided at budget update() or run() time.
+    Any included ifo will be assigned as an attribute to the returned
+    Budget object.
 
     """
     ifo = None
@@ -85,6 +115,9 @@ def load_budget(name_or_path, freq=None):
     logger.info("loading module {}...".format(modname))
     mod, modpath = load_module(modname)
     Budget = getattr(mod, bname)
+    if freq is None:
+        freq = getattr(Budget, 'freq', None)
+    freq = freq_from_spec(freq)
     ifopath = os.path.join(modpath, 'ifo.yaml')
     if not ifo and os.path.exists(ifopath):
         ifo = Struct.from_file(ifopath)

gwinc/__main__.py

@@ -3,9 +3,15 @@ import os
 import signal
 import logging
 import argparse
-import numpy as np
 
-from . import IFOS, load_budget, plot_budget, logger
+from . import (
+    IFOS,
+    DEFAULT_FREQ,
+    freq_from_spec,
+    load_budget,
+    plot_budget,
+    logger,
+)
 from . import io
 
 logger.setLevel(os.getenv('LOG_LEVEL', 'WARNING').upper())
@@ -49,10 +55,11 @@ e.g.:
   gwinc --fom range:m1=20,m2=20 ...
 
 See the inspiral_range package documentation for details.
+NOTE: The range will be calculated with the supplied frequency array,
+and may therefore be inaccurate if a restricted array is specified.
 """
 
 IFO = 'aLIGO'
-FREQ = '5:3000:6000'
 FOM = 'range:m1=1.4,m2=1.4'
 
 parser = argparse.ArgumentParser(
@@ -61,7 +68,7 @@ parser = argparse.ArgumentParser(
     formatter_class=argparse.RawDescriptionHelpFormatter)
 parser.add_argument(
     '--freq', '-f', metavar='FLO:[NPOINTS:]FHI',
-    help="logarithmic frequency array specification in Hz [{}]".format(FREQ))
+    help="logarithmic frequency array specification in Hz [{}]".format(DEFAULT_FREQ))
 parser.add_argument(
     '--ifo', '-o', metavar='PARAM=VAL',
     #nargs='+', action='extend',
@@ -97,17 +104,6 @@ parser.add_argument(
     help="IFO name or path")
 
 
-def freq_from_spec(spec):
-    fspec = spec.split(':')
-    if len(fspec) == 2:
-        fspec = fspec[0], FREQ.split(':')[1], fspec[1]
-    return np.logspace(
-        np.log10(float(fspec[0])),
-        np.log10(float(fspec[2])),
-        int(fspec[1]),
-    )
-
-
 def main():
     signal.signal(signal.SIGINT, signal.SIG_DFL)
 
@@ -129,15 +125,12 @@ def main():
         plot_style = trace.plot_style
 
     else:
-        budget = load_budget(args.IFO)
+        try:
+            freq = freq_from_spec(args.freq)
+        except IndexError:
+            parser.exit(2, "Improper frequency specification: {}\n".format(args.freq))
+        budget = load_budget(args.IFO, freq=freq)
         ifo = budget.ifo
-        if args.freq:
-            try:
-                freq = freq_from_spec(args.freq)
-            except IndexError:
-                parser.exit(2, "Improper frequency specification: {}\n".format(args.freq))
-        else:
-            freq = getattr(budget, 'freq', freq_from_spec(FREQ))
         plot_style = getattr(budget, 'plot_style', {})
         trace = None
 

gwinc/ifo/Aplus/__init__.py

@@ -2,6 +2,26 @@ from gwinc.ifo.noises import *
 from gwinc.ifo import PLOT_STYLE
 
 
+class QuantumVacuum(nb.Budget):
+    """Quantum Vacuum
+
+    """
+    style = dict(
+        label='Quantum Vacuum',
+        color='#ad03de',
+    )
+
+    noises = [
+        QuantumVacuumAS,
+        QuantumVacuumArm,
+        QuantumVacuumSEC,
+        QuantumVacuumFilterCavity,
+        QuantumVacuumInjection,
+        QuantumVacuumReadout,
+        QuantumVacuumQuadraturePhase,
+    ]
+
+
 class Aplus(nb.Budget):
 
     name = 'A+'

gwinc/ifo/CE1/__init__.py

@@ -2,6 +2,26 @@ from gwinc.ifo.noises import *
 from gwinc.ifo import PLOT_STYLE
 
 
+class QuantumVacuum(nb.Budget):
+    """Quantum Vacuum
+
+    """
+    style = dict(
+        label='Quantum Vacuum',
+        color='#ad03de',
+    )
+
+    noises = [
+        QuantumVacuumAS,
+        QuantumVacuumArm,
+        QuantumVacuumSEC,
+        QuantumVacuumFilterCavity,
+        QuantumVacuumInjection,
+        QuantumVacuumReadout,
+        QuantumVacuumQuadraturePhase,
+    ]
+
+
 class Newtonian(nb.Budget):
     """Newtonian Gravity
 

gwinc/ifo/CE2/__init__.py

@@ -2,6 +2,26 @@ from gwinc.ifo.noises import *
 from gwinc.ifo import PLOT_STYLE
 
 
+class QuantumVacuum(nb.Budget):
+    """Quantum Vacuum
+
+    """
+    style = dict(
+        label='Quantum Vacuum',
+        color='#ad03de',
+    )
+
+    noises = [
+        QuantumVacuumAS,
+        QuantumVacuumArm,
+        QuantumVacuumSEC,
+        QuantumVacuumFilterCavity,
+        QuantumVacuumInjection,
+        QuantumVacuumReadout,
+        QuantumVacuumQuadraturePhase,
+    ]
+
+
 class Newtonian(nb.Budget):
     """Newtonian Gravity
 

gwinc/ifo/CE2/ifo.yaml

@@ -321,6 +321,7 @@ Squeezer:
   AmplitudedB: 15                  # SQZ amplitude [dB]
   InjectionLoss: 0.02              # power loss to sqz
   SQZAngle: 0                      # SQZ phase [radians]
+  LOAngleRMS: 10e-3               # quadrature noise [radians]
 
   # Parameters for frequency dependent squeezing
   FilterCavity:

gwinc/ifo/Voyager/__init__.py

@@ -2,6 +2,25 @@ from gwinc.ifo.noises import *
 from gwinc.ifo import PLOT_STYLE
 
 
+class QuantumVacuum(nb.Budget):
+    """Quantum Vacuum
+
+    """
+    style = dict(
+        label='Quantum Vacuum',
+        color='#ad03de',
+    )
+
+    noises = [
+        QuantumVacuumAS,
+        QuantumVacuumArm,
+        QuantumVacuumSEC,
+        QuantumVacuumFilterCavity,
+        QuantumVacuumInjection,
+        QuantumVacuumReadout,
+    ]
+
+
 class Voyager(nb.Budget):
 
     name = 'Voyager'

gwinc/ifo/aLIGO/__init__.py

@@ -2,6 +2,23 @@ from gwinc.ifo.noises import *
 from gwinc.ifo import PLOT_STYLE
 
 
+class QuantumVacuum(nb.Budget):
+    """Quantum Vacuum
+
+    """
+    style = dict(
+        label='Quantum Vacuum',
+        color='#ad03de',
+    )
+
+    noises = [
+        QuantumVacuumAS,
+        QuantumVacuumArm,
+        QuantumVacuumSEC,
+        QuantumVacuumReadout,
+    ]
+
+
 class aLIGO(nb.Budget):
 
     name = 'Advanced LIGO'

gwinc/ifo/noises.py

@@ -9,16 +9,10 @@ from .. import nb
 from .. import noise
 from .. import suspension
 
-##################################################
-
 
-def coating_thickness(materials, optic):
-    if 'CoatLayerOpticalThickness' in optic:
-        return optic['CoatLayerOpticalThickness']
-    T = optic.Transmittance
-    dL = optic.CoatingThicknessLown
-    dCap = optic.CoatingThicknessCap
-    return noise.coatingthermal.getCoatDopt(materials, T, dL, dCap=dCap)
+############################################################
+# helper functions
+############################################################
 
 
 def mirror_struct(ifo, tm):
@@ -33,7 +27,13 @@ def mirror_struct(ifo, tm):
     # it would otherwise be stored by reference)
     mirror = copy.deepcopy(ifo.Materials)
     optic = ifo.Optics.get(tm)
-    mirror.Coating.dOpt = coating_thickness(ifo.Materials, optic)
+    if 'CoatLayerOpticalThickness' in optic:
+        mirror.Coating.dOpt = optic['CoatLayerOpticalThickness']
+    else:
+        T = optic.Transmittance
+        dL = optic.CoatingThicknessLown
+        dCap = optic.CoatingThicknessCap
+        mirror.Coating.dOpt = noise.coatingthermal.getCoatDopt(mirror, T, dL, dCap=dCap)
     mirror.update(optic)
     mirror.MassVolume = pi * mirror.MassRadius**2 * mirror.MassThickness
     mirror.MirrorMass = mirror.MassVolume * mirror.Substrate.MassDensity
@@ -64,7 +64,14 @@ def arm_cavity(ifo):
     z2 = L * g1 * (1 - g2) / gden
 
     # waist, input, output
-    return w0, w1, w2
+    cavity = Struct()
+    cavity.w0 = w0
+    cavity.wBeam_ITM = w1
+    cavity.wBeam_ETM = w2
+    cavity.zr = zr
+    cavity.zBeam_ITM = z1
+    cavity.zBeam_ETM = z2
+    return cavity
 
 
 def ifo_power(ifo, PRfixed=True):
@@ -108,7 +115,14 @@ def ifo_power(ifo, PRfixed=True):
     if pbs > pbsl:
         logger.warning('P_BS exceeds BS Thermal limit!')
 
-    return pbs, parm, finesse, prfactor, Tpr
+    power = Struct()
+    power.pbs = pbs
+    power.parm = parm
+    power.finesse = finesse
+    power.gPhase = finesse * 2/np.pi
+    power.prfactor = prfactor
+    power.Tpr = Tpr
+    return power
 
 
 def dhdl(f, armlen):
@@ -153,40 +167,80 @@ def precomp_suspension(f, ifo):
     return pc
 
 
-def precomp_mirror(f, ifo):
+
+def precomp_quantum(f, ifo):
     pc = Struct()
-    pc.ITM = mirror_struct(ifo, 'ITM')
-    pc.ETM = mirror_struct(ifo, 'ETM')
+    mirror_mass = mirror_struct(ifo, 'ETM').MirrorMass
+    power = ifo_power(ifo)
+    noise_dict = noise.quantum.shotrad(f, ifo, mirror_mass, power)
+    pc.ASvac = noise_dict['ASvac']
+    pc.SEC = noise_dict['SEC']
+    pc.Arm = noise_dict['arm']
+    pc.Injection = noise_dict['injection']
+    pc.PD = noise_dict['pd']
+
+    # FC0 are the noises from the filter cavity losses and FC0_unsqzd_back
+    # are noises from the unsqueezed vacuum injected at the back mirror
+    # Right now there are at most one filter cavity in all the models;
+    # if there were two, there would also be FC1 and FC1_unsqzd_back, etc.
+    # keys = list(noise_dict.keys())
+    fc_keys = [key for key in noise_dict.keys() if 'FC' in key]
+    if fc_keys:
+        pc.FC = np.zeros_like(pc.ASvac)
+        for key in fc_keys:
+            pc.FC += noise_dict[key]
+
+    if 'phase' in noise_dict.keys():
+        pc.Phase = noise_dict['phase']
+
+    if 'ofc' in noise_dict.keys():
+        pc.OFC = noise_dict['OFC']
+
     return pc
 
 
-def precomp_power(f, ifo):
-    pc = Struct()
-    pbs, parm, finesse, prfactor, Tpr = ifo_power(ifo)
-    pc.pbs = pbs
-    pc.parm = parm
-    pc.finesse = finesse
-    pc.gPhase = finesse * 2/np.pi
-    pc.prfactor = prfactor
-    return pc
+############################################################
+# calibration
+############################################################
 
+def dhdl(f, armlen):
+    """Strain to length conversion for noise power spetra
 
-def precomp_cavity(f, ifo):
-    pc = Struct()
-    w0, wBeam_ITM, wBeam_ETM = arm_cavity(ifo)
-    pc.w0 = w0
-    pc.wBeam_ITM = wBeam_ITM
-    pc.wBeam_ETM = wBeam_ETM
-    return pc
+    This takes into account the GW wavelength and is only important
+    when this is comparable to the detector arm length.
+
+    From R. Schilling, CQG 14 (1997) 1513-1519, equation 5,
+    with n = 1, nu = 0.05, ignoring overall phase and cos(nu)^2.
+    A small value of nu is used instead of zero to avoid infinities.
+
+    Returns the square of the dh/dL function, and the same divided by
+    the arm length squared.
+
+    """
+    c = const.c
+    nu_small = 15*pi/180
+    omega_arm = pi * f * armlen / c
+    omega_arm_f = (1 - sin(nu_small)) * pi * f * armlen / c
+    omega_arm_b = (1 + sin(nu_small)) * pi * f * armlen / c
+    sinc_sqr = 4 / abs(sin(omega_arm_f) * exp(-1j * omega_arm) / omega_arm_f
+                       + sin(omega_arm_b) * exp(1j * omega_arm) / omega_arm_b)**2
+    dhdl_sqr = sinc_sqr / armlen**2
+    return dhdl_sqr, sinc_sqr
 
-##################################################
 
 class Strain(nb.Calibration):
     def calc(self):
         dhdl_sqr, sinc_sqr = dhdl(self.freq, self.ifo.Infrastructure.Length)
         return dhdl_sqr
 
-##################################################
+
+############################################################
+# noise sources
+############################################################
+
+#########################
+# quantum
+#########################
 
 class QuantumVacuum(nb.Noise):
     """Quantum Vacuum
@@ -197,10 +251,109 @@ class QuantumVacuum(nb.Noise):
         color='#ad03de',
     )
 
-    @nb.precomp(mirror=precomp_mirror)
-    @nb.precomp(power=precomp_power)
-    def calc(self, mirror, power):
-        return noise.quantum.shotrad(self.freq, self.ifo, mirror.ETM.MirrorMass, power)
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        total = np.zeros_like(quantum.ASvac)
+        for nn in quantum.values():
+            total += nn
+        return total
+
+
+class QuantumVacuumAS(nb.Noise):
+    """Quantum vacuum from the AS port
+
+    """
+    style = dict(
+        label='AS Port Vacuum',
+        color='xkcd:emerald green'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.ASvac
+
+
+class QuantumVacuumArm(nb.Noise):
+    """Quantum vacuum due to arm cavity loss
+
+    """
+    style = dict(
+        label='Arm Loss',
+        color='xkcd:orange brown'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.Arm
+
+
+class QuantumVacuumSEC(nb.Noise):
+    """Quantum vacuum due to SEC loss
+
+    """
+    style = dict(
+        label='SEC Loss',
+        color='xkcd:cerulean'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.SEC
+
+
+class QuantumVacuumFilterCavity(nb.Noise):
+    """Quantum vacuum due to filter cavity loss
+
+    """
+    style = dict(
+        label='Filter Cavity Loss',
+        color='xkcd:goldenrod'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.FC
+
+
+class QuantumVacuumInjection(nb.Noise):
+    """Quantum vacuum due to injection loss
+
+    """
+    style = dict(
+        label='Injection Loss',
+        color='xkcd:fuchsia'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.Injection
+
+
+class QuantumVacuumReadout(nb.Noise):
+    """Quantum vacuum due to readout loss
+
+    """
+    style = dict(
+        label='Readout Loss',
+        color='xkcd:mahogany'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.PD
+
+
+class QuantumVacuumQuadraturePhase(nb.Noise):
+    """Quantum vacuum noise due to quadrature phase noise
+    """
+    style = dict(
+        label='Quadrature Phase',
+        color='xkcd:slate'
+    )
+
+    @nb.precomp(quantum=precomp_quantum)
+    def calc(self, quantum):
+        return quantum.Phase
 
 
 class StandardQuantumLimit(nb.Noise):
@@ -213,10 +366,13 @@ class StandardQuantumLimit(nb.Noise):
         linestyle=":",
     )
 
-    @nb.precomp(mirror=precomp_mirror)
-    def calc(self, mirror):
-        return 8 * const.hbar / (mirror.ETM.MirrorMass * (2 * np.pi * self.freq) ** 2)
+    def calc(self):
+        ETM = mirror_struct(self.ifo, 'ETM')
+        return 8 * const.hbar / (ETM.MirrorMass * (2 * np.pi * self.freq) ** 2)
 
+#########################
+# seismic
+#########################
 
 class Seismic(nb.Noise):
     """Seismic
@@ -243,6 +399,10 @@ class Seismic(nb.Noise):
         return n * 4
 
 
+#########################
+# Newtonian
+#########################
+
 class Newtonian(nb.Noise):
     """Newtonian Gravity
 
@@ -300,6 +460,10 @@ class NewtonianInfrasound(nb.Noise):
         return n * 2
 
 
+#########################
+# suspension thermal
+#########################
+
 class SuspensionThermal(nb.Noise):
     """Suspension Thermal
 
@@ -316,6 +480,10 @@ class SuspensionThermal(nb.Noise):
         return n * 4
 
 
+#########################
+# coating thermal
+#########################
+
 class CoatingBrownian(nb.Noise):
     """Coating Brownian
 
@@ -325,16 +493,16 @@ class CoatingBrownian(nb.Noise):
         color='#fe0002',
     )
 
-    @nb.precomp(mirror=precomp_mirror)
-    @nb.precomp(cavity=precomp_cavity)
-    @nb.precomp(power=precomp_power)
-    def calc(self, mirror, cavity, power):
+    def calc(self):
+        ITM = mirror_struct(self.ifo, 'ITM')
+        ETM = mirror_struct(self.ifo, 'ETM')
+        cavity = arm_cavity(self.ifo)
         wavelength = self.ifo.Laser.Wavelength
         nITM = noise.coatingthermal.coating_brownian(
-            self.freq, mirror.ITM, wavelength, cavity.wBeam_ITM, power.parm,
+            self.freq, ITM, wavelength, cavity.wBeam_ITM
         )
         nETM = noise.coatingthermal.coating_brownian(
-            self.freq, mirror.ETM, wavelength, cavity.wBeam_ETM, power.parm,
+            self.freq, ETM, wavelength, cavity.wBeam_ETM
         )
         return (nITM + nETM) * 2
 
@@ -349,20 +517,25 @@ class CoatingThermoOptic(nb.Noise):
         linestyle='--',
     )
 
-    @nb.precomp(mirror=precomp_mirror)
-    @nb.precomp(cavity=precomp_cavity)
-    def calc(self, mirror, cavity):
+    def calc(self):
         wavelength = self.ifo.Laser.Wavelength
         materials = self.ifo.Materials
+        ITM = mirror_struct(self.ifo, 'ITM')
+        ETM = mirror_struct(self.ifo, 'ETM')
+        cavity = arm_cavity(self.ifo)
         nITM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic(
-            self.freq, mirror.ITM, wavelength, cavity.wBeam_ITM,
+            self.freq, ITM, wavelength, cavity.wBeam_ITM,
         )
         nETM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic(
-            self.freq, mirror.ETM, wavelength, cavity.wBeam_ETM,
+            self.freq, ETM, wavelength, cavity.wBeam_ETM,
         )
         return (nITM + nETM) * 2
 
 
+#########################
+# substrate thermal
+#########################
+
 class ITMThermoRefractive(nb.Noise):
     """ITM Thermo-Refractive
 
@@ -373,12 +546,13 @@ class ITMThermoRefractive(nb.Noise):
         linestyle='--',
     )
 
-    @nb.precomp(cavity=precomp_cavity)
-    @nb.precomp(power=precomp_power)
-    def calc(self, cavity, power):
+    def calc(self):
+        power = ifo_power(self.ifo)
+        gPhase = power.finesse * 2/np.pi
+        cavity = arm_cavity(self.ifo)
         n = noise.substratethermal.substrate_thermorefractive(
             self.freq, self.ifo.Materials, cavity.wBeam_ITM)
-        return n * 2 / power.gPhase**2
+        return n * 2 / gPhase**2
 
 
 class SubstrateBrownian(nb.Noise):
@@ -391,8 +565,8 @@ class SubstrateBrownian(nb.Noise):
         linestyle='--',
     )
 
-    @nb.precomp(cavity=precomp_cavity)
-    def calc(self, cavity):
+    def calc(self):
+        cavity = arm_cavity(self.ifo)
         nITM = noise.substratethermal.substrate_brownian(
             self.freq, self.ifo.Materials, cavity.wBeam_ITM)
         nETM = noise.substratethermal.substrate_brownian(
@@ -410,8 +584,8 @@ class SubstrateThermoElastic(nb.Noise):
         linestyle='--',
     )
 
-    @nb.precomp(cavity=precomp_cavity)
-    def calc(self, cavity):
+    def calc(self):
+        cavity = arm_cavity(self.ifo)
         nITM = noise.substratethermal.substrate_thermoelastic(
             self.freq, self.ifo.Materials, cavity.wBeam_ITM)
         nETM = noise.substratethermal.substrate_thermoelastic(
@@ -419,6 +593,11 @@ class SubstrateThermoElastic(nb.Noise):
         return (nITM + nETM) * 2
 
 
+#########################
+# residual gas
+#########################
+
+
 class ExcessGas(nb.Noise):
     """Excess Gas
 

gwinc/noise/coatingthermal.py

@@ -10,7 +10,7 @@ from ..const import BESSEL_ZEROS as zeta
 from ..const import J0M as j0m
 
 
-def coating_brownian(f, mirror, wavelength, wBeam, power):
+def coating_brownian(f, mirror, wavelength, wBeam, power=None):
     """Coating brownian noise for a given collection of coating layers
 
     This function calculates Coating Brownian noise using
@@ -32,22 +32,24 @@ def coating_brownian(f, mirror, wavelength, wBeam, power):
          wBeam = beam radius (at 1 / e**2 power)
          power = laser power falling on the mirror (W)
 
-    If the mirror.Coating.IncCoatBrAmpNoise parameter is present and
-    evaluates to True the amplitude noise due to coating brownian
-    noise will be calculated and its effect on the phase noise will be
-    added.  In that case the following new optional arguments should
-    be made available in the Materials object to provide separate Bulk
-    and Shear loss angles and to observe effect of photoelasticity:*
-     lossBlown = Coating Bulk Loss Angle of Low Refrac.Index layer @ 100Hz
-     lossSlown = Coating Shear Loss Angle of Low Refrac. Index layer @ 100Hz
+    If the power argument is present and is not None, the amplitude noise due
+    to coating brownian noise will be calculated and its effect on the phase
+    noise will be added (assuming the susceptibility is that of a free mass)
+
+    ***The following parameters are experimental and unsupported as yet***
+    The following optional parameters are available in the Materials object
+    to provide separate Bulk and Shear loss angles and to include the effect
+    of photoelasticity:
+    lossBlown = Coating Bulk Loss Angle of Low Refrac.Index layer @ 100Hz
+    lossSlown = Coating Shear Loss Angle of Low Refrac. Index layer @ 100Hz
     lossBhighn = Coating Bulk Loss Angle of High Refrac. Index layer @ 100Hz
     lossShighn = Coating Shear Loss Angle of High Refrac. Index layer @ 100Hz
-  lossBlown_slope = Coating Bulk Loss Angle Slope of Low Refrac. Index layer
-  lossSlown_slope = Coating Shear Loss Angle Slope of Low Refrac. Index layer
-  lossBhighn_slope = Coating Bulk Loss Angle Slope of High Refrac. Index layer
-  lossShighn_slope = Coating Shear Loss Angle Slope of High Refrac. Index layer
-       PETlown = Relevant component of Photoelastic Tensor of High n layer*
-      PEThighn = Relevant component of Photoelastic Tensor of Low n layer*
+    lossBlown_slope = Coating Bulk Loss Angle Slope of Low Refrac. Index layer
+    lossSlown_slope = Coating Shear Loss Angle Slope of Low Refrac. Index layer
+    lossBhighn_slope = Coating Bulk Loss Angle Slope of High Refrac. Index layer
+    lossShighn_slope = Coating Shear Loss Angle Slope of High Refrac. Index layer
+    PETlown = Relevant component of Photoelastic Tensor of High n layer*
+    PEThighn = Relevant component of Photoelastic Tensor of Low n layer*
 
     Returns:
       SbrZ = Brownian noise spectra for one mirror in m**2 / Hz
@@ -216,7 +218,7 @@ def coating_brownian(f, mirror, wavelength, wBeam, power):
 
     # From Sec II.E. Eq.(41)
     # Conversion of brownian amplitude noise to displacement noise
-    if coat.get('IncCoatBrAmpNoise'):
+    if power is not None:
 
         # get/calculate optic transmittance
         mTi = mirror.get('Transmittance', 1-np.abs(rho)**2)

gwinc/squeeze.py

 from numpy import pi, sqrt
 
 from . import const
+from .ifo.noises import ifo_power
 
 
 def sql(ifo):
     """Computer standard quantum limit (SQL) for IFO"""
     c = const.c
-    Parm = ifo.gwinc.parm
+    power = ifo_power(ifo)
     w0 = 2 * pi * c / ifo.Laser.Wavelength
-    m = ifo.Materials.MirrorMass
+    rho = ifo.Materials.Substrate.MassDensity
+    m = pi * ifo.Materials.MassRadius**2 * ifo.Materials.MassThickness * rho
     Titm = ifo.Optics.ITM.Transmittance
     Tsrm = ifo.Optics.SRM.Transmittance
     tSR = sqrt(Tsrm)
     rSR = sqrt(1 - Tsrm)
-    fSQL = (1/(2*pi))*(8/c)*sqrt((Parm*w0)/(m*Titm))*(tSR/(1+rSR))
+    fSQL = (1/(2*pi))*(8/c)*sqrt((power.parm*w0)/(m*Titm))*(tSR/(1+rSR))
     return fSQL
 
 
-def computeFCParams(ifo, fcParams):
+def computeFCParams(ifo):
     """Compute ideal filter cavity Tin, detuning [Hz] and bandwidth [Hz]
 
     """
     # FC parameters
+    fcParams = ifo.Squeezer.FilterCavity
     c = const.c
     fsrFC = c / (2 * fcParams.L)
     lossFC = fcParams.Lrt + fcParams.Te
@@ -40,7 +43,8 @@ def computeFCParams(ifo, fcParams):
     # input mirror transmission
     TinFC = 4 * pi * gammaFC / fsrFC - lossFC
     if TinFC < lossFC:
-        raise RuntimeError('IFC: Losses are too high! %.1f ppm max.' % 1e6 * gammaFC / fsrFC)
+        raise RuntimeError(
+            'IFC: Losses are too high! {:0.1f} ppm max.'.format(1e6 * gammaFC / fsrFC))
 
     # Add to fcParams structure
     fcParams.Ti = TinFC