diff --git a/README.md b/README.md index ca2693bf5d5442c35846a4971d835f23cd5abcd0..3f261c18eb862cadd48d225eb57bdc85a26143ee 100644 --- a/README.md +++ b/README.md @@ -113,13 +113,22 @@ For custom plotting, parameter optimization, etc. all functionality can be accessed directly through the `gwinc` library interface: ```python >>> import gwinc ->>> import numpy as np ->>> freq = np.logspace(1, 3, 1000) ->>> budget = gwinc.load_budget('aLIGO', freq) +>>> budget = gwinc.load_budget('aLIGO') >>> trace = budget.run() >>> fig = gwinc.plot_budget(trace) >>> fig.show() ``` +A default frequency array is used, but alternative frequencies can be +provided to `load_budget()` either in the form of a numpy array: +```python +>>> import numpy as np +>>> freq = np.logspace(1, 3, 1000) +>>> budget = gwinc.load_budget('aLIGO', freq=freq) +``` +or frequency specification string ('FLO:[NPOINTS:]FHI'): +``` +>>> budget = gwinc.load_budget('aLIGO', freq='10:1000:1000') +``` The `load_budget()` function takes most of the same inputs as the command line interface (e.g. IFO names, budget module paths, YAML @@ -134,6 +143,15 @@ properties. The budget sub-traces are available through a dictionary (`trace['QuantumVacuum']`) interface and via attributes (`trace.QuantumVacumm`). +The budget `freq` and `ifo` attributes can be updated at run time by +passing them as keyword arguments to the `run()` method: +```python +>>> budget = load_budget('aLIGO') +>>> freq = np.logspace(1, 3, 1000) +>>> ifo = Struct.from_file('/path/to/ifo_alt.yaml') +>>> trace = budget.run(freq=freq, ifo=ifo) +``` + ## noise functions @@ -340,12 +358,6 @@ If a budget module defined as a package includes an `ifo.yaml` [parameter file](#parameter-files) in the package directory, the `load_budget()` function will automatically load the YAML data into an `ifo` `gwinc.Struct` and assign it to the `budget.ifo` attribute. -Alternate ifos can be specified at run time: -```python -budget = load_budget('/path/to/MyBudget', freq) -ifo = Struct.from_file('/path/to/MyBudget.ifo') -trace = budget.run(ifo=ifo) -``` The IFOs included in `gwinc.ifo` provide examples of the use of the budget interface (e.g. [gwinc.ifo.aLIGO](gwinc/ifo/aLIGO)). diff --git a/gwinc/__init__.py b/gwinc/__init__.py index 134bc3740319cc201aafaf5a2c1d53e6d7355fc3..d3baf384962c0657fdd561ff0266330195433341 100644 --- a/gwinc/__init__.py +++ b/gwinc/__init__.py @@ -14,7 +14,32 @@ from .plot import plot_noise logger = logging.getLogger('gwinc') -def _load_module(name_or_path): +DEFAULT_FREQ = '5:3000:6000' + + +def freq_from_spec(spec=None): + """logarithmicly spaced frequency array, based on specification string + + Specification string should be of form 'START:[NPOINTS:]STOP'. If + `spec` is an array, then the array is returned as-is, and if it's + None the DEFAULT_FREQ specification is used. + + """ + if isinstance(spec, np.ndarray): + return spec + elif spec is None: + spec = DEFAULT_FREQ + fspec = spec.split(':') + if len(fspec) == 2: + fspec = fspec[0], DEFAULT_FREQ.split(':')[1], fspec[1] + return np.logspace( + np.log10(float(fspec[0])), + np.log10(float(fspec[2])), + int(fspec[1]), + ) + + +def load_module(name_or_path): """Load module from name or path. Return loaded module and module path. @@ -43,19 +68,24 @@ def load_budget(name_or_path, freq=None): Accepts either the name of a built-in canonical budget (see gwinc.IFOS), the path to a budget package (directory) or module - (ending in .py), or the path to an IFO struct definition file (see + (ending in .py), or the path to an IFO Struct definition file (see gwinc.Struct). If the budget is a package directory which includes an 'ifo.yaml' file the ifo Struct will be loaded from that file and assigned to - the budget.ifo attribute. If a struct definition file is provided - the base aLIGO budget definition and ifo will be assumed. + the budget.ifo attribute. If a Struct definition file is provided + the base aLIGO budget definition will be assumed. + + Returns an instantiated Budget object. If a frequency array or + frequency specification string (see `freq_from_spec()`) is + provided, the budget will be instantiated with the provided array. + If a frequency array is not provided and the Budget class + definition includes a `freq` attribute defining either an array or + specification string, then that array will be used. Otherwise a + default array will be provided (see DEFAULT_FREQ). - Returns a Budget object instantiated with the provided frequency - array (if specified), and with any included ifo assigned as an - object attribute. If a frequency array is not provided and the - budget class does not define it's own, the frequency array must be - provided at budget update() or run() time. + Any included ifo will be assigned as an attribute to the returned + Budget object. """ ifo = None @@ -83,8 +113,11 @@ def load_budget(name_or_path, freq=None): modname = 'gwinc.ifo.'+name_or_path logger.info("loading module {}...".format(modname)) - mod, modpath = _load_module(modname) + mod, modpath = load_module(modname) Budget = getattr(mod, bname) + if freq is None: + freq = getattr(Budget, 'freq', None) + freq = freq_from_spec(freq) ifopath = os.path.join(modpath, 'ifo.yaml') if not ifo and os.path.exists(ifopath): ifo = Struct.from_file(ifopath) diff --git a/gwinc/__main__.py b/gwinc/__main__.py index e2ab196698075e130bc2b1cb62016a124e06dfcb..9f3d071b49d56708fb5f5a3a12b8dd930df5048d 100644 --- a/gwinc/__main__.py +++ b/gwinc/__main__.py @@ -3,9 +3,16 @@ import os import signal import logging import argparse -import numpy as np -from . import IFOS, load_budget, plot_budget, logger +from . import ( + IFOS, + DEFAULT_FREQ, + freq_from_spec, + load_budget, + plot_budget, + logger, +) +from . import io logger.setLevel(os.getenv('LOG_LEVEL', 'WARNING').upper()) formatter = logging.Formatter('%(message)s') @@ -48,12 +55,12 @@ e.g.: gwinc --fom range:m1=20,m2=20 ... See the inspiral_range package documentation for details. +NOTE: The range will be calculated with the supplied frequency array, +and may therefore be inaccurate if a restricted array is specified. """ IFO = 'aLIGO' -FREQ = '5:3000:6000' FOM = 'range:m1=1.4,m2=1.4' -DATA_SAVE_FORMATS = ['.hdf5', '.h5'] parser = argparse.ArgumentParser( prog='gwinc', @@ -61,7 +68,7 @@ parser = argparse.ArgumentParser( formatter_class=argparse.RawDescriptionHelpFormatter) parser.add_argument( '--freq', '-f', metavar='FLO:[NPOINTS:]FHI', - help="logarithmic frequency array specification in Hz [{}]".format(FREQ)) + help="logarithmic frequency array specification in Hz [{}]".format(DEFAULT_FREQ)) parser.add_argument( '--ifo', '-o', metavar='PARAM=VAL', #nargs='+', action='extend', @@ -97,17 +104,6 @@ parser.add_argument( help="IFO name or path") -def freq_from_spec(spec): - fspec = spec.split(':') - if len(fspec) == 2: - fspec = fspec[0], FREQ.split(':')[1], fspec[1] - return np.logspace( - np.log10(float(fspec[0])), - np.log10(float(fspec[2])), - int(fspec[1]), - ) - - def main(): signal.signal(signal.SIGINT, signal.SIG_DFL) @@ -116,7 +112,7 @@ def main(): ########## # initial arg processing - if os.path.splitext(os.path.basename(args.IFO))[1] in DATA_SAVE_FORMATS: + if os.path.splitext(os.path.basename(args.IFO))[1] in io.DATA_SAVE_FORMATS: if args.freq: parser.exit(2, "Frequency specification not allowed when loading traces from file.\n") if args.ifo: @@ -129,15 +125,12 @@ def main(): plot_style = trace.plot_style else: - budget = load_budget(args.IFO) + try: + freq = freq_from_spec(args.freq) + except IndexError: + parser.exit(2, "Improper frequency specification: {}\n".format(args.freq)) + budget = load_budget(args.IFO, freq=freq) ifo = budget.ifo - if args.freq: - try: - freq = freq_from_spec(args.freq) - except IndexError: - parser.exit(2, "Improper frequency specification: {}\n".format(args.freq)) - else: - freq = getattr(budget, 'freq', freq_from_spec(FREQ)) plot_style = getattr(budget, 'plot_style', {}) trace = None @@ -178,7 +171,7 @@ def main(): if args.save: args.plot = False for path in args.save: - if os.path.splitext(path)[1] in DATA_SAVE_FORMATS: + if os.path.splitext(path)[1] in io.DATA_SAVE_FORMATS: out_data_files.add(path) out_plot_files = set(args.save) - out_data_files @@ -303,12 +296,11 @@ In [.]: plt.savefig("foo.pdf") # save noise trace to HDF5 file if out_data_files: - from .io import save_hdf5 for path in out_data_files: logger.info("saving budget trace: {}".format(path)) - save_hdf5( - path=path, + io.save_hdf5( trace=trace, + path=path, ifo=ifo, plot_style=plot_style, ) diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py index 014c6a8d52d43d3876a883e51f7c8dff60c46921..bd8d37c88303f226309653254a151d17d1342cbf 100644 --- a/gwinc/ifo/noises.py +++ b/gwinc/ifo/noises.py @@ -14,15 +14,6 @@ from .. import suspension # helper functions ############################################################ -def coating_thickness(materials, optic): - if 'CoatLayerOpticalThickness' in optic: - return optic['CoatLayerOpticalThickness'] - T = optic.Transmittance - dL = optic.CoatingThicknessLown - dCap = optic.CoatingThicknessCap - return noise.coatingthermal.getCoatDopt(materials, T, dL, dCap=dCap) - - def mirror_struct(ifo, tm): """Create a "mirror" Struct for a LIGO core optic @@ -35,7 +26,13 @@ def mirror_struct(ifo, tm): # it would otherwise be stored by reference) mirror = copy.deepcopy(ifo.Materials) optic = ifo.Optics.get(tm) - mirror.Coating.dOpt = coating_thickness(ifo.Materials, optic) + if 'CoatLayerOpticalThickness' in optic: + mirror.Coating.dOpt = optic['CoatLayerOpticalThickness'] + else: + T = optic.Transmittance + dL = optic.CoatingThicknessLown + dCap = optic.CoatingThicknessCap + mirror.Coating.dOpt = noise.coatingthermal.getCoatDopt(mirror, T, dL, dCap=dCap) mirror.update(optic) mirror.MassVolume = pi * mirror.MassRadius**2 * mirror.MassThickness mirror.MirrorMass = mirror.MassVolume * mirror.Substrate.MassDensity @@ -66,7 +63,14 @@ def arm_cavity(ifo): z2 = L * g1 * (1 - g2) / gden # waist, input, output - return w0, w1, w2 + cavity = Struct() + cavity.w0 = w0 + cavity.wBeam_ITM = w1 + cavity.wBeam_ETM = w2 + cavity.zr = zr + cavity.zBeam_ITM = z1 + cavity.zBeam_ETM = z2 + return cavity def ifo_power(ifo, PRfixed=True): @@ -110,7 +114,14 @@ def ifo_power(ifo, PRfixed=True): if pbs > pbsl: logger.warning('P_BS exceeds BS Thermal limit!') - return pbs, parm, finesse, prfactor, Tpr + power = Struct() + power.pbs = pbs + power.parm = parm + power.finesse = finesse + power.gPhase = finesse * 2/np.pi + power.prfactor = prfactor + power.Tpr = Tpr + return power ############################################################ @@ -132,32 +143,7 @@ def precomp_suspension(f, ifo): return pc -def precomp_mirror(f, ifo): - pc = Struct() - pc.ITM = mirror_struct(ifo, 'ITM') - pc.ETM = mirror_struct(ifo, 'ETM') - return pc - - -def precomp_power(f, ifo): - pc = Struct() - pbs, parm, finesse, prfactor, Tpr = ifo_power(ifo) - pc.pbs = pbs - pc.parm = parm - pc.finesse = finesse - pc.gPhase = finesse * 2/np.pi - pc.prfactor = prfactor - return pc - - -def precomp_cavity(f, ifo): - pc = Struct() - w0, wBeam_ITM, wBeam_ETM = arm_cavity(ifo) - pc.w0 = w0 - pc.wBeam_ITM = wBeam_ITM - pc.wBeam_ETM = wBeam_ETM - return pc - +################################################## def precomp_quantum(f, ifo): pc = Struct() @@ -243,12 +229,10 @@ class QuantumVacuum(nb.Noise): color='#ad03de', ) - @nb.precomp(quantum=precomp_quantum) - def calc(self, quantum): - total = np.zeros_like(quantum.ASvac) - for nn in quantum.values(): - total += nn - return total + def calc(self): + ETM = mirror_struct(self.ifo, 'ETM') + power = ifo_power(self.ifo) + return noise.quantum.shotrad(self.freq, self.ifo, ETM.MirrorMass, power) class QuantumVacuumAS(nb.Noise): @@ -358,9 +342,9 @@ class StandardQuantumLimit(nb.Noise): linestyle=":", ) - @nb.precomp(mirror=precomp_mirror) - def calc(self, mirror): - return 8 * const.hbar / (mirror.ETM.MirrorMass * (2 * np.pi * self.freq) ** 2) + def calc(self): + ETM = mirror_struct(self.ifo, 'ETM') + return 8 * const.hbar / (ETM.MirrorMass * (2 * np.pi * self.freq) ** 2) ######################### @@ -486,16 +470,17 @@ class CoatingBrownian(nb.Noise): color='#fe0002', ) - @nb.precomp(mirror=precomp_mirror) - @nb.precomp(cavity=precomp_cavity) - @nb.precomp(power=precomp_power) - def calc(self, mirror, cavity, power): + def calc(self): + ITM = mirror_struct(self.ifo, 'ITM') + ETM = mirror_struct(self.ifo, 'ETM') + power = ifo_power(self.ifo) + cavity = arm_cavity(self.ifo) wavelength = self.ifo.Laser.Wavelength nITM = noise.coatingthermal.coating_brownian( - self.freq, mirror.ITM, wavelength, cavity.wBeam_ITM, power.parm, + self.freq, ITM, wavelength, cavity.wBeam_ITM, power.parm, ) nETM = noise.coatingthermal.coating_brownian( - self.freq, mirror.ETM, wavelength, cavity.wBeam_ETM, power.parm, + self.freq, ETM, wavelength, cavity.wBeam_ETM, power.parm, ) return (nITM + nETM) * 2 @@ -510,16 +495,17 @@ class CoatingThermoOptic(nb.Noise): linestyle='--', ) - @nb.precomp(mirror=precomp_mirror) - @nb.precomp(cavity=precomp_cavity) - def calc(self, mirror, cavity): + def calc(self): wavelength = self.ifo.Laser.Wavelength materials = self.ifo.Materials + ITM = mirror_struct(self.ifo, 'ITM') + ETM = mirror_struct(self.ifo, 'ETM') + cavity = arm_cavity(self.ifo) nITM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic( - self.freq, mirror.ITM, wavelength, cavity.wBeam_ITM, + self.freq, ITM, wavelength, cavity.wBeam_ITM, ) nETM, junk1, junk2, junk3 = noise.coatingthermal.coating_thermooptic( - self.freq, mirror.ETM, wavelength, cavity.wBeam_ETM, + self.freq, ETM, wavelength, cavity.wBeam_ETM, ) return (nITM + nETM) * 2 @@ -538,12 +524,13 @@ class ITMThermoRefractive(nb.Noise): linestyle='--', ) - @nb.precomp(cavity=precomp_cavity) - @nb.precomp(power=precomp_power) - def calc(self, cavity, power): + def calc(self): + power = ifo_power(self.ifo) + gPhase = power.finesse * 2/np.pi + cavity = arm_cavity(self.ifo) n = noise.substratethermal.substrate_thermorefractive( self.freq, self.ifo.Materials, cavity.wBeam_ITM) - return n * 2 / power.gPhase**2 + return n * 2 / gPhase**2 class SubstrateBrownian(nb.Noise): @@ -556,8 +543,8 @@ class SubstrateBrownian(nb.Noise): linestyle='--', ) - @nb.precomp(cavity=precomp_cavity) - def calc(self, cavity): + def calc(self): + cavity = arm_cavity(self.ifo) nITM = noise.substratethermal.substrate_brownian( self.freq, self.ifo.Materials, cavity.wBeam_ITM) nETM = noise.substratethermal.substrate_brownian( @@ -575,8 +562,8 @@ class SubstrateThermoElastic(nb.Noise): linestyle='--', ) - @nb.precomp(cavity=precomp_cavity) - def calc(self, cavity): + def calc(self): + cavity = arm_cavity(self.ifo) nITM = noise.substratethermal.substrate_thermoelastic( self.freq, self.ifo.Materials, cavity.wBeam_ITM) nETM = noise.substratethermal.substrate_thermoelastic( diff --git a/gwinc/io.py b/gwinc/io.py index 98b900e5bd02b4099d982fe9dc8f3b2aa1e2a65d..13d765007055fef027650b7fc8786e79dfd87201 100644 --- a/gwinc/io.py +++ b/gwinc/io.py @@ -8,6 +8,7 @@ from .trace import BudgetTrace SCHEMA = 'GWINC noise budget' +DATA_SAVE_FORMATS = ['.hdf5', '.h5'] def _write_trace_recursive(f, trace): @@ -20,16 +21,17 @@ def _write_trace_recursive(f, trace): _write_trace_recursive(tgrp, t) -def save_hdf5(path, trace, ifo=None, plot_style=None, **kwargs): +def save_hdf5(trace, path, ifo=None, plot_style=None, **kwargs): """Save GWINC budget data to an HDF5 file. - The `freq` argument should be the frequency array, and `traces` - should be the traces (recursive) dictionary. Keyword arguments - are stored in the HDF5 top level 'attrs' key-value store. If an - 'ifo' keyword arg is supplied, it is assumed to be a Struct and - will be serialized to YAML for storage. + `trace` should be a BudgetTrace object (see budget.run()), and + `path` should be the path to the output HDF5 file. Keyword + arguments are stored in the HDF5 top level 'attrs' key-value + store. If an 'ifo' keyword arg is supplied, it is assumed to be a + Struct and will be serialized to YAML for storage. - See HDF5_SCHEMA. + See HDF5_SCHEMA for information on the BudgetTrace/HDF5 storage + format. """ with h5py.File(path, 'w') as f: @@ -87,27 +89,28 @@ def _load_hdf5_v1(f): def _load_hdf5_v2(f): - def read_recursive(element): + def read_recursive(element, freq): psd = element['PSD'][:] style = yaml.safe_load(element.attrs['style']) budget = [] if 'budget' in element: for name, item in element['budget'].items(): - trace = read_recursive(item) + trace = read_recursive(item, freq) trace.name = name budget.append(trace) return BudgetTrace( + freq=freq, psd=psd, budget=budget, style=style, ) + freq = f['Freq'][:] for name, item in f.items(): if name == 'Freq': - freq = item[:] + continue else: - trace = read_recursive(item) + trace = read_recursive(item, freq) trace.name = name - trace._freq = freq trace.ifo = None attrs = dict(f.attrs) ifo = attrs.get('ifo') @@ -117,18 +120,19 @@ def _load_hdf5_v2(f): del attrs['ifo'] except yaml.constructor.ConstructorError: logger.warning("HDF5 load warning: Could not de-serialize 'ifo' YAML attribute.") - trace.plot_style = yaml.safe_load(attrs['plot_style']) + trace.plot_style = yaml.safe_load(attrs.get('plot_style', '')) return trace def load_hdf5(path): """Load GWINC budget data from an HDF5 file. - Returns BudgetTrace. An 'ifo' attr will be de-serialized from - YAML into a Struct object and assigned as an attribute to the - BudgetTrace. + Returns BudgetTrace for trace data stored in the HDF5 file at + `path`. An 'ifo' attr will be de-serialized from YAML into a + Struct object and assigned as an attribute to the BudgetTrace. - See HDF5_SCHEMA. + See HDF5_SCHEMA for information on the BudgetTrace/HDF5 storage + format. """ loaders = { diff --git a/gwinc/nb.py b/gwinc/nb.py index 54391b881c460c34bed8cb4bbbaf7fa6f8dda6ad..9069798058dde1d11599caded16488626d4dc3be 100644 --- a/gwinc/nb.py +++ b/gwinc/nb.py @@ -151,6 +151,7 @@ class BudgetItem: if freq is not None: assert isinstance(freq, np.ndarray) self.freq = freq + self.ifo = None for key, val in kwargs.items(): setattr(self, key, val) self._loaded = False @@ -551,7 +552,7 @@ class Budget(Noise): _precomp=_precomp, ) budget.append(trace) - total = quadsum([trace.psd for trace in budget]) + total = quadsum([trace.psd for trace in budget if trace.name in self._budget_noises]) return self._make_trace( psd=total, budget=budget ) diff --git a/gwinc/plot.py b/gwinc/plot.py index 1753108d36102db91ee83d897b03853754e27fa3..884dbd57d8c3446b30e5744fca8af7d8044c8fe0 100644 --- a/gwinc/plot.py +++ b/gwinc/plot.py @@ -34,7 +34,7 @@ def plot_budget( for name, trace in budget.items(): style = trace.style if 'label' not in style: - style['label'] = budget.name + style['label'] = name if 'linewidth' in style: style['lw'] = style['linewidth'] elif 'lw' not in style: diff --git a/gwinc/squeeze.py b/gwinc/squeeze.py index 5cbd9dcbfc4d66f47e630ee364bd7c4f91a2359a..1dc304f85e39919436cfd9a86e719119b7f6aebf 100644 --- a/gwinc/squeeze.py +++ b/gwinc/squeeze.py @@ -1,27 +1,30 @@ from numpy import pi, sqrt from . import const +from .ifo.noises import ifo_power def sql(ifo): """Computer standard quantum limit (SQL) for IFO""" c = const.c - Parm = ifo.gwinc.parm + power = ifo_power(ifo) w0 = 2 * pi * c / ifo.Laser.Wavelength - m = ifo.Materials.MirrorMass + rho = ifo.Materials.Substrate.MassDensity + m = pi * ifo.Materials.MassRadius**2 * ifo.Materials.MassThickness * rho Titm = ifo.Optics.ITM.Transmittance Tsrm = ifo.Optics.SRM.Transmittance tSR = sqrt(Tsrm) rSR = sqrt(1 - Tsrm) - fSQL = (1/(2*pi))*(8/c)*sqrt((Parm*w0)/(m*Titm))*(tSR/(1+rSR)) + fSQL = (1/(2*pi))*(8/c)*sqrt((power.parm*w0)/(m*Titm))*(tSR/(1+rSR)) return fSQL -def computeFCParams(ifo, fcParams): +def computeFCParams(ifo): """Compute ideal filter cavity Tin, detuning [Hz] and bandwidth [Hz] """ # FC parameters + fcParams = ifo.Squeezer.FilterCavity c = const.c fsrFC = c / (2 * fcParams.L) lossFC = fcParams.Lrt + fcParams.Te @@ -40,7 +43,8 @@ def computeFCParams(ifo, fcParams): # input mirror transmission TinFC = 4 * pi * gammaFC / fsrFC - lossFC if TinFC < lossFC: - raise RuntimeError('IFC: Losses are too high! %.1f ppm max.' % 1e6 * gammaFC / fsrFC) + raise RuntimeError( + 'IFC: Losses are too high! {:0.1f} ppm max.'.format(1e6 * gammaFC / fsrFC)) # Add to fcParams structure fcParams.Ti = TinFC diff --git a/gwinc/struct.py b/gwinc/struct.py index 2de5efef366cee6ad04ac379eea92500c9db0a36..f6aac55ff7603b3d2381fbc343813a8a2cc74a0d 100644 --- a/gwinc/struct.py +++ b/gwinc/struct.py @@ -376,7 +376,7 @@ class Struct(object): base, ext = os.path.splitext(path) if ext == '.m': - from ..gwinc_matlab import Matlab + from .gwinc_matlab import Matlab matlab = Matlab() matlab.addpath(os.path.dirname(path)) func_name = os.path.basename(base)