remove precomp

This finally gets rid of the whole precompIFO function.  It does this
by breaking up precomp into various smaller functions, moving them
into the common IFO noise definition module, and using them
appropriately when needed.  This does result in some of the functions
being called multiple times in the full budget calculation for the
aLIGO-like budgets, but the cost should be minor given the convenience
of getting rid of the whole precomp thing, thereby opening up
everything to non-aLIGO-like configurations.

Even though some functions are calculated multiple times, we're still
orders of magnitude faster than matgwinc, so...  Should still revisit
down the line though.

closes #40

README.md

@@ -17,7 +17,6 @@ description is loaded, the noise budget can be calculated and plotted:
 >>> import numpy as np
 >>> freq = np.logspace(1, 3, 1000)
 >>> Budget, ifo, freq_, plot_style = gwinc.load_ifo('aLIGO')
->>> ifo = gwinc.precompIFO(freq, ifo)
 >>> traces = Budget(freq, ifo=ifo).calc_trace()
 >>> fig = gwinc.plot_noise(freq, traces, **plot_style)
 >>> fig.show()

gwinc/__init__.py

 from .ifo import available_ifos, load_ifo
-from .precomp import precompIFO
 from .struct import Struct
 from .plot import plot_noise
 from .io import load_hdf5, save_hdf5

gwinc/__main__.py

@@ -10,7 +10,6 @@ logging.basicConfig(format='%(message)s',
                     level=os.getenv('LOG_LEVEL', logging.INFO))
 
 from .ifo import available_ifos, load_ifo
-from .precomp import precompIFO
 from . import plot_noise
 from . import io
 
@@ -174,9 +173,6 @@ def main():
     # main calculations
 
     if not traces:
-        if ifo:
-            logging.info("precomputing ifo...")
-            ifo = precompIFO(freq, ifo)
         logging.info("calculating budget...")
         budget = Budget(freq=freq, ifo=ifo)
         budget.load()

gwinc/gwinc.py

@@ -5,12 +5,11 @@ from collections import OrderedDict
 import logging
 
 from . import load_ifo
-from .precomp import precompIFO
 from . import noise
 from .plot import plot_noise
 
 
-def gwinc(freq, ifoin, source=None, plot=False, PRfixed=True):
+def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
     """Calculate strain noise budget for a specified interferometer model.
 
     Argument `freq` is the frequency array for which the noises will
@@ -31,10 +30,6 @@ def gwinc(freq, ifoin, source=None, plot=False, PRfixed=True):
     # from just ifo description, without having to specify full budget
     Budget, ifo_, freq_, plot_style = load_ifo('aLIGO')
 
-    # add some precomputed info to the ifo struct
-    # this implicitly deepcopies and the return value is the copy
-    ifo = precompIFO(freq, ifoin, PRfixed)
-
     traces = Budget(freq, ifo=ifo).calc_trace()
 
     # construct matgwinc-compatible noises structure

gwinc/ifo/noises.py

@@ -2,11 +2,95 @@ import numpy as np
 from numpy import pi, sin, exp, sqrt
 
 from .. import const
+from ..struct import Struct
 from .. import nb
 from .. import noise
+from .. import suspension
 
 ##################################################
 
+def coating_thickness(ifo, optic):
+    optic = ifo.Optics.get(optic)
+    if 'CoatLayerOpticalThickness' in optic:
+        return optic['CoatLayerOpticalThickness']
+    T = optic.Transmittance
+    dL = optic.CoatingThicknessLown
+    dCap = optic.CoatingThicknessCap
+    return noise.coatingthermal.getCoatDopt(ifo.Materials, T, dL, dCap=dCap)
+
+
+def arm_cavity(ifo):
+    L = ifo.Infrastructure.Length
+
+    g1 = 1 - L / ifo.Optics.Curvature.ITM
+    g2 = 1 - L / ifo.Optics.Curvature.ETM
+    gcav = sqrt(g1 * g2 * (1 - g1 * g2))
+    gden = g1 - 2 * g1 * g2 + g2
+
+    if (g1 * g2 * (1 - g1 * g2)) <= 0:
+        raise Exception('Unstable arm cavity g-factors.  Change ifo.Optics.Curvature')
+    elif gcav < 1e-3:
+        logging.warning('Nearly unstable arm cavity g-factors.  Reconsider ifo.Optics.Curvature')
+
+    ws = sqrt(L * ifo.Laser.Wavelength / pi)
+    w1 = ws * sqrt(abs(g2) / gcav)
+    w2 = ws * sqrt(abs(g1) / gcav)
+
+    # waist size
+    w0 = ws * sqrt(gcav / abs(gden))
+    zr = pi * w0**2 / ifo.Laser.Wavelength
+    z1 = L * g2 * (1 - g1) / gden
+    z2 = L * g1 * (1 - g2) / gden
+
+    # waist, input, output
+    return w0, w1, w2
+
+
+def ifo_power(ifo, PRfixed=True):
+    """Compute power on beamsplitter, finesse, and power recycling factor.
+
+    """
+    c = const.c
+    pin = ifo.Laser.Power
+    t1 = sqrt(ifo.Optics.ITM.Transmittance)
+    r1 = sqrt(1 - ifo.Optics.ITM.Transmittance)
+    r2 = sqrt(1 - ifo.Optics.ETM.Transmittance)
+    t5 = sqrt(ifo.Optics.PRM.Transmittance)
+    r5 = sqrt(1 - ifo.Optics.PRM.Transmittance)
+    loss = ifo.Optics.Loss  # single TM loss
+    bsloss = ifo.Optics.BSLoss
+    acoat = ifo.Optics.ITM.CoatingAbsorption
+    pcrit = ifo.Optics.pcrit
+
+    # Finesse, effective number of bounces in cavity, power recycling factor
+    finesse = 2*pi / (t1**2 + 2*loss)  # arm cavity finesse
+    neff    = 2 * finesse / pi
+
+    # Arm cavity reflectivity with finite loss
+    garm = t1 / (1 - r1*r2*sqrt(1-2*loss))  # amplitude gain wrt input field
+    rarm = r1 - t1 * r2 * sqrt(1-2*loss) * garm
+
+    if PRfixed:
+        Tpr = ifo.Optics.PRM.Transmittance  # use given value
+    else:
+        Tpr = 1-(rarm*sqrt(1-bsloss))**2  # optimal recycling mirror transmission
+        t5 = sqrt(Tpr)
+        r5 = sqrt(1 - Tpr)
+    prfactor = t5**2 / (1 + r5 * rarm * sqrt(1-bsloss))**2
+
+    pbs  = pin * prfactor  # BS power from input power
+    parm = pbs * garm**2 / 2  # arm power from BS power
+
+    thickness = ifo.Optics.ITM.get('Thickness', ifo.Materials.MassThickness)
+    asub = 1.3 * 2 * thickness * ifo.Optics.SubstrateAbsorption
+    pbsl = 2 *pcrit / (asub+acoat*neff)  # bs power limited by thermal lensing
+
+    if pbs > pbsl:
+        logging.warning('P_BS exceeds BS Thermal limit!')
+
+    return pbs, parm, finesse, prfactor, Tpr
+
+
 def dhdl(f, armlen):
     """Strain to length conversion for noise power spetra
 
@@ -35,6 +119,44 @@ def dhdl(f, armlen):
 
 ##################################################
 
+def precomp_mirror(f, ifo):
+    ifo.Materials.MirrorVolume = \
+        pi * ifo.Materials.MassRadius**2 \
+        * ifo.Materials.MassThickness
+    ifo.Materials.MirrorMass = \
+        ifo.Materials.MirrorVolume \
+        * ifo.Materials.Substrate.MassDensity
+
+def precomp_gwinc(f, ifo):
+    ifo.gwinc = Struct()
+    pbs, parm, finesse, prfactor, Tpr = ifo_power(ifo)
+    ifo.gwinc.pbs = pbs
+    ifo.gwinc.parm = parm
+    ifo.gwinc.finesse = finesse
+    ifo.gwinc.prfactor = prfactor
+
+
+def precomp_suspension(f, ifo):
+    if 'VHCoupling' not in ifo.Suspension:
+        ifo.Suspension.VHCoupling = Struct()
+        ifo.Suspension.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
+
+    hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension)
+    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
+    ifo.Suspension.hTable = hTable
+    ifo.Suspension.vTable = vTable
+
+##################################################
+
 class Strain(nb.Calibration):
     def calc(self):
         dhdl_sqr, sinc_sqr = dhdl(self.freq, self.ifo.Infrastructure.Length)
@@ -52,6 +174,8 @@ class QuantumVacuum(nb.Noise):
     )
 
     def calc(self):
+        precomp_mirror(self.freq, self.ifo)
+        precomp_gwinc(self.freq, self.ifo)
         return noise.quantum.shotrad(self.freq, self.ifo)
 
 
@@ -65,6 +189,7 @@ class Seismic(nb.Noise):
     )
 
     def calc(self):
+        precomp_suspension(self.freq, self.ifo)
         if 'PlatformMotion' in self.ifo.Seismic:
             if self.ifo.Seismic.PlatformMotion == 'BSC':
                 nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
@@ -146,6 +271,7 @@ class SuspensionThermal(nb.Noise):
     )
 
     def calc(self):
+        precomp_suspension(self.freq, self.ifo)
         n = noise.suspensionthermal.suspension_thermal(self.freq, self.ifo.Suspension)
         return n * 4
 
@@ -162,10 +288,9 @@ class CoatingBrownian(nb.Noise):
     def calc(self):
         wavelength = self.ifo.Laser.Wavelength
         materials = self.ifo.Materials
-        wBeam_ITM = self.ifo.Optics.ITM.BeamRadius
-        wBeam_ETM = self.ifo.Optics.ETM.BeamRadius
-        dOpt_ITM = self.ifo.Optics.ITM.CoatLayerOpticalThickness
-        dOpt_ETM = self.ifo.Optics.ETM.CoatLayerOpticalThickness
+        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
+        dOpt_ITM = coating_thickness(self.ifo, 'ITM')
+        dOpt_ETM = coating_thickness(self.ifo, 'ETM')
         nITM = noise.coatingthermal.coating_brownian(
             self.freq, materials, wavelength, wBeam_ITM, dOpt_ITM)
         nETM = noise.coatingthermal.coating_brownian(
@@ -186,10 +311,9 @@ class CoatingThermoOptic(nb.Noise):
     def calc(self):
         wavelength = self.ifo.Laser.Wavelength
         materials = self.ifo.Materials
-        wBeam_ITM = self.ifo.Optics.ITM.BeamRadius
-        wBeam_ETM = self.ifo.Optics.ETM.BeamRadius
-        dOpt_ITM = self.ifo.Optics.ITM.CoatLayerOpticalThickness
-        dOpt_ETM = self.ifo.Optics.ETM.CoatLayerOpticalThickness
+        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
+        dOpt_ITM = coating_thickness(self.ifo, 'ITM')
+        dOpt_ETM = coating_thickness(self.ifo, 'ETM')
         nITM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic(
             self.freq, materials, wavelength, wBeam_ITM, dOpt_ITM[:])
         nETM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic(
@@ -209,9 +333,11 @@ class ITMThermoRefractive(nb.Noise):
 
     def calc(self):
         L = self.ifo.Infrastructure.Length
-        gPhase = self.ifo.gwinc.finesse * 2/np.pi
+        finesse = ifo_power(self.ifo)[2]
+        gPhase =  finesse * 2/np.pi
+        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
         n = noise.substratethermal.substrate_thermorefractive(
-            self.freq, self.ifo.Materials, self.ifo.Optics.ITM.BeamRadius)
+            self.freq, self.ifo.Materials, wBeam_ITM)
         return n * 2 / (gPhase*L)**2
 
 
@@ -227,9 +353,11 @@ class ITMCarrierDensity(nb.Noise):
 
     def calc(self):
         L = self.ifo.Infrastructure.Length
-        gPhase = self.ifo.gwinc.finesse * 2/np.pi
+        finesse = ifo_power(self.ifo)[2]
+        gPhase = finesse * 2/np.pi
+        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
         n = noise.substratethermal.substrate_carrierdensity(
-            self.freq, self.ifo.Materials, self.ifo.Optics.ITM.BeamRadius)
+            self.freq, self.ifo.Materials, wBeam_ITM)
         return n * 2 / (gPhase*L)**2
 
 
@@ -244,10 +372,11 @@ class SubstrateBrownian(nb.Noise):
     )
 
     def calc(self):
+        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
         nITM = noise.substratethermal.substrate_brownian(
-            self.freq, self.ifo.Materials, self.ifo.Optics.ITM.BeamRadius)
+            self.freq, self.ifo.Materials, wBeam_ITM)
         nETM = noise.substratethermal.substrate_brownian(
-            self.freq, self.ifo.Materials, self.ifo.Optics.ETM.BeamRadius)
+            self.freq, self.ifo.Materials, wBeam_ETM)
         return (nITM + nETM) * 2
 
 
@@ -262,10 +391,11 @@ class SubstrateThermoElastic(nb.Noise):
     )
 
     def calc(self):
+        w0, wBeam_ITM, wBeam_ETM = arm_cavity(self.ifo)
         nITM = noise.substratethermal.substrate_thermoelastic(
-            self.freq, self.ifo.Materials, self.ifo.Optics.ITM.BeamRadius)
+            self.freq, self.ifo.Materials, wBeam_ITM)
         nETM = noise.substratethermal.substrate_thermoelastic(
-            self.freq, self.ifo.Materials, self.ifo.Optics.ETM.BeamRadius)
+            self.freq, self.ifo.Materials, wBeam_ETM)
         return (nITM + nETM) * 2
 
 

gwinc/precomp.py

-from __future__ import division
-from numpy import pi, sqrt, sin, exp
-from scipy.io import loadmat
-import logging
-import copy
-
-from . import const
-from . import suspension
-from .struct import Struct
-from .noise.coatingthermal import getCoatDopt
-
-
-def precompIFO(f, ifoin, PRfixed=True):
-    """Add precomputed data to the IFO model.
-
-    To prevent recomputation of these precomputed data, if the
-    ifo argument contains ifo.gwinc.PRfixed, and this matches
-    the argument PRfixed, no changes are made.
-
-    This function DOES NOT modify IFO. Its return value is a copy of
-    IFO with precomputations filled out.
-
-    """
-    ifo = copy.deepcopy(ifoin)
-
-    if 'gwinc' not in ifo:
-        ifo.gwinc = Struct()
-
-    ifo.gwinc.PRfixed = PRfixed
-
-    ##############################
-    # derived temp
-
-    if 'Temp' not in ifo.Materials.Substrate:
-        ifo.Materials.Substrate.Temp = ifo.Infrastructure.Temp
-
-    ##############################
-    # suspension vertical-horizontal coupling
-
-    if 'VHCoupling' not in ifo.Suspension:
-        ifo.Suspension.VHCoupling = Struct()
-        ifo.Suspension.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
-
-    ##############################
-    # optics values
-
-    # calculate optics' parameters
-    ifo.Materials.MirrorVolume = pi*ifo.Materials.MassRadius**2 * \
-                                 ifo.Materials.MassThickness
-    ifo.Materials.MirrorMass = ifo.Materials.MirrorVolume * \
-                               ifo.Materials.Substrate.MassDensity
-    ifo.Optics.ITM.Thickness = ifo.Materials.MassThickness
-
-    # coating layer optical thicknesses - mevans 2 May 2008
-    if 'CoatLayerOpticalThickness' not in ifo.Optics.ITM:
-        T = ifo.Optics.ITM.Transmittance
-        dL = ifo.Optics.ITM.CoatingThicknessLown
-        dCap = ifo.Optics.ITM.CoatingThicknessCap
-        ifo.Optics.ITM.CoatLayerOpticalThickness = getCoatDopt(ifo.Materials, T, dL, dCap=dCap)
-        T = ifo.Optics.ETM.Transmittance
-        dL = ifo.Optics.ETM.CoatingThicknessLown
-        dCap = ifo.Optics.ETM.CoatingThicknessCap
-        ifo.Optics.ETM.CoatLayerOpticalThickness = getCoatDopt(ifo.Materials, T, dL, dCap=dCap)
-
-    ##############################
-    # beam parameters
-
-    armlen = ifo.Infrastructure.Length
-
-    # g-factors
-    g1 = 1 - armlen / ifo.Optics.Curvature.ITM
-    g2 = 1 - armlen / ifo.Optics.Curvature.ETM
-    gcav = sqrt(g1 * g2 * (1 - g1 * g2))
-    gden = g1 - 2 * g1 * g2 + g2
-
-    if (g1 * g2 * (1 - g1 * g2)) <= 0:
-        raise Exception('Unstable arm cavity g-factors.  Change ifo.Optics.Curvature')
-    elif gcav < 1e-3:
-        logging.warning('Nearly unstable arm cavity g-factors.  Reconsider ifo.Optics.Curvature')
-
-    ws = sqrt(armlen * ifo.Laser.Wavelength / pi)
-    w1 = ws * sqrt(abs(g2) / gcav)
-    w2 = ws * sqrt(abs(g1) / gcav)
-
-    # waist size
-    w0 = ws * sqrt(gcav / abs(gden))
-    zr = pi * w0**2 / ifo.Laser.Wavelength
-    z1 = armlen * g2 * (1 - g1) / gden
-    z2 = armlen * g1 * (1 - g2) / gden
-
-    ifo.Optics.ITM.BeamRadius = w1
-    ifo.Optics.ETM.BeamRadius = w2
-
-    ##############################
-    # calc power and IFO parameters
-
-    pbs, parm, finesse, prfactor, Tpr = precompPower(ifo, PRfixed)
-
-    ifo.gwinc.pbs = pbs
-    ifo.gwinc.parm = parm
-    ifo.gwinc.finesse = finesse
-    ifo.gwinc.prfactor = prfactor
-    ifo.gwinc.gITM = g1
-    ifo.gwinc.gETM = g2
-    ifo.gwinc.BeamWaist = w0
-    ifo.gwinc.BeamRayleighRange = zr
-    ifo.gwinc.BeamWaistToITM = z1
-    ifo.gwinc.BeamWaistToETM = z2
-    ifo.Optics.PRM.Transmittance = Tpr
-
-    ##############################
-    # saved seismic spectrum
-
-    if 'darmSeiSusFile' in ifo.Seismic and ifo.Seismic.darmSeiSusFile:
-        darmsei = loadmat(ifo.Seismic.darmSeiSusFile)
-        ifo.Seismic.darmseis_f = darmsei['darmseis_f'][0]
-        ifo.Seismic.darmseis_x = darmsei['darmseis_x'][0]
-
-    # --------------------------------------------------------
-    # Suspension
-    # if the suspension code supports different temps for the stages
-    fname = eval('suspension.susp{}'.format(ifo.Suspension.Type))
-    hForce, vForce, hTable, vTable = fname(f, ifo.Suspension)
-
-    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
-
-    ifo.Suspension.hTable = hTable
-    ifo.Suspension.vTable = vTable
-
-    return ifo
-
-
-def precompPower(ifo, PRfixed=True):
-    """Compute power on beamsplitter, finesse, and power recycling factor.
-
-    """
-    c       = const.c
-    pin     = ifo.Laser.Power
-    t1      = sqrt(ifo.Optics.ITM.Transmittance)
-    r1      = sqrt(1 - ifo.Optics.ITM.Transmittance)
-    r2      = sqrt(1 - ifo.Optics.ETM.Transmittance)
-    t5      = sqrt(ifo.Optics.PRM.Transmittance)
-    r5      = sqrt(1 - ifo.Optics.PRM.Transmittance)
-    loss    = ifo.Optics.Loss                          # single TM loss
-    bsloss  = ifo.Optics.BSLoss
-    acoat   = ifo.Optics.ITM.CoatingAbsorption
-    pcrit   = ifo.Optics.pcrit
-
-    # Finesse, effective number of bounces in cavity, power recycling factor
-    finesse = 2*pi / (t1**2 + 2*loss)        # arm cavity finesse
-    neff    = 2 * finesse / pi
-
-    # Arm cavity reflectivity with finite loss
-    garm = t1 / (1 - r1*r2*sqrt(1-2*loss))  # amplitude gain wrt input field
-    rarm = r1 - t1 * r2 * sqrt(1-2*loss) * garm
-
-    if PRfixed:
-        Tpr = ifo.Optics.PRM.Transmittance  # use given value
-    else:
-        Tpr = 1-(rarm*sqrt(1-bsloss))**2 # optimal recycling mirror transmission
-        t5 = sqrt(Tpr)
-        r5 = sqrt(1 - Tpr)
-    prfactor = t5**2 / (1 + r5 * rarm * sqrt(1-bsloss))**2
-
-    pbs  = pin * prfactor          # BS power from input power
-    parm = pbs * garm**2 / 2       # arm power from BS power
-
-    asub = 1.3*2*ifo.Optics.ITM.Thickness*ifo.Optics.SubstrateAbsorption
-    pbsl = 2*pcrit/(asub+acoat*neff) # bs power limited by thermal lensing
-
-    if pbs > pbsl:
-        logging.warning('P_BS exceeds BS Thermal limit!')
-
-    return pbs, parm, finesse, prfactor, Tpr