diff --git a/README.md b/README.md
index 0c29370458d1d68583963655f6b4b8157fac1a96..fb66b0c24b32eb9d9b20e90328989cf2e9e15998 100644
--- a/README.md
+++ b/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()
diff --git a/gwinc/__init__.py b/gwinc/__init__.py
index 9511e40e3f0ffe39e688e5d8f5932220fdd8b214..71e5bf03ecade9be17c7c5078708c2bd2a2cd5c4 100644
--- a/gwinc/__init__.py
+++ b/gwinc/__init__.py
@@ -1,5 +1,4 @@
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
diff --git a/gwinc/__main__.py b/gwinc/__main__.py
index b20919da799a681221cac83e851e2ac67b4aedc0..29c2586909109373a647d06812a16fdf47822036 100644
--- a/gwinc/__main__.py
+++ b/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()
diff --git a/gwinc/gwinc.py b/gwinc/gwinc.py
index 93b46d958a744bc85904f6e6f7ae62cdf9239152..e4fe68902960d74b61e7aee4e34d1209115a1f79 100644
--- a/gwinc/gwinc.py
+++ b/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
diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py
index 9dc2212c40d7d4bf6640f6f9ad941d244d05617d..b55f5650f162852b32ec8d97a008f7f198999967 100644
--- a/gwinc/ifo/noises.py
+++ b/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
diff --git a/gwinc/precomp.py b/gwinc/precomp.py
deleted file mode 100644
index 554834c23e038eeb05e8b397a8daeb40b650687b..0000000000000000000000000000000000000000
--- a/gwinc/precomp.py
+++ /dev/null
@@ -1,183 +0,0 @@
-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