diff --git a/README.md b/README.md
index bb04aa78fa18890846e377ab7bf34db91e3ff6a9..af9fdbee01bf17541ae87fe9be4690cb75a02570 100644
--- a/README.md
+++ b/README.md
@@ -17,8 +17,7 @@ description is loaded, the noise budget can be calculated and plotted:
>>> import numpy as np
>>> freq = np.logspace(1, 3, 1000)
>>> Budget = gwinc.load_budget('aLIGO')
->>> ifo = gwinc.precompIFO(freq, Budget.ifo)
->>> traces = Budget(freq, ifo=ifo).run()
+>>> traces = Budget(freq).run()
>>> fig = gwinc.plot_noise(freq, traces)
>>> fig.show()
```
diff --git a/gwinc/__init__.py b/gwinc/__init__.py
index 3347d489b628bd714fda268127e032d11db8031b..efa72a1cd02b15dd4279a56261d7668d79a3e054 100644
--- a/gwinc/__init__.py
+++ b/gwinc/__init__.py
@@ -6,7 +6,6 @@ import numpy as np
from .ifo import available_ifos
from .struct import Struct
-from .precomp import precompIFO
from .plot import plot_noise
from .io import load_hdf5, save_hdf5
@@ -112,7 +111,6 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
# FIXME: how do we allow adding arbitrary addtional noise sources
# from just ifo description, without having to specify full budget
Budget = load_budget('aLIGO')
- ifo = precompIFO(freq, ifo, PRfixed)
traces = Budget(freq, ifo=ifo).run()
plot_style = getattr(Budget, 'plot_style', {})
diff --git a/gwinc/__main__.py b/gwinc/__main__.py
index d6b220f38d919c0f1f76e288dbf1c11ca3b746a2..beb810e52a1ddbd781f3e9e50e77dfefb95ff5f7 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 . import available_ifos, load_budget, plot_noise
-from .precomp import precompIFO
from . import io
##################################################
@@ -175,9 +174,6 @@ def main():
# main calculations
if not traces:
- if ifo:
- logging.info("precomputing ifo...")
- ifo = precompIFO(freq, ifo)
logging.info("calculating budget...")
traces = Budget(freq=freq, ifo=ifo).run()
diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py
index 5c0b901596b6049e057773a8cf6918ab7e946abb..56e06b1adc577508966f2d3db03f3a08675432f4 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,49 @@ 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
+
+ if ifo.Suspension.Type == 'BQuad':
+ material = 'Silicon'
+ else:
+ material = 'Silica'
+ hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension, material)
+ 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 +179,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 +194,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 +276,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 +293,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 +316,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(
@@ -208,9 +337,11 @@ class ITMThermoRefractive(nb.Noise):
)
def calc(self):
- 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)**2
@@ -225,9 +356,11 @@ class ITMCarrierDensity(nb.Noise):
)
def calc(self):
- 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)**2
@@ -242,10 +375,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
@@ -260,10 +394,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
diff --git a/gwinc/test/__main__.py b/gwinc/test/__main__.py
index cfac045c4f110d111aa6af225b197f8b076f0446..41c684f492244b56ab83b9143f2c36fe9cca3b3b 100644
--- a/gwinc/test/__main__.py
+++ b/gwinc/test/__main__.py
@@ -15,7 +15,6 @@ logging.basicConfig(format='%(message)s',
level=os.getenv('LOG_LEVEL', logging.INFO))
from .. import available_ifos, load_budget
-from .. import precompIFO
from .. import load_hdf5, save_hdf5
try:
@@ -180,8 +179,7 @@ def main():
freq = np.logspace(np.log10(5), np.log10(6000), 3000)
for name in available_ifos():
Budget = load_budget(name)
- ifo = precompIFO(freq, Budget.ifo)
- traces = Budget(freq, ifo=ifo).run()
+ traces = Budget(freq).run()
path = os.path.join(args.cache, name+'.h5')
save_hdf5(path, freq, traces)
return
@@ -219,8 +217,7 @@ def main():
freq, tracesA, attrs = load_hdf5(path)
Budget = load_budget(name)
- ifo = precompIFO(freq, Budget.ifo)
- tracesB = Budget(freq, ifo=ifo).run()
+ tracesB = Budget(freq).run()
if inspiral_range:
totalA = tracesA['Total'][0]