From c96bdbd6828752898ed22ec04083d92331c39dc1 Mon Sep 17 00:00:00 2001
From: Jameson Graef Rollins <jrollins@finestructure.net>
Date: Fri, 31 Jan 2020 11:38:10 -0800
Subject: [PATCH] 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 | 3 +-
gwinc/__init__.py | 2 -
gwinc/__main__.py | 4 -
gwinc/ifo/noises.py | 167 +++++++++++++++++++++++++++++++++----
gwinc/precomp.py | 183 -----------------------------------------
gwinc/test/__main__.py | 7 +-
6 files changed, 154 insertions(+), 212 deletions(-)
delete mode 100644 gwinc/precomp.py
diff --git a/README.md b/README.md
index bb04aa78..af9fdbee 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 3347d489..efa72a1c 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 d6b220f3..beb810e5 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 5c0b9015..56e06b1a 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 554834c2..00000000
--- 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 cfac045c..41c684f4 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]
--
GitLab