From 334bf5d7a54387c28780ff33b11c888a6f8ae189 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 | 1 -
gwinc/__init__.py | 1 -
gwinc/__main__.py | 4 -
gwinc/gwinc.py | 7 +-
gwinc/ifo/noises.py | 159 ++++++++++++++++++++++++++++++++++----
gwinc/precomp.py | 183 --------------------------------------------
6 files changed, 145 insertions(+), 210 deletions(-)
delete mode 100644 gwinc/precomp.py
diff --git a/README.md b/README.md
index 0c293704..fb66b0c2 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 9511e40e..71e5bf03 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 b20919da..29c25869 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 93b46d95..e4fe6890 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 9dc2212c..11d23029 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
@@ -228,8 +354,9 @@ class ITMCarrierDensity(nb.Noise):
def calc(self):
L = self.ifo.Infrastructure.Length
gPhase = self.ifo.gwinc.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 +371,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 +390,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
--
GitLab