diff --git a/HDF5_SCHEMATA b/HDF5_SCHEMATA
deleted file mode 100644
index f7989feee01cdef64904d7bdae023299de3dad75..0000000000000000000000000000000000000000
--- a/HDF5_SCHEMATA
+++ /dev/null
@@ -1,107 +0,0 @@
-# HDF5 Schema for GWINC noise trace storage
-
-
-This file describes a schemata for HDF5 storage of noise trace data
-and plot styling GWINC noise budgets. Python functions for writing
-budget data to and reading budget data from this format are included
-in the `gwinc.io` module.
-
-HDF5 is a heirarchical, structured data storage format [0]. Content
-is organized into a heirarchical folder-like structure, with two
-types of named objects:
-
- * groups: holder of other objects (like folders)
- * datasets: holder of data arrays (like files)
-
-Objects can also have attributes as (almost) arbitrary key:value
-pairs.
-
-Bindings are available for most platforms including Python [1] and
-Matlab.
-
-[0] https://en.wikipedia.org/wiki/Hierarchical_Data_Format
-[1] http://www.h5py.org/
-
-
-## version history
-
-v1
-- first versioned schema release
-
-
-## schema
-
-The following describes the noise budget schema. Specific strings are
-enclosed in single quotes (''), and variables are described in
-brackets (<>). Group objects are indicated by a closing '/'
-separator, data set are indicated by a closing ':' followed by a
-specification of their length and type (e.g. "(N),float"), and
-attributes are specified in the .attrs[] dictionary format. Optional
-elements are enclosed in parentheses.
-
-A single trace is a length N array (with optional plot style specified
-in attributes:
-```
-/<trace>: (N),float
-(/<trace>.attrs['label'] = <label>)
-(/<trace>.attrs['color] = <color>)
-...
-```
-
-A budget item, i.e. a collection of noises is structured as follows:
-```
-/<budget>/
- /'Total': (N),float
- /<trace_0>: (N),float
- (/<trace_1>: (N),float)
-```
-<!-- ``` -->
-<!-- /'noises'/ -->
-<!-- /*<noise_0>: (N),float -->
-<!-- ... -->
-<!-- /\*'references'/ -->
-<!-- /*<ref_0>: (N),float -->
-<!-- ... -->
-<!-- ``` -->
-
-
-## Top-level Budget
-
-The following two root attributes expected: a string describing the schema,
-and an int schema version number:
-```
-/.attrs['SCHEMA'] = 'GWINC Noise Budget'
-/.attrs['SCHEMA_VERSION'] = 1
-```
-
-Top-level attributes can be used for storing any relevant metadata,
-such as "ifo" parameters in YAML format or overall plot styling.
-Typically at least the following root attributes are defined:
-```
-/.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
-/.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>
-```
-
-The budget frequency array is defined in a 'Freq' dataset:
-```
-/'Freq': (N), float
-```
-
-The budget traces are defined a traces group. The overall structure
-looks something like this:
-```
-/.attrs['SCHEMA'] = 'GWINC Noise Budget'
-/.attrs['SCHEMA_VERSION'] = 1
-/.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
-/.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>
-/'Freq': (N), float
-/traces/
- /'Total': (N),float
- /<noise_0>: (N),float
- /<noise_1>: (N),float
- /<noise_2>/
- /'Total': (N),float
- /<noise_3>: (N),float
- /<noise_4>: (N),float
- ...
-```
diff --git a/HDF5_SCHEMATA.md b/HDF5_SCHEMATA.md
new file mode 100644
index 0000000000000000000000000000000000000000..107980f3c8ec0aa0c4f96daf047dca2bc732fe95
--- /dev/null
+++ b/HDF5_SCHEMATA.md
@@ -0,0 +1,122 @@
+# HDF5 Schema for GWINC noise trace storage
+
+This file describes a schemata for HDF5 storage of noise trace data
+and plot styling GWINC noise budgets. Python functions for writing
+budget data to and reading budget data from this format are included
+in the `gwinc.io` module.
+
+HDF5 is a heirarchical, structured data storage format [0]. Content
+is organized into a heirarchical folder-like structure, with two
+types of named objects:
+
+ * groups: holder of other objects (like folders)
+ * datasets: holder of data arrays (like files)
+
+Objects can also have attributes as (almost) arbitrary key:value
+pairs.
+
+Bindings are available for most platforms including Python [1] and
+Matlab.
+
+ * [0] https://en.wikipedia.org/wiki/Hierarchical_Data_Format
+ * [1] http://www.h5py.org/
+
+
+## schema
+
+The following describes the noise budget schema. Specific strings are
+enclosed in single quotes (''), and variables are described in
+brackets (<>). Group objects are indicated by a closing '/', data
+sets are indicated by a closing ':' followed by a specification of
+their length and type (e.g. "(N),float"), and attributes are specified
+in the .attrs[] dictionary format. Optional elements are enclosed in
+parentheses.
+
+
+## top-level attributes
+
+The following two root attributes expected: a string describing the schema,
+and an int schema version number:
+```
+/.attrs['SCHEMA'] = 'GWINC Noise Budget'
+/.attrs['SCHEMA_VERSION'] = 1
+```
+
+The following root root attributes are defined:
+```
+/.attrs['title'] = <experiment description string (e.g. 'H1 Strain Budget')>
+/.attrs['date'] = <ISO-formatted string (e.g. '2015-10-24T20:30:00.000000Z')>
+/.attrs['ifo'] = <IFO Struct object, YAML serialized>
+```
+
+The remaining root level attributes usually pertain to the plot style.
+
+The budget frequency array is defined in a top level 'Freq' dataset:
+
+```
+/'Freq': (N),float
+```
+
+## version history
+
+### v1
+
+A single trace is a length N array (with optional plot style specified
+in attributes:
+```
+/<trace>: (N),float
+/<trace>.attrs['label'] = <label>
+/<trace>.attrs['color] = <color>
+...
+```
+
+A budget item, i.e. a collection of noises is structured as follows:
+```
+/<budget>/
+/<budget>/Total': (N),float
+/<budget>/<trace_0>: (N),float
+/<budget>/...
+```
+
+
+The budget traces are defined a traces group. The overall structure
+looks something like this:
+```
+/traces/
+/traces/'Total': (N),float
+/traces/<noise_0>: (N),float
+/traces/<noise_1>: (N),float
+/traces/<noise_2>/
+/traces/<noise_2>/'Total': (N),float
+/traces/<noise_2>/<noise_3>: (N),float
+/traces/<noise_2>/...
+/traces/...
+```
+
+
+### v2
+
+Each trace is given the following structure:
+```
+/<trace>/
+/<trace>.attrs['style'] = <YAML trace style dict>
+/<trace>/'PSD': (N),float
+/<trace>/budget/
+/<trace>/budget/<subtrace_0>/
+/<trace>/budget/<subtrace_1>/
+/<trace>/budget/...
+```
+
+The overall structure is:
+```
+/<budget>/
+/<budget>/.attrs['plot_style'] = <YAML plot style dict>
+/<budget>/.attrs['style'] = <YAML "Total" trace style dict>
+/<budget>/.attrs[*] = <arbitrary data>
+/<budget>/'Freq': (N),float
+/<budget>/'PSD': (N),float
+/<budget>/budget/
+/<budget>/budget/<trace_0>/...
+/<budget>/budget/<trace_1>/...
+/<budget>/budget/...
+```
diff --git a/README.md b/README.md
index c1869f1f1d92d035fac176f475e81d0e7af5f2e2..ca2693bf5d5442c35846a4971d835f23cd5abcd0 100644
--- a/README.md
+++ b/README.md
@@ -39,8 +39,8 @@ data.
The [`inspiral_range`](https://git.ligo.org/gwinc/inspiral-range)
package can be used to calculate various common "inspiral range"
-figures of merit for gravitational wave detector budgets. See
-[figures of merit](#figures-of-merit) section below.
+figures of merit for gravitational wave detector budgets. See the
+[inspiral range](#inspiral-range) section below.
## usage
@@ -48,12 +48,12 @@ figures of merit for gravitational wave detector budgets. See
### command line interface
`pygwinc` provides a command line interface that can be used to
-calculate and plot noise budgets for generic noise budgets or the
-various canonical IFOs described above, save/plot hdf5 trace data, and
-dump budget IFO parameters:
+calculate and plot the various canonical IFO noise budgets described
+above, as well as custom noise budgets (see below):
```shell
$ python3 -m gwinc aLIGO
```
+Budget plots save/plot hdf5 trace data, and dump budget IFO parameters:
You can play with IFO parameters and see the effects on the budget by
dumping the pre-defined parameters to a [YAML-formatted parameter
@@ -74,7 +74,7 @@ command line:
$ python3 -m gwinc aLIGO --ifo Optics.SRM.Tunephase=3.14
```
-Stand-alone YAML files will always assume the nominal ['aLIGO' budget
+Stand-alone YAML files assume the nominal ['aLIGO' budget
description](gwinc/ifo/aLIGO).
[Custom budgets](#budget-interface) may also be processed by providing
@@ -84,18 +84,19 @@ $ python3 -m gwinc path/to/mybudget
```
The command line interface also includes an "interactive" mode which
-provides an [IPython](https://ipython.org/) shell for interacting with a processed budget:
+provides an [IPython](https://ipython.org/) shell for interacting with
+a processed budget:
```shell
$ python3 -m gwinc -i Aplus
GWINC interactive shell
-The 'ifo' Struct and 'traces' data are available for inspection.
+The 'ifo' Struct and 'trace' data are available for inspection.
Use the 'whos' command to view the workspace.
You may interact with the plot using the 'plt' functions, e.g.:
-In [.]: plt.title("foo")
-In [.]: plt.savefig("foo.pdf")
+In [.]: plt.title("My Special Budget")
+In [.]: plt.savefig("mybudget.pdf")
In [1]:
```
@@ -106,7 +107,7 @@ $ python3 -m gwinc -h
```
-### python library
+### library interface
For custom plotting, parameter optimization, etc. all functionality can be
accessed directly through the `gwinc` library interface:
@@ -114,31 +115,32 @@ accessed directly through the `gwinc` library interface:
>>> import gwinc
>>> import numpy as np
>>> freq = np.logspace(1, 3, 1000)
->>> Budget = gwinc.load_budget('aLIGO')
->>> traces = Budget(freq).run()
->>> fig = gwinc.plot_noise(freq, traces)
+>>> budget = gwinc.load_budget('aLIGO', freq)
+>>> trace = budget.run()
+>>> fig = gwinc.plot_budget(trace)
>>> fig.show()
```
The `load_budget()` function takes most of the same inputs as the
command line interface (e.g. IFO names, budget module paths, YAML
-parameter files), and returns the un-instantiated `Budget` class
-defined in the specified budget module (see [budget
-interface](#budget-interface) below).
+parameter files), and returns the instantiated `Budget` object defined
+in the specified budget module (see [budget
+interface](#budget-interface) below). The budget `ifo` `gwinc.Struct`
+is available in the `budget.ifo` attribute.
-The budget `run()` method is a convenience method that calculates all
-budget noises and the noise total and returns a (possibly) nested
-dictionary of a noise data, in the form of a `(data, style)` tuple
-where 'data' is the PSD data and 'style' is a plot style dictionary
-for the trace. The dictionary will be nested if the budget includes
-any sub-budgets.
+The budget `run()` method calculates all budget noises and the noise
+total and returns a `BudgetTrace` object with `freq`, `psd`, and `asd`
+properties. The budget sub-traces are available through a dictionary
+(`trace['QuantumVacuum']`) interface and via attributes
+(`trace.QuantumVacumm`).
## noise functions
-`pygwinc` noise functions are available in the `gwinc.noise` package.
-This package includes multiple sub-modules for the different types of
-noises, e.g. `suspensionthermal`, `coatingthermal`, `quantum`, etc.)
+The `pygwinc` analytical noise functions are available in the
+`gwinc.noise` package. This package includes multiple sub-modules for
+the different types of noises, e.g. `suspensionthermal`,
+`coatingthermal`, `quantum`, etc.)
The various noise functions need many different parameters to
calculate their noise outputs. Many parameters are expected to be in
@@ -161,7 +163,7 @@ def coating_brownian(f, materials, wavelength, wBeam, dOpt):
```
-### `gwinc.Struct` objects
+## `gwinc.Struct` objects
`pygwinc` provides a `Struct` class that can hold parameters in
attributes and additionally acts like a dictionary, for passing to the
@@ -232,8 +234,12 @@ are:
* `update(**kwargs)`: update data/attributes
* `calc()`: return final data array
-See the built-in documentation for more info (e.g. `pydoc3
-gwinc.nb.BudgetItem`)
+Generally these methods are not called directly. Instead, the `Noise`
+and `Budget` classes include a `run` method that smartly executes them
+in sequence and returns a `BudgetTrace` object for the budget.
+
+See the built-in `BudgetItem` documentation for more info
+(e.g. `pydoc3 gwinc.nb.BudgetItem`)
### budget module definition
@@ -313,11 +319,10 @@ The `style` attributes of the various `Noise` classes define plot
style for the noise.
This budget can be loaded with the `gwinc.load_budget()` function, and
-processed with the `Budget.run()` method:
+processed as usual with the `Budget.run()` method:
```python
-Budget = load_budget('/path/to/MyBudget')
-budget = Budget(freq)
-traces = budget.run()
+budget = load_budget('/path/to/MyBudget', freq)
+trace = budget.run()
```
Other than the necessary `freq` initialization argument that defines
@@ -333,63 +338,71 @@ suspension Struct is extracted from the `self.ifo` Struct at
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 a
-`gwinc.Struct` and include it as an `Budget.ifo` attribute in the
-returned `Budget` class. This would provide the `self.ifo` needed in
-the `SuspensionThermal` Noise class above and is therefore a
-convenient way to provide parameter structures in budget packages.
-Otherwise it would need to be created/loaded in some other way and
-passed to the budget at instantiation, e.g.:
+`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')
+budget = load_budget('/path/to/MyBudget', freq)
ifo = Struct.from_file('/path/to/MyBudget.ifo')
-budget = Budget(freq, ifo=ifo)
-traces = budget.run()
+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)).
-### "precomp" functions
+### the "precomp" decorator
-`BudgetItems` support "precomp" functions that are calculated during
-the `update()` method call and can be used to cache common derived
-values in the `ifo` struct for later use. The execution state of
-these functions is cached at the Budget level, to prevent needlessly
-re-calculating them in multiple BudgetItems. When they are calculated
-they are provided with the `freq` and `ifo` attributes as arguments.
-For example, take the following budget definition:
+The `BudgetItem` supports "precomp" functions that can be used to
+calculate common derived values needed in multiple `BudgetItems`.
+They are specified using the `nb.precomp` decorator applied to the
+`BudgetItem.calc()` method. These functions are executed during the
+`update()` method call, supplied with the budget `freq` and `ifo`
+attributes as input arguments, and are expected to update the `ifo`
+struct. The execution state of the precomp functions is cached at the
+Budget level, to prevent needlessly re-calculating them multiple
+times. For example:
```python
+from gwinc import nb
-class Noise0(Noise):
- precomp = [
- precomp_foo,
- ]
+
+def precomp_foo(freq, ifo):
+ pc = Struct()
...
+ return pc
-class Noise1(Noise):
- precomp = [
- precomp_foo,
- precomp_bar,
- ]
+
+def precomp_bar(freq, ifo):
+ pc = Struct()
...
+ return pc
-class MyBudget(Budget):
+
+class Noise0(nb.Noise):
+ @nb.precomp(foo=precomp_foo)
+ def calc(self, foo):
+ ...
+
+class Noise1(nb.Noise):
+ @nb.precomp(foo=precomp_foo)
+ @nb.precomp(bar=precomp_bar)
+ def calc(self, foo, bar):
+ ...
+
+class MyBudget(nb.Budget):
noises = [
Noise0,
Noise1,
]
```
-When the `MyBudget.update()` method is called, all the `precomp`
-functions will be execute:
+When `MyBudget.run()` is called, all the `precomp` functions will be
+executed, e.g.:
```python
precomp_foo(self.freq, self.ifo)
precomp_bar(self.freq, self.ifo)
```
-But the `precomp_foo` function will only be calculated once per budget
-update() call, even though it's specified as needed by both `Noise0`
-and `Noise1`.
+But the `precomp_foo` function will only be calculated once even
+though it's specified as needed by both `Noise0` and `Noise1`.
### extracting single noise terms
@@ -399,10 +412,9 @@ interface. The most straightforward way is to run the full budget,
and extract the noise data the output traces dictionary:
```python
-Budget = load_budget('/path/to/MyBudget')
-budget = Budget(freq)
-traces = budget.calc_traces()
-data, plot_style = traces['QuantumVacuum']
+budget = load_budget('/path/to/MyBudget', freq)
+trace = budget.run()
+quantum_trace = trace['QuantumVacuum']
```
You can also calculate the final calibrated output noise for just a
@@ -420,7 +432,7 @@ data = budget['QuantumVacuum'].calc()
```
-# figures of merit
+# inspiral range
The [`inspiral_range`](https://git.ligo.org/gwinc/inspiral-range)
package can be used to calculate various common "inspiral range"
@@ -433,18 +445,15 @@ import inspiral_range
import numpy as np
freq = np.logspace(1, 3, 1000)
-Budget = gwinc.load_budget('Aplus')
-traces = Budget(freq).run()
+budget = gwinc.load_budget('Aplus', freq)
+trace = budget.run()
range = inspiral_range.range(
- freq, traces['Total'][0],
+ freq, trace.psd,
m1=30, m2=30,
)
```
-Note you need to extract the zeroth element of the `traces['Total']`
-tuple, which is the actual PSD data.
-
See the [`inspiral_range`](https://git.ligo.org/gwinc/inspiral-range)
package for more details.
diff --git a/gwinc/__init__.py b/gwinc/__init__.py
index 0c1592dae9a73522424c4640b8a575109639959b..134bc3740319cc201aafaf5a2c1d53e6d7355fc3 100644
--- a/gwinc/__init__.py
+++ b/gwinc/__init__.py
@@ -7,6 +7,7 @@ import numpy as np
from .ifo import IFOS
from .struct import Struct
+from .plot import plot_budget
from .plot import plot_noise
@@ -37,24 +38,24 @@ def _load_module(name_or_path):
return mod, path
-def load_budget(name_or_path):
- """Load GWINC IFO Budget by name or from path.
+def load_budget(name_or_path, freq=None):
+ """Load GWINC Budget
- Named IFOs should correspond to one of the IFOs available in the
- gwinc package (see gwinc.IFOS). If a path is provided it should
- either be a budget package (directory) or module (ending in .py),
- or an IFO struct definition (see gwinc.Struct).
+ 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
+ gwinc.Struct).
- If a budget package path is provided, and the package includes an
- 'ifo.yaml' file, that file will be loaded into a Struct and
- assigned as an attribute to the returned Budget class.
+ 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.
- If a bare struct is provided the base aLIGO budget definition will
- be assumed.
-
- Returns a Budget class with 'ifo', 'freq', and 'plot_style', ifo
- structure, frequency array, and plot style dictionary, with the
- last three being None if they are not defined in the budget.
+ 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.
"""
ifo = None
@@ -68,11 +69,9 @@ def load_budget(name_or_path):
ifo = Struct.from_file(path)
bname = 'aLIGO'
modname = 'gwinc.ifo.aLIGO'
- logger.info("loading budget {}...".format(modname))
else:
modname = path
- logger.info("loading module path {}...".format(modname))
else:
if name_or_path not in IFOS:
@@ -82,17 +81,14 @@ def load_budget(name_or_path):
))
bname = name_or_path
modname = 'gwinc.ifo.'+name_or_path
- logger.info("loading module {}...".format(modname))
+ logger.info("loading module {}...".format(modname))
mod, modpath = _load_module(modname)
-
Budget = getattr(mod, bname)
ifopath = os.path.join(modpath, 'ifo.yaml')
if not ifo and os.path.exists(ifopath):
ifo = Struct.from_file(ifopath)
- Budget.ifo = ifo
-
- return Budget
+ return Budget(freq=freq, ifo=ifo)
def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
@@ -114,15 +110,16 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
# assume generic aLIGO configuration
# FIXME: how do we allow adding arbitrary addtional noise sources
# from just ifo description, without having to specify full budget
- Budget = load_budget('aLIGO')
- traces = Budget(freq, ifo=ifo).run()
- plot_style = getattr(Budget, 'plot_style', {})
+ budget = load_budget('aLIGO', freq)
+ traces = budget.run()
+ plot_style = getattr(budget, 'plot_style', {})
# construct matgwinc-compatible noises structure
noises = {}
- for name, (data, style) in traces.items():
- noises[style.get('label', name)] = data
- noises['Freq'] = freq
+ for name, trace in traces.items():
+ noises[name] = trace.psd
+ noises['Total'] = traces.psd
+ noises['Freq'] = traces.freq
pbs = ifo.gwinc.pbs
parm = ifo.gwinc.parm
@@ -169,6 +166,6 @@ def gwinc(freq, ifo, source=None, plot=False, PRfixed=True):
logger.info('BBH Inspiral Range: ' + str(score.effr0bh) + ' Mpc/ z = ' + str(score.zHorizonBH))
logger.info('Stochastic Omega: %4.1g Universes' % score.Omega)
- plot_noise(ifo, traces, **plot_style)
+ plot_budget(traces, **plot_style)
return score, noises, ifo
diff --git a/gwinc/__main__.py b/gwinc/__main__.py
index 6b6142bedb1b5b58bc6edc81c596f1e6116cf1d3..e2ab196698075e130bc2b1cb62016a124e06dfcb 100644
--- a/gwinc/__main__.py
+++ b/gwinc/__main__.py
@@ -5,10 +5,10 @@ import logging
import argparse
import numpy as np
-from . import IFOS, load_budget, plot_noise, logger
+from . import IFOS, load_budget, plot_budget, logger
logger.setLevel(os.getenv('LOG_LEVEL', 'WARNING').upper())
-formatter = logging.Formatter('%(name)s: %(message)s')
+formatter = logging.Formatter('%(message)s')
handler = logging.StreamHandler()
handler.setFormatter(formatter)
logger.addHandler(handler)
@@ -61,7 +61,7 @@ parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument(
'--freq', '-f', metavar='FLO:[NPOINTS:]FHI',
- help="frequency array specification in Hz [{}]".format(FREQ))
+ help="logarithmic frequency array specification in Hz [{}]".format(FREQ))
parser.add_argument(
'--ifo', '-o', metavar='PARAM=VAL',
#nargs='+', action='extend',
@@ -71,7 +71,7 @@ parser.add_argument(
'--title', '-t',
help="plot title")
parser.add_argument(
- '--fom', metavar='FUNC[:PARAM=VAL,...]',
+ '--fom', metavar='FUNC[:PARAM=VAL,...]', default=FOM,
help="use inspiral_range.FUNC to calculate range figure-of-merit on resultant spectrum")
group = parser.add_mutually_exclusive_group()
group.add_argument(
@@ -118,29 +118,28 @@ def main():
if os.path.splitext(os.path.basename(args.IFO))[1] in DATA_SAVE_FORMATS:
if args.freq:
- parser.exit(2, "Can not specify frequency array when loading traces from file.\n")
+ parser.exit(2, "Frequency specification not allowed when loading traces from file.\n")
if args.ifo:
- parser.exit(2, "Can not override ifo parameters when loading traces from file.\n")
+ parser.exit(2, "IFO parameter specification not allowed when loading traces from file.\n")
from .io import load_hdf5
- Budget = None
- freq, traces, attrs = load_hdf5(args.IFO)
- ifo = attrs.get('ifo')
- # FIXME: deprecate 'IFO'
- ifo = attrs.get('IFO', ifo)
- plot_style = attrs
+ budget = None
+ trace = load_hdf5(args.IFO)
+ freq = trace.freq
+ ifo = trace.ifo
+ plot_style = trace.plot_style
else:
- Budget = load_budget(args.IFO)
- ifo = Budget.ifo
+ budget = load_budget(args.IFO)
+ ifo = budget.ifo
if args.freq:
try:
freq = freq_from_spec(args.freq)
except IndexError:
- parser.error("improper frequency specification '{}'".format(args.freq))
+ 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', {})
- traces = None
+ freq = getattr(budget, 'freq', freq_from_spec(FREQ))
+ plot_style = getattr(budget, 'plot_style', {})
+ trace = None
if args.ifo:
for paramval in args.ifo:
@@ -149,19 +148,19 @@ def main():
if args.yaml:
if not ifo:
- parser.exit(2, "No 'ifo' structure available.\n")
+ parser.exit(2, "IFO structure not provided.\n")
print(ifo.to_yaml(), end='')
return
if args.text:
if not ifo:
- parser.exit(2, "No 'ifo' structure available.\n")
+ parser.exit(2, "IFO structure not provided.\n")
print(ifo.to_txt(), end='')
return
if args.diff:
if not ifo:
- parser.exit(2, "No 'ifo' structure available.\n")
- Budget = load_budget(args.diff)
- diffs = ifo.diff(Budget.ifo)
+ parser.exit(2, "IFO structure not provided.\n")
+ dbudget = load_budget(args.diff)
+ diffs = ifo.diff(dbudget.ifo)
if diffs:
w = max([len(d[0]) for d in diffs])
fmt = '{{:{}}} {{:>20}} {{:>20}}'.format(w)
@@ -204,22 +203,21 @@ def main():
import inspiral_range
logger_ir = logging.getLogger('inspiral_range')
logger_ir.setLevel(logger.getEffectiveLevel())
- logger_ir.removeHandler(logger_ir.handlers[0])
- logger_ir.addHandler(logger.handlers[0])
- if not args.fom:
- args.fom = FOM
+ handler = logging.StreamHandler()
+ handler.setFormatter(logging.Formatter('%(name)s: %(message)s'))
+ logger_ir.addHandler(handler)
+
except ModuleNotFoundError:
- if args.fom:
- raise
- logger.warning("inspiral_range package not available, figure of merit will not be calculated...")
+ logger.warning("WARNING: inspiral_range package not available, figure of merit will not be calculated.")
+ args.fom = None
if args.fom:
try:
- range_func, fargs = args.fom.split(':')
+ range_func, range_func_args = args.fom.split(':')
except ValueError:
range_func = args.fom
- fargs = ''
+ range_func_args = ''
range_params = {}
- for param in fargs.split(','):
+ for param in range_func_args.split(','):
if not param:
continue
p, v = param.split('=')
@@ -234,37 +232,33 @@ def main():
##########
# main calculations
- if not traces:
+ if not trace:
logger.info("calculating budget...")
- traces = Budget(freq=freq, ifo=ifo).run()
-
- # logger.info('recycling factor: {: >0.3f}'.format(ifo.gwinc.prfactor))
- # logger.info('BS power: {: >0.3f} W'.format(ifo.gwinc.pbs))
- # logger.info('arm finesse: {: >0.3f}'.format(ifo.gwinc.finesse))
- # logger.info('arm power: {: >0.3f} kW'.format(ifo.gwinc.parm/1000))
+ trace = budget.run(freq=freq)
if args.title:
plot_style['title'] = args.title
elif 'title' in plot_style:
pass
- elif Budget:
- plot_style['title'] = "GWINC Noise Budget: {}".format(Budget.name)
+ elif budget:
+ plot_style['title'] = "GWINC Noise Budget: {}".format(budget.name)
else:
plot_style['title'] = "GWINC Noise Budget: {}".format(args.IFO)
if args.fom:
logger.info("calculating inspiral {}...".format(range_func))
H = inspiral_range.CBCWaveform(freq, **range_params)
- logger.debug("params: {}".format(H.params))
- fom = eval('inspiral_range.{}'.format(range_func))(freq, traces['Total'][0], H=H)
- logger.info("{}({}) = {:.2f} Mpc".format(range_func, fargs, fom))
- fom_title = 'inspiral {func} {m1}/{m2} $\mathrm{{M}}_\odot$: {fom:.0f} Mpc'.format(
+ logger.debug("waveform params: {}".format(H.params))
+ fom = eval('inspiral_range.{}'.format(range_func))(freq, trace.psd, H=H)
+ logger.info("{}({}) = {:.2f} Mpc".format(range_func, range_func_args, fom))
+ subtitle = 'inspiral {func} {m1}/{m2} $\mathrm{{M}}_\odot$: {fom:.0f} Mpc'.format(
func=range_func,
m1=H.params['m1'],
m2=H.params['m2'],
fom=fom,
)
- plot_style['title'] += '\n{}'.format(fom_title)
+ else:
+ subtitle = None
##########
# interactive
@@ -273,14 +267,14 @@ def main():
if args.interactive:
banner = """GWINC interactive shell
-The 'ifo' Struct and 'traces' data are available for inspection.
-Use the 'whos' command to view the workspace.
+The 'ifo' Struct and 'budget' trace data objects are available for
+inspection. Use the 'whos' command to view the workspace.
"""
if not args.plot:
banner += """
-You may plot the budget using the 'plot_noise()' function:
+You may plot the budget using the 'plot_budget()' function:
-In [.]: plot_noise(freq, traces, **plot_style)
+In [.]: plot_budget(budget, **plot_style)
"""
banner += """
You may interact with the plot using the 'plt' functions, e.g.:
@@ -292,34 +286,31 @@ In [.]: plt.savefig("foo.pdf")
ipshell = InteractiveShellEmbed(
banner1=banner,
user_ns={
- 'freq': freq,
- 'traces': traces,
+ 'budget': trace,
'ifo': ifo,
'plot_style': plot_style,
- 'plot_noise': plot_noise,
+ 'plot_budget': plot_budget,
},
)
ipshell.enable_pylab()
if args.plot:
- ipshell.ex("fig = plot_noise(freq, traces, **plot_style)")
- ipshell.ex("plt.title('{}')".format(plot_style['title']))
+ ipshell.ex("fig = plot_budget(budget, **plot_style)")
+ ipshell.ex("plt.title(plot_style['title'])")
ipshell()
##########
# output
- # save noise traces to HDF5 file
+ # save noise trace to HDF5 file
if out_data_files:
from .io import save_hdf5
- attrs = dict(plot_style)
for path in out_data_files:
- logger.info("saving budget traces: {}".format(path))
+ logger.info("saving budget trace: {}".format(path))
save_hdf5(
path=path,
- freq=freq,
- traces=traces,
+ trace=trace,
ifo=ifo,
- **attrs,
+ plot_style=plot_style,
)
# standard plotting
@@ -327,10 +318,10 @@ In [.]: plt.savefig("foo.pdf")
logger.debug("plotting noises...")
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
- plot_noise(
- freq,
- traces,
+ plot_budget(
+ trace,
ax=ax,
+ subtitle=subtitle,
**plot_style
)
fig.tight_layout()
diff --git a/gwinc/ifo/CE2/__init__.py b/gwinc/ifo/CE2/__init__.py
index 0d0f11eea2e9fc10b5053cdc9d12e6cf53b88ab0..f4e0e69d0008982b6b01d346141f2ac4a10b9fc7 100644
--- a/gwinc/ifo/CE2/__init__.py
+++ b/gwinc/ifo/CE2/__init__.py
@@ -73,7 +73,6 @@ class Substrate(nb.Budget):
noises = [
ITMThermoRefractive,
- ITMCarrierDensity,
SubstrateBrownian,
SubstrateThermoElastic,
]
diff --git a/gwinc/ifo/CE2/ifo.yaml b/gwinc/ifo/CE2/ifo.yaml
index 3c65f15b46e88ac9294a243267eb85dd7e76ec59..97f5f0d70a4c2ab7b8ce3c230acd42fcff23f373 100644
--- a/gwinc/ifo/CE2/ifo.yaml
+++ b/gwinc/ifo/CE2/ifo.yaml
@@ -227,15 +227,6 @@ Materials:
RefractiveIndex: 3.5 # 3.38 * (1 + 4e-5 * T) (ioffe)
dndT: 1e-4 # ~123K & 1900 nm : http://arxiv.org/abs/physics/0606168
Temp: 123 # mirror temperature [K]
- ## parameters for semiconductor optics
- isSemiConductor: True # we are doing semiconductor optics
- CarrierDensity: 1e19 # 1/m^3; carrier density for phosphorous-doped silicon
- ElectronDiffusion: 9.7e-3 # m^2/s; electron diffusion coefficient for silicon at 120 K
- HoleDiffusion: 3.5e-3 # m^2/s; hole diffusion coefficient for silicon at 120 K
- ElectronEffMass: 9.747e-31 # kg; effective mass of each electron 1.07*m_e
- HoleEffMass: 8.016e-31 # kg; effective mass of each hole 0.88*m_e
- ElectronIndexGamma: -8.8e-28 # m**3; dependence of index of refraction on electron carrier density
- HoleIndexGamma: -1.02e-27 # m**3; dependence of index of refraction on hole carrier density
MassRadius: 0.4 # m; 80 cm mCZ silicon
MassThickness: 0.286
@@ -257,7 +248,6 @@ Optics:
# zz = loadmat('CryogenicLIGO/Sensitivity/GWINC/optRuns/' + qopt_mat)
ITM:
SubstrateAbsorption: 1e-3 # 1/m; 10 ppm/cm for MCZ Si
- BeamRadius: 0.141 # m; 1/e^2 power radius w1
CoatingAbsorption: 0.5e-6 # absorption of ITM
Transmittance: 0.014 # zz['x'][0][3]
#CoatingThicknessLown: 0.308
@@ -278,7 +268,6 @@ Optics:
- 0.3867382
- 0.08814237
ETM:
- BeamRadius: 0.141 # m; 1/e^2 power radius w2
Transmittance: 5e-6 # Transmittance of ETM
#CoatingThicknessLown: 0.27
#CoatingThicknessCap: 0.5
diff --git a/gwinc/ifo/Voyager/__init__.py b/gwinc/ifo/Voyager/__init__.py
index 933a7828b72c7dd0cd8de8fb5dcb7414b93ce5b1..eaa433506e907bb55a355c5642bcdfed6e6e60c8 100644
--- a/gwinc/ifo/Voyager/__init__.py
+++ b/gwinc/ifo/Voyager/__init__.py
@@ -14,7 +14,6 @@ class Voyager(nb.Budget):
CoatingBrownian,
CoatingThermoOptic,
ITMThermoRefractive,
- ITMCarrierDensity,
SubstrateBrownian,
SubstrateThermoElastic,
ExcessGas,
diff --git a/gwinc/ifo/Voyager/ifo.yaml b/gwinc/ifo/Voyager/ifo.yaml
index 67f06660e5b11575ac3fd36a6d4a950ebb2a3fdf..81125d1594b1c0169a9fbe77f45ad56a8cc18d9a 100644
--- a/gwinc/ifo/Voyager/ifo.yaml
+++ b/gwinc/ifo/Voyager/ifo.yaml
@@ -252,15 +252,6 @@ Materials:
RefractiveIndex: 3.5 # 3.38 * (1 + 4e-5 * T) (ioffe)
dndT: 1e-4 # ~123K & 1900 nm : http://arxiv.org/abs/physics/0606168
Temp: 123 # mirror temperature [K]
- ## parameters for semiconductor optics
- isSemiConductor: True # we are doing semiconductor optics
- CarrierDensity: 1e19 # 1/m^3; carrier density for phosphorous-doped silicon
- ElectronDiffusion: 9.7e-3 # m^2/s; electron diffusion coefficient for silicon at 120 K
- HoleDiffusion: 3.5e-3 # m^2/s; hole diffusion coefficient for silicon at 120 K
- ElectronEffMass: 9.747e-31 # kg; effective mass of each electron 1.07*m_e
- HoleEffMass: 8.016e-31 # kg; effective mass of each hole 0.88*m_e
- ElectronIndexGamma: -8.8e-28 # m**3; dependence of index of refraction on electron carrier density
- HoleIndexGamma: -1.02e-27 # m**3; dependence of index of refraction on hole carrier density
MassRadius: 0.225 # m; 45 cm mCZ silicon
MassThickness: 0.55
@@ -282,7 +273,6 @@ Optics:
# zz = loadmat('CryogenicLIGO/Sensitivity/GWINC/optRuns/' + qopt_mat)
ITM:
SubstrateAbsorption: 1e-3 # 1/m; 10 ppm/cm for MCZ Si
- BeamRadius: 0.0585 # m; 1/e^2 power radius w1
CoatingAbsorption: 1e-6 # absorption of ITM
Transmittance: 2e-3 # zz['x'][0][3]
#CoatingThicknessLown: 0.308
@@ -303,7 +293,6 @@ Optics:
- 0.386738199718736
- 0.088142365969070
ETM:
- BeamRadius: 0.0835 # m; 1/e^2 power radius w2
Transmittance: 5e-6 # Transmittance of ETM
#CoatingThicknessLown: 0.27
#CoatingThicknessCap: 0.5
diff --git a/gwinc/ifo/__main__.py b/gwinc/ifo/__main__.py
index 1da5c7c0a65ce61ce665c1d3af6b6d40324c01fd..a049331ecd84f152e59a39caeeea618b26a443d0 100644
--- a/gwinc/ifo/__main__.py
+++ b/gwinc/ifo/__main__.py
@@ -29,18 +29,17 @@ def main():
budgets = {}
range_pad = 0
for ifo in IFOS:
- budget = load_budget(ifo)(freq)
+ budget = load_budget(ifo, freq)
name = budget.name
budgets[name] = budget
range_pad = max(len(name), range_pad)
for name, budget in budgets.items():
- budget.update()
- data = budget.calc()
+ trace = budget.run()
try:
import inspiral_range
label_range = ' {:>6.0f} Mpc'.format(
- inspiral_range.range(freq, data),
+ inspiral_range.range(freq, trace.psd),
)
except ModuleNotFoundError:
label_range = ''
@@ -50,7 +49,7 @@ def main():
pad=range_pad,
range=label_range,
)
- ax.loglog(freq, np.sqrt(data), label=label, lw=2)
+ ax.loglog(freq, trace.asd, label=label, lw=2)
ax.grid(
True,
diff --git a/gwinc/ifo/aLIGO/ifo.yaml b/gwinc/ifo/aLIGO/ifo.yaml
index 8005a969d50db68839d9a2a02f8aabdf7cea912b..9ffde7798870a2c68454dae4972bcd34d39deaa0 100644
--- a/gwinc/ifo/aLIGO/ifo.yaml
+++ b/gwinc/ifo/aLIGO/ifo.yaml
@@ -188,7 +188,7 @@ Suspension:
## Optic Material -------------------------------------------------------
Materials:
- MassRadius: 0.17 # m;
+ MassRadius: 0.17 # m;
MassThickness: 0.200 # m; Peter F 8/11/2005
## Dielectric coating material parameters----------------------------------
@@ -251,13 +251,11 @@ Optics:
dc: 1.5707963 # pi/2 # demod/detection/homodyne phase
ITM:
- # BeamRadius: 0.055 # m, 1/e^2 power radius, now in precompIFO
Transmittance: 0.014
CoatingThicknessLown: 0.308
CoatingThicknessCap: 0.5
CoatingAbsorption: 0.5e-6
ETM:
- # BeamRadius: 0.062 # m, 1/e^2 power radius, now in precompIFO
Transmittance: 5e-6
CoatingThicknessLown: 0.27
CoatingThicknessCap: 0.5
diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py
index b760e7ae01ea9535a11d00895567724e6a3d9fb8..f78821ce1bc4cab96ebaf2d629daad0cf49433cd 100644
--- a/gwinc/ifo/noises.py
+++ b/gwinc/ifo/noises.py
@@ -1,3 +1,4 @@
+import copy
import numpy as np
from numpy import pi, sin, exp, sqrt
@@ -20,7 +21,26 @@ def coating_thickness(ifo, optic):
T = optic.Transmittance
dL = optic.CoatingThicknessLown
dCap = optic.CoatingThicknessCap
- return noise.coatingthermal.getCoatDopt(ifo.Materials, T, dL, dCap=dCap)
+ return noise.coatingthermal.getCoatDopt(materials, T, dL, dCap=dCap)
+
+
+def mirror_struct(ifo, tm):
+ """Create a "mirror" Struct for a LIGO core optic
+
+ This is a copy of the ifo.Materials Struct, containing Substrate
+ and Coating sub-Structs, as well as some basic geometrical
+ properties of the optic.
+
+ """
+ # NOTE: we deepcopy this struct since we'll be modifying it (and
+ # 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)
+ mirror.update(optic)
+ mirror.MassVolume = pi * mirror.MassRadius**2 * mirror.MassThickness
+ mirror.MirrorMass = mirror.MassVolume * mirror.Substrate.MassDensity
+ return mirror
def arm_cavity(ifo):
@@ -384,7 +404,7 @@ class Seismic(nb.Noise):
else:
nt, nr = noise.seismic.seismic_BSC_ISI(self.freq)
n, nh, nv = noise.seismic.seismic_suspension_fitered(
- self.ifo.Suspension, nt)
+ self.ifo.Suspension, sustf, nt)
return n * 4
diff --git a/gwinc/io.py b/gwinc/io.py
index fb407486d18a48d133240555c0fa1bcfe99e99ab..98b900e5bd02b4099d982fe9dc8f3b2aa1e2a65d 100644
--- a/gwinc/io.py
+++ b/gwinc/io.py
@@ -4,25 +4,23 @@ import datetime
from . import logger
from . import Struct
+from .trace import BudgetTrace
SCHEMA = 'GWINC noise budget'
-SCHEMA_VERSION = 1
-def _write_trace_recursive(grp, traces):
- for name, trace in traces.items():
- if isinstance(trace, dict):
- tgrp = grp.create_group(name)
- _write_trace_recursive(tgrp, trace)
- else:
- data, style = trace
- dset = grp.create_dataset(name, data=data)
- for key, val in style.items():
- dset.attrs[key] = val
+def _write_trace_recursive(f, trace):
+ f.attrs['style'] = yaml.dump(trace.style)
+ f.create_dataset('PSD', data=trace.psd)
+ if trace.budget:
+ bgrp = f.create_group('budget')
+ for t in trace:
+ tgrp = bgrp.create_group(t.name)
+ _write_trace_recursive(tgrp, t)
-def save_hdf5(path, freq, traces, **kwargs):
+def save_hdf5(path, trace, ifo=None, plot_style=None, **kwargs):
"""Save GWINC budget data to an HDF5 file.
The `freq` argument should be the frequency array, and `traces`
@@ -36,46 +34,107 @@ def save_hdf5(path, freq, traces, **kwargs):
"""
with h5py.File(path, 'w') as f:
f.attrs['SCHEMA'] = SCHEMA
- f.attrs['SCHEMA_VERSION'] = SCHEMA_VERSION
+ f.attrs['SCHEMA_VERSION'] = 2
# FIXME: add budget code hash or something
f.attrs['date'] = datetime.datetime.now().isoformat()
+ if ifo:
+ f.attrs['ifo'] = ifo.to_yaml()
+ if plot_style:
+ f.attrs['plot_style'] = yaml.dump(plot_style)
for key, val in kwargs.items():
- if key == 'ifo':
- f.attrs['ifo'] = val.to_yaml()
+ f.attrs[key] = val
+ f.create_dataset('Freq', data=trace.freq)
+ tgrp = f.create_group(trace.name)
+ _write_trace_recursive(tgrp, trace)
+
+
+def _load_hdf5_v1(f):
+ def read_recursive(element):
+ budget = []
+ for name, item in element.items():
+ if name == 'Total':
+ total = item[:]
+ continue
+ if isinstance(item, h5py.Group):
+ trace = read_recursive(item)
+ trace.name = name
else:
- f.attrs[key] = val
- f.create_dataset('Freq', data=freq)
- tgrp = f.create_group('traces')
- _write_trace_recursive(tgrp, traces)
+ trace = BudgetTrace(
+ name=name,
+ style=dict(item.attrs),
+ psd=item[:],
+ budget=[],
+ )
+ budget.append(trace)
+ return BudgetTrace(
+ psd=total,
+ budget=budget,
+ )
+ trace = read_recursive(f['/traces'])
+ trace.name = 'Total'
+ trace._freq = f['Freq'][:]
+ attrs = dict(f.attrs)
+ if 'ifo' in attrs:
+ try:
+ ifo = Struct.from_yaml(attrs['ifo'])
+ del attrs['ifo']
+ except yaml.constructor.ConstructorError:
+ logger.warning("HDF5 load warning: Could not de-serialize 'ifo' YAML attribute.")
+ ifo = None
+ trace.style = attrs
+ trace.ifo = ifo
+ return trace
-def _read_trace_recursive(element):
- trace = {}
- for name, item in element.items():
- if isinstance(item, h5py.Group):
- trace[name] = _read_trace_recursive(item)
+def _load_hdf5_v2(f):
+ def read_recursive(element):
+ 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.name = name
+ budget.append(trace)
+ return BudgetTrace(
+ psd=psd,
+ budget=budget,
+ style=style,
+ )
+ for name, item in f.items():
+ if name == 'Freq':
+ freq = item[:]
else:
- trace[name] = item[:], dict(item.attrs.items())
+ trace = read_recursive(item)
+ trace.name = name
+ trace._freq = freq
+ trace.ifo = None
+ attrs = dict(f.attrs)
+ ifo = attrs.get('ifo')
+ if ifo:
+ try:
+ trace.ifo = Struct.from_yaml(ifo)
+ 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'])
return trace
def load_hdf5(path):
"""Load GWINC budget data from an HDF5 file.
- Returns (freq, traces, attrs). An 'ifo' attr will be
- de-serialized from YAML into a Struct object.
+ Returns BudgetTrace. An 'ifo' attr will be de-serialized from
+ YAML into a Struct object and assigned as an attribute to the
+ BudgetTrace.
See HDF5_SCHEMA.
"""
+ loaders = {
+ 1: _load_hdf5_v1,
+ 2: _load_hdf5_v2,
+ }
with h5py.File(path, 'r') as f:
- # FIXME: check SCHEMA name/version
- freq = f['Freq'][:]
- traces = _read_trace_recursive(f['/traces'])
- attrs = dict(f.attrs)
- if 'ifo' in attrs:
- try:
- attrs['ifo'] = Struct.from_yaml(attrs['ifo'])
- except yaml.constructor.ConstructorError:
- logger.warning("HDF5 load warning: Could not de-serialize 'ifo' YAML attribute.")
- return freq, traces, attrs
+ version = f.attrs.get('SCHEMA_VERSION', 1)
+ return loaders[version](f)
diff --git a/gwinc/nb.py b/gwinc/nb.py
index 2680862aceeb7c2e751672a0f352674c0d14e749..54391b881c460c34bed8cb4bbbaf7fa6f8dda6ad 100644
--- a/gwinc/nb.py
+++ b/gwinc/nb.py
@@ -1,11 +1,77 @@
-import operator
import itertools
import collections
import numpy as np
import scipy.interpolate
-from functools import reduce
from . import logger
+from .trace import BudgetTrace
+
+
+def precomp(*precomp_funcs, **precomp_fmaps):
+ """BudgetItem.calc decorator to add pre-computed functions
+
+ This is intended to decorate BudgetItem.calc() methods with common
+ functions whose return value are cached for later use. The
+ functions are supplied with the `freq` array and `ifo` Struct
+ attributes as arguments, and their return values are provided as
+ keyword arguments to calc().
+
+ For example, if a calc method is defined as:
+
+ @precomp(foo=precomp_foo)
+ @precomp(bar=precomp_bar)
+ def calc(self, foo, bar):
+ ...
+
+ then when the budget is run the precomp functions will be executed
+ before calc(), roughly the equivalent of:
+
+ foo = precomp_foo(self.freq, self.ifo)
+ bar = precomp_bar(self.freq, self.ifo)
+ psd = calc(foo=foo, bar=bar)
+
+ The execution state of each precomp function is cached so that the
+ same function is not needlessly executed multiple times.
+
+ """
+ def decorator(func):
+ if precomp_funcs:
+ try:
+ func._precomp_list.extend(precomp_funcs)
+ except AttributeError:
+ func._precomp_list = list(precomp_funcs)
+
+ if precomp_fmaps:
+ try:
+ func._precomp_mapped.update(precomp_fmaps)
+ except AttributeError:
+ func._precomp_mapped = dict(precomp_fmaps)
+ return func
+ return decorator
+
+
+def _precomp_recurse_mapping(func, freq, ifo, _precomp):
+ """Recursively execute @precomp decorator functions
+
+ Recurses down functions which may themselves have precomp
+ decorators, building a **kwargs mapping to pass to the wrapped
+ function call.
+
+ """
+ for pc_func in itertools.chain(
+ getattr(func, '_precomp_list', []),
+ getattr(func, '_precomp_mapped', {}).values(),
+ ):
+ if pc_func in _precomp:
+ continue
+ pc_map = _precomp_recurse_mapping(pc_func, freq, ifo, _precomp=_precomp)
+ logger.debug("precomp {}".format(pc_func))
+ _precomp[pc_func] = pc_func(freq, ifo, **pc_map)
+
+ kwargs = {
+ name: _precomp[pc_func] for name, pc_func in getattr(func, '_precomp_mapped', {}).items()
+ }
+ return kwargs
def quadsum(data):
@@ -20,7 +86,6 @@ def quadsum(data):
return np.nansum(data, 0)
-
class BudgetItem:
"""GWINC BudgetItem class
@@ -38,25 +103,10 @@ class BudgetItem:
provided are written directly as attribute variables to the
class (as with __init__).
- Second, after attribute update, all functions listed in the
- `precomp` attribute are executed, supplied with the `freq` and
- `ifo` attributes. So if the `precomp` attribute is defined
- as:
-
- precomp = [precomp_foo, precomp_bar]
-
- then this method will execute the following after attribute
- update:
+ When overloading this method it is recommended to execute the
+ method from the base class with e.g.:
- precomp_foo(self.freq, self.ifo)
- precomp_bar(self.freq, self.ifo)
-
- These functions are intended to update the `ifo` Struct
- attribute with derived values cached for later usage. If the
- `_precomp` argument is provided to this method then it is
- assumed to be a set of previously executed precomp functions,
- and any function already listed there will not be re-executed,
- and it will be updated with any newly executed functions.
+ super().update(**kwargs)
"""
for key, val in kwargs.items():
@@ -79,6 +129,12 @@ class BudgetItem:
"""
return None
+ def _calc(self, _precomp=None):
+ """internal function to call calc with precomp evaluation"""
+ pcmap = _precomp_recurse_mapping(self.calc, self.freq, self.ifo, _precomp)
+ logger.debug("calc {}".format(self.name))
+ return self.calc(**pcmap)
+
##########
def __init__(self, freq=None, **kwargs):
@@ -95,10 +151,10 @@ class BudgetItem:
if freq is not None:
assert isinstance(freq, np.ndarray)
self.freq = freq
- elif not hasattr(self, 'freq'):
- raise AttributeError("Frequency array not provided or defined.")
for key, val in kwargs.items():
setattr(self, key, val)
+ self._loaded = False
+ self._precomp = dict()
@property
def name(self):
@@ -144,7 +200,7 @@ class Calibration(BudgetItem):
e.g. point-by-point product of data and calibration arrays.
"""
- cal = self.calc()
+ cal = self._calc()
assert data.shape == cal.shape, \
"data shape does not match calibration ({} != {})".format(data.shape, cal.shape)
return data * cal
@@ -163,38 +219,75 @@ class Noise(BudgetItem):
style = {}
"""Trace plot style dictionary"""
- def calc_trace(self, calibrations=None, calc=True):
- """Returns noise (PSD, style) tuple.
+ def _make_trace(self, psd=None, budget=None):
+ return BudgetTrace(
+ name=self.name,
+ style=self.style,
+ freq=self.freq,
+ psd=psd,
+ budget=budget,
+ )
+
+ def calc_trace(self, calibration=1, calc=True, _precomp=None):
+ """Calculate noise and return BudgetTrace object
- If `calibrations` is not None it is assumed to be a list of
- len(self.freq) arrays that will be multiplied to the output
- PSD.
+ `calibration` should either be a scalar or a len(self.freq)
+ array that will be multiplied to the output PSD of the budget
+ and all sub noises.
- If calc=False, the noise will not be calculated and the PSD
- will be None. This is useful for just getting the style.
+ If `calc` is False, the noise will not be calculated and the
+ trace PSD will be None. This is useful for just getting the
+ trace style info.
"""
+ if _precomp is None:
+ _precomp = dict()
+
+ total = None
if calc:
- data = self.calc()
- if calibrations:
- data *= reduce(
- operator.mul,
- calibrations,
- )
- else:
- data = None
- return data, self.style
+ total = self._calc(_precomp) * calibration
+
+ return self._make_trace(psd=total)
def run(self, **kwargs):
- """Convenience method to load, update, and return calc traces.
+ """Run full budget and return BudgetTrace object.
+
+ Roughly the equivalent running load(), update(), and
+ calc_trace() in sequence. Keyword arguments are passed to the
+ update() method.
- Equivalent of load(), update(), calc_traces() run in sequence.
- Keyword arguments are passed to update().
+ NOTE: The load status is cached such that subsequent calls to
+ this method will not re-execute the load() method.
+
+ NOTE: The update() method is only run if keyword arguments
+ (`kwargs`) are supplied, or if the `ifo` attribute has
+ changed.
"""
- self.load()
- self.update(**kwargs)
- return self.calc_trace()
+ if not self._loaded:
+ self.load()
+ self._loaded = True
+
+ ifo = kwargs.get('ifo', getattr(self, 'ifo'))
+ if ifo:
+ if not hasattr(ifo, '_orig_keys'):
+ logger.debug("new ifo detected")
+ ifo._orig_keys = tuple(k for k, v in ifo.walk())
+ ifo_hash = ifo.hash()
+ kwargs['ifo'] = ifo
+ else:
+ ifo_hash = ifo.hash(ifo._orig_keys)
+ if ifo_hash != getattr(self, '_ifo_hash', 0):
+ logger.debug("ifo hash change")
+ kwargs['ifo'] = self.ifo
+ self._ifo_hash = ifo_hash
+
+ if kwargs:
+ self.update(**kwargs)
+ # clear precomp cache
+ self._precomp = dict()
+
+ return self.calc_trace(_precomp=self._precomp)
class Budget(Noise):
@@ -391,98 +484,74 @@ class Budget(Noise):
item.load()
def update(self, _precomp=None, **kwargs):
- """Update all noise and cal objects with supplied kwargs."""
- if _precomp is None:
- _precomp = set()
+ """Recursively update all noise and cal objects with supplied kwargs.
+
+ See BudgetItem.update() for more info.
+
+ """
+ for key, val in kwargs.items():
+ setattr(self, key, val)
+
for name, item in itertools.chain(
self._cal_objs.items(),
self._noise_objs.items()):
logger.debug("update {}".format(item))
item.update(_precomp=_precomp, **kwargs)
- def calc_noise(self, name):
- """Return calibrated individual noise term.
+ def calc_noise(self, name, calibration=1, calc=True, _cals=None, _precomp=None):
+ """Return calibrated individual noise BudgetTrace.
- The noise PSD and calibration transfer functions are
- calculated, and the calibrated noise array is returned.
+ The noise and calibration transfer functions are calculated,
+ and the calibrated noise BudgetTrace is returned.
+ `calibration` is an overall calculated calibration to apply to
+ the noise.
"""
+ if _precomp is None:
+ _precomp = dict()
+ for cal in self._noise_cals[name]:
+ if _cals:
+ calibration *= _cals[cal]
+ else:
+ obj = self._cal_objs[cal]
+ calibration *= obj._calc(_precomp)
noise = self._noise_objs[name]
- cals = [self._cal_objs[cal].calc() for cal in self._noise_cals[name]]
- data = noise.calc_trace(cals)
- if isinstance(data, dict):
- return data['Total'][0]
- else:
- return data[0]
-
- def calc(self):
- """Calculate sum of all noises.
-
- """
- data = [self.calc_noise(name) for name in self._noise_objs.keys() if name in self._budget_noises]
- return quadsum(data)
-
- def calc_trace(self, calibrations=None, calc=True):
- """Returns a dictionary of noises traces, keyed by noise names.
+ return noise.calc_trace(
+ calibration=calibration,
+ calc=calc,
+ _precomp=_precomp,
+ )
- Values are (data, style) trace tuples (see Noise.calc_trace).
- The key of the budget sum total is 'Total'. The values of sub
- budgets are themselves dictionaries returned from
- calc_trace() of the sub budget.
+ def calc_trace(self, calibration=1, calc=True, _precomp=None):
+ """Calculate all budget noises and return BudgetTrace object
- If `calibrations` is not None it is assumed to be a list of
- len(self.freq) array that will be multiplied to the output PSD
- of the budget and all sub noises.
+ `calibration` should either be a scalar or a len(self.freq)
+ array that will be multiplied to the output PSD of the budget
+ and all sub noises.
- If calc=False, the noise will not be calculated and the PSD
- will be None. This is useful for just getting style the
- style.
+ If `calc` is False, the noise will not be calculated and the
+ trace PSD will be None. This is useful for just getting the
+ trace style info.
"""
- # start by creating an empty OrderedDict used for outputing trace data
- # or style info with the following order:
- # references
- # total
- # constituents
- d = collections.OrderedDict()
- # allocate references
- for name, noise in self._noise_objs.items():
- if name in self._budget_noises:
- continue
- d[name] = noise.calc_trace(calc=False)
- # allocate total
- if self._budget_noises:
- d['Total'] = None, self.style
- # allocate constituent
- for name, noise in self._noise_objs.items():
- if name not in self._budget_noises:
- continue
- d[name] = noise.calc_trace(calc=False)
- # if we're not calc'ing, just return the dict now
- if not calc:
- return d
-
- # calc all calibrations
- c = {}
- for name, cal in self._cal_objs.items():
- logger.debug("calc {}".format(name))
- c[name] = cal.calc()
- # calc all noises
- for name, noise in self._noise_objs.items():
- cals = [c[cal] for cal in self._noise_cals[name]]
- if calibrations:
- cals += calibrations
- logger.debug("calc {}".format(name))
- d[name] = noise.calc_trace(
- calibrations=cals,
+ if _precomp is None:
+ _precomp = dict()
+
+ _cals = {}
+ if calc:
+ for name, cal in self._cal_objs.items():
+ _cals[name] = cal._calc(_precomp)
+ budget = []
+ for name in self._noise_objs:
+ trace = self.calc_noise(
+ name,
+ calibration=calibration,
+ calc=calc,
+ _cals=_cals,
+ _precomp=_precomp,
)
- # calc budget total
- constituent_data = []
- for name in self._budget_noises:
- if isinstance(d[name], dict):
- data = d[name]['Total'][0]
- else:
- data = d[name][0]
- constituent_data.append(data)
- d['Total'] = quadsum(constituent_data), self.style
- return d
+ budget.append(trace)
+ total = quadsum([trace.psd for trace in budget])
+ return self._make_trace(
+ psd=total, budget=budget
+ )
diff --git a/gwinc/noise/coatingthermal.py b/gwinc/noise/coatingthermal.py
index 4de77536050d5b7eaae55b1f21282a91767d7de6..a8c191f5cd38876cfdd3470032f11b5fa46ecc21 100644
--- a/gwinc/noise/coatingthermal.py
+++ b/gwinc/noise/coatingthermal.py
@@ -10,57 +10,104 @@ from ..const import BESSEL_ZEROS as zeta
from ..const import J0M as j0m
-def coating_brownian(f, materials, wavelength, wBeam, dOpt):
- """Optical coating Brownian thermal displacement noise spectrum
+def coating_brownian(f, mirror, wavelength, wBeam, power):
+ """Coating brownian noise for a given collection of coating layers
- :f: frequency array in Hz
- :materials: gwinc optic materials structure
- :wavelength: laser wavelength
- :wBeam: beam radius (at 1 / e^2 power)
- :dOpt: coating layer thickness array (Nlayer x 1)
+ This function calculates Coating Brownian noise using
+ Hong et al . PRD 87, 082001 (2013).
+ All references to 'the paper', 'Eq' and 'Sec' are to this paper.
+
+ ***Important Note***
+ Inside this function phi is used for denoting the phase shift suffered
+ by light in one way propagation through a layer. This is in conflict
+ with present nomenclature everywhere else where it is used as loss angle.
The layers are assumed to be alernating low-n high-n layers, with
low-n first.
- :returns: displacement noise power spectrum at :f:
-
- Added by G Harry 8/3/02 from work by Nakagawa, Gretarsson, et al.
- Expanded to reduce approximation, GMH 8/03
- Modified to return strain noise, PF 4/07
- Modified to accept coating with non-quater-wave layers, mevans 25 Apr 2008
+ Inputs:
+ f = frequency vector in Hz
+ mirror = mirror properties Struct
+ wavelength = laser wavelength
+ wBeam = beam radius (at 1 / e**2 power)
+ power = laser power falling on the mirror (W)
+
+ If the mirror.Coating.IncCoatBrAmpNoise parameter is present and
+ evaluates to True the amplitude noise due to coating brownian
+ noise will be calculated and its effect on the phase noise will be
+ added. In that case the following new optional arguments should
+ be made available in the Materials object to provide separate Bulk
+ and Shear loss angles and to observe effect of photoelasticity:*
+ lossBlown = Coating Bulk Loss Angle of Low Refrac.Index layer @ 100Hz
+ lossSlown = Coating Shear Loss Angle of Low Refrac. Index layer @ 100Hz
+ lossBhighn = Coating Bulk Loss Angle of High Refrac. Index layer @ 100Hz
+ lossShighn = Coating Shear Loss Angle of High Refrac. Index layer @ 100Hz
+ lossBlown_slope = Coating Bulk Loss Angle Slope of Low Refrac. Index layer
+ lossSlown_slope = Coating Shear Loss Angle Slope of Low Refrac. Index layer
+ lossBhighn_slope = Coating Bulk Loss Angle Slope of High Refrac. Index layer
+ lossShighn_slope = Coating Shear Loss Angle Slope of High Refrac. Index layer
+ PETlown = Relevant component of Photoelastic Tensor of High n layer*
+ PEThighn = Relevant component of Photoelastic Tensor of Low n layer*
+
+ Returns:
+ SbrZ = Brownian noise spectra for one mirror in m**2 / Hz
+
+ *
+ Choice of PETlown and PEThighn can be inspired from sec. A.1. of the paper.
+ There, values are chosen to get the longitudnal coefficent of
+ photoelasticity as -0.5 for Tantala and -0.27 for Silica.
+ These values also need to be added in Materials object.
+ *
+ If the optional arguments are not present, Phihighn and Philown will be
+ used as both Bulk and Shear loss angles and PET coefficients will be set
+ to 0.
"""
# extract substructures
- sub = materials.Substrate
- coat = materials.Coating
+ sub = mirror.Substrate
+ coat = mirror.Coating
# Constants
kBT = const.kB * sub.Temp
+ c = const.c
# substrate properties
- Ysub = sub.MirrorY # Young's Modulous
- pratsub = sub.MirrorSigma # Poisson Ratio
+ Ysub = sub.MirrorY # Young's Modulous
+ pratsub = sub.MirrorSigma # Poisson Ratio
+ nsub = sub.RefractiveIndex # Refractive Index
# coating properties
+ dOpt = mirror.Coating.dOpt
Yhighn = coat.Yhighn
sigmahighn = coat.Sigmahighn
- phihighn = coat.Phihighn
nH = coat.Indexhighn
+ if 'PEThighn' in coat:
+ PEThighn = coat.PEThighn
+ else:
+ PEThighn = 0.0
Ylown = coat.Ylown
sigmalown = coat.Sigmalown
- philown = coat.Philown
nL = coat.Indexlown
+ if 'PETlown' in coat:
+ PETlown = coat.PETlown
+ else:
+ PETlown = 0.0
+
+ # Loss Angles
+ lossBhighn, lossShighn, lossBlown, lossSlown = interpretLossAngles(coat)
# make vectors of material properties
- nN = np.zeros(len(dOpt))
- yN = np.zeros(len(dOpt))
- pratN = np.zeros(len(dOpt))
- phiN = np.zeros(len(dOpt))
+ nol = len(dOpt) # Number of layers
+ nN = np.zeros(nol)
+ yN = np.zeros(nol)
+ pratN = np.zeros(nol)
+ CPE = np.zeros(nol)
# make simple alternating structure (low-index, high-index doublets)
# (adapted from the more general calculation in
- # Yam, W., Gras, S., & Evans, M. Multimaterial coatings with reduced thermal noise.
+ # Yam, W., Gras, S., & Evans, M. Multimaterial coatings with reduced
+ # thermal noise.
# Physical Review D, 91(4), 042002. (2015).
# http://doi.org/10.1103/PhysRevD.91.042002 )
nN[::2] = nL
@@ -72,59 +119,147 @@ def coating_brownian(f, materials, wavelength, wBeam, dOpt):
pratN[::2] = sigmalown
pratN[1::2] = sigmahighn
- phiN[::2] = philown
- phiN[1::2] = phihighn
+ # Coefficient of photoelastic effect
+ # PRD 87, 082001 Eq (A6)
+ CPE[::2] = -0.5*PETlown*nL**3
+ CPE[1::2] = -0.5*PEThighn*nH**3
# geometrical thickness of each layer and total
dGeo = wavelength * np.asarray(dOpt) / nN
- #dCoat = np.sum(dGeo)
-
- ###################################
- # Brownian
- ###################################
- # coating reflectivity
- dcdp = getCoatRefl(materials, dOpt)[1]
-
- # Z-dir (1 => away from the substrate, -1 => into the substrate)
- zdir = -1
- dcdp_z = zdir * dcdp # z-dir only matters here
-
- # layer contributions, b_j (eq 1) from doi:10.1103/PhysRevD.91.042002, errors corrected
- brLayer = ( 1/(1-pratN) *
- ( (1-nN*dcdp_z)**2 * (1-2*pratN)*(1+pratN)*Ysub / ((1-2*pratsub)*(1+pratsub)*yN) +
- (1-2*pratsub)*(1+pratsub)*yN / ((1+pratN)*Ysub) ) )
- # sum them up for total
- w = 2 * pi * f
-
- is_low_slope = 'Philown_slope' in coat
- is_high_slope = 'Phihighn_slope' in coat
-
- if (not is_low_slope) and (not is_high_slope):
- # this is the old code for frequency independent loss
- SbrZ = (4 * kBT / (pi * wBeam**2 * w)) * \
- sum(dGeo * brLayer * phiN) * (1 - pratsub - 2 * pratsub**2) / Ysub
+ # WaveLength of light in each layer
+ lambdaN = wavelength / nN
+
+ # Calculate rho and derivatives of rho
+ # with respect to both phi_k and r_j
+ rho, dLogRho_dPhik, dLogRho_dRk, r = getCoatReflAndDer(nN, nsub, dOpt)
+
+ # Define the function epsilon as per Eq (25)
+ # Split epsilon function as:
+ # Epsilon = Ep1 - Ep2 * cos(2k0n(z-zjp1)) - Ep3 * sin(2k0n(z-zjp1))
+
+ # Part 1 of epsilon function
+ Ep1 = (nN + CPE) * dLogRho_dPhik[:-1]
+ # Part 2 of epsilon function (Prefactor of cosine term)
+ Ep2 = CPE * (dLogRho_dPhik[:-1] * (1 - r[:-1]**2) / (2*r[:-1])
+ - (dLogRho_dPhik[1:] * (1 + r[:-1]**2) / (2 * r[:-1])))
+ # Part 3 of epsilon function (Prefactor of sine term)
+ Ep3 = (1 - r[:-1]**2) * CPE * dLogRho_dRk[:-1]
+
+ # Define (1 - Im(epsilon)/2)
+ Ip1 = 1 - imag(Ep1) / 2 # First part of (1 - Im(epsilon)/2)
+ Ip2 = imag(Ep2) / 2 # Prefactor to cosine in (1 - Im(epsilon)/2)
+ Ip3 = imag(Ep3) / 2 # Prefactor to sine in (1 - Im(epsilon)/2)
+
+ # Define transfer functions from bulk and shear noise fields to layer
+ # thickness and surface height as per Table I in paper
+ C_B = np.sqrt(0.5*(1+pratN))
+ C_SA = np.sqrt(1 - 2*pratN)
+ D_B = ((1 - pratsub - 2*pratsub**2)*yN/(np.sqrt(2*(1+pratN))*Ysub))
+ D_SA = -((1 - pratsub - 2*pratsub**2)*yN/(2*np.sqrt(1-2*pratN)*Ysub))
+ D_SB = (np.sqrt(3)*(1-pratN)*(1 - pratsub - 2*pratsub**2)*yN
+ / (2*np.sqrt(1-2*pratN)*(1+pratN)*Ysub))
+
+ # Calculating effective beam area on each layer
+ # Assuming the beam radius does not change significantly through the
+ # depth of the mirror.
+ Aeff = pi*(wBeam**2)
+
+ # PSD at single layer with thickness equal to WaveLength
+ # in the medium Eq (96)
+ S_Bk = np.zeros((nol, len(f)))
+ S_Bk[::2, :] = (4 * kBT * lambdaN[0] * lossBlown(f)
+ * (1 - pratN[0] - 2 * pratN[0]**2)
+ / (3 * pi * f * yN[0] * ((1 - pratN[0])**2) * Aeff))
+ S_Bk[1::2, :] = (4 * kBT * lambdaN[1] * lossBhighn(f)
+ * (1 - pratN[1] - 2 * pratN[1]**2)
+ / (3 * pi * f * yN[1] * ((1 - pratN[1])**2) * Aeff))
+
+ # Shear
+ S_Sk = np.zeros((nol, len(f)))
+ S_Sk[::2, :] = (4 * kBT * lambdaN[0] * lossSlown(f)
+ * (1 - pratN[0] - 2 * pratN[0]**2)
+ / (3 * pi * f * yN[0] * ((1 - pratN[0])**2) * Aeff))
+ S_Sk[1::2, :] = (4 * kBT * lambdaN[1] * lossShighn(f)
+ * (1 - pratN[1] - 2 * pratN[1]**2)
+ / (3 * pi * f * yN[1] * ((1 - pratN[1])**2) * Aeff))
+
+ # Coefficients q_j from Eq (94
+ # See https://dcc.ligo.org/T2000552 for derivation
+ k0 = 2 * pi / wavelength
+ q_Bk = (+ 8 * C_B * (D_B + C_B * Ip1) * Ip3
+ + 2 * C_B**2 * Ip2 * Ip3
+ + 4 * (2 * D_B**2 + 4 * C_B * D_B * Ip1
+ + C_B**2 * (2 * Ip1**2 + Ip2**2 + Ip3**2)
+ ) * k0 * nN * dGeo
+ - 8 * C_B * (D_B + C_B * Ip1) * Ip3 * np.cos(2 * k0 * nN * dGeo)
+ - 2 * C_B**2 * Ip2 * Ip3 * np.cos(4 * k0 * nN * dGeo)
+ + 8 * C_B * (D_B + C_B * Ip1) * Ip2 * np.sin(2 * k0 * nN * dGeo)
+ + C_B**2 * (Ip2 - Ip3) * (Ip2 + Ip3) * np.sin(4 * k0 * nN * dGeo)
+ ) / (8 * k0 * lambdaN * nN)
+
+ q_Sk = (D_SB**2 * 8 * k0 * nN * dGeo
+ + 8 * C_SA * (D_SA + C_SA * Ip1) * Ip3
+ + 2 * C_SA**2 * Ip2 * Ip3
+ + 4 * (2 * D_SA**2 + 4 * C_SA * D_SA * Ip1
+ + C_SA**2 * (2 * Ip1**2 + Ip2**2 + Ip3**2)
+ ) * k0 * nN * dGeo
+ - 8 * C_SA * (D_SA + C_SA * Ip1) * Ip3 * np.cos(2 * k0 * nN * dGeo)
+ - 2 * C_SA**2 * Ip2 * Ip3 * np.cos(4 * k0 * nN * dGeo)
+ + 8 * C_SA * (D_SA + C_SA * Ip1) * Ip2 * np.sin(2 * k0 * nN * dGeo)
+ + C_SA**2 * (Ip2 - Ip3) * (Ip2 + Ip3) * np.sin(4 * k0 * nN * dGeo)
+ ) / (8 * k0 * lambdaN * nN)
+
+ # S_Xi as per Eq(94)
+ S_Xi = (np.tensordot(q_Bk, S_Bk, axes=1)
+ + np.tensordot(q_Sk, S_Sk, axes=1))
+
+ # From Sec II.E. Eq.(41)
+ # Conversion of brownian amplitude noise to displacement noise
+ if coat.get('IncCoatBrAmpNoise'):
+
+ # get/calculate optic transmittance
+ mTi = mirror.get('Transmittance', 1-np.abs(rho)**2)
+
+ # Define Re(epsilon)/2
+ Rp1 = np.real(Ep1) / 2 # First part of Re(epsilon)/2
+ Rp2 = -np.real(Ep2) / 2 # Prefactor to cosine in Re(epsilon)/2
+ Rp3 = -np.real(Ep3) / 2 # Prefactor to sine in Re(epsilon)/2
+ # Coefficients p_j from Eq (95)
+ # See https://dcc.ligo.org/T2000552 for derivation
+ p_BkbyC = (+ 8 * Rp1 * Rp3
+ + 2 * Rp2 * Rp3
+ + 4 * (2 * Rp1**2 + Rp2**2 + Rp3**2) * k0 * nN * dGeo
+ - 8 * Rp1 * Rp3 * np.cos(2 * k0 * nN * dGeo)
+ - 2 * Rp2 * Rp3 * np.cos(4 * k0 * nN * dGeo)
+ + 8 * Rp1 * Rp2 * np.sin(2 * k0 * nN * dGeo)
+ + (Rp2 - Rp3) * (Rp2 + Rp3) * np.sin(4 * k0 * nN * dGeo)
+ ) / (8 * k0 * lambdaN * nN)
+ p_Bk = p_BkbyC * C_B**2
+ p_Sk = p_BkbyC * C_SA**2
+
+ # S_Zeta as per Eq(95)
+ S_Zeta = (np.tensordot(p_Bk, S_Bk, axes=1)
+ + np.tensordot(p_Sk, S_Sk, axes=1))
+
+ AmpToDispConvFac = ((32 * power)
+ / (mirror.MirrorMass * wavelength
+ * f**2 * c * 2 * pi * sqrt(mTi)))
+ # Adding the two pathways of noise contribution as correlated noise
+ SbrZ = (sqrt(S_Xi) + AmpToDispConvFac * sqrt(S_Zeta))**2
else:
- SbrZ = (4 * kBT / (pi * wBeam**2 * w)) * \
- (1 - pratsub - 2 * pratsub**2) / Ysub
- SbrZ *= (np.sum(dGeo[::2] * brLayer[::2] * phiN[::2]) * (f / 100)**(coat.Philown_slope) +
- np.sum(dGeo[1::2] * brLayer[1::2] * phiN[1::2]) * (f / 100)**(coat.Phihighn_slope))
-
- # for the record: evaluate summation in eq 1 of PhysRevD.91.042002
- # normalized by total coating thickness to make it unitless
- #cCoat = np.sum(dGeo * brLayer * phiN) / dCoat
+ SbrZ = S_Xi
return SbrZ
-def coating_thermooptic(f, materials, wavelength, wBeam, dOpt):
+def coating_thermooptic(f, mirror, wavelength, wBeam):
"""Optical coating thermo-optic displacement noise spectrum
:f: frequency array in Hz
- :materials: gwinc optic materials structure
+ :mirror: mirror parameter Struct
:wavelength: laser wavelength
- :wBeam: beam radius (at 1 / e^2 power)
- :dOpt: coating layer thickness array (Nlayer x 1)
+ :wBeam: beam radius (at 1 / e**2 power)
:returns: tuple of:
StoZ = displacement noise power spectrum at :f:
@@ -134,15 +269,15 @@ def coating_thermooptic(f, materials, wavelength, wBeam, dOpt):
"""
# compute coefficients
- dTO, dTR, dTE, T, junk = getCoatTOPos(materials, wavelength, wBeam, dOpt)
+ dTO, dTR, dTE, T, junk = getCoatTOPos(mirror, wavelength, wBeam)
# compute correction factors
- gTO = getCoatThickCorr(f, materials, wavelength, dOpt, dTE, dTR)
- gTE = getCoatThickCorr(f, materials, wavelength, dOpt, dTE, 0)
- gTR = getCoatThickCorr(f, materials, wavelength, dOpt, 0, dTR)
+ gTO = getCoatThickCorr(f, mirror, wavelength, dTE, dTR)
+ gTE = getCoatThickCorr(f, mirror, wavelength, dTE, 0)
+ gTR = getCoatThickCorr(f, mirror, wavelength, 0, dTR)
# compute thermal source spectrum
- SsurfT, junk = getCoatThermal(f, materials, wBeam)
+ SsurfT, junk = getCoatThermal(f, mirror, wBeam)
StoZ = SsurfT * gTO * dTO**2
SteZ = SsurfT * gTE * dTE**2
@@ -151,13 +286,12 @@ def coating_thermooptic(f, materials, wavelength, wBeam, dOpt):
return (StoZ, SteZ, StrZ, T)
-def getCoatTOPos(materials, wavelength, wBeam, dOpt):
+def getCoatTOPos(mirror, wavelength, wBeam):
"""Mirror position derivative wrt thermal fluctuations
- :materials: gwinc optic materials structure
+ :mirror: mirror parameter Struct
:wavelength: laser wavelength
- :wBeam: beam radius (at 1 / e^2 power)
- :dOpt: coating layer thickness array (Nlayer x 1)
+ :wBeam: beam radius (at 1 / e**2 power)
:returns: tuple of:
dTO = total thermo-optic dz/dT
@@ -172,13 +306,14 @@ def getCoatTOPos(materials, wavelength, wBeam, dOpt):
"""
# parameters
- nS = materials.Substrate.RefractiveIndex
+ nS = mirror.Substrate.RefractiveIndex
+ dOpt = mirror.Coating.dOpt
# compute refractive index, effective alpha and beta
- nLayer, aLayer, bLayer, dLayer, sLayer = getCoatLayers(materials, wavelength, dOpt)
+ nLayer, aLayer, bLayer, dLayer, sLayer = getCoatLayers(mirror, wavelength)
# compute coating average parameters
- dc, Cc, Kc, aSub = getCoatAvg(materials, wavelength, dOpt)
+ dc, Cc, Kc, aSub = getCoatAvg(mirror, wavelength)
# compute reflectivity and parameters
dphi_dT, dphi_TE, dphi_TR, rCoat = getCoatTOPhase(1, nS, nLayer, dOpt, aLayer, bLayer, sLayer)
@@ -193,7 +328,7 @@ def getCoatTOPos(materials, wavelength, wBeam, dOpt):
dTE = dphi_TE * wavelength / (4 * pi) - aSub * dc
# mirror finite size correction
- Cfsm = getCoatFiniteCorr(materials, wavelength, wBeam, dOpt)
+ Cfsm = getCoatFiniteCorr(mirror, wavelength, wBeam)
dTE = dTE * Cfsm
# add TE and TR effects (sign is already included)
@@ -202,14 +337,13 @@ def getCoatTOPos(materials, wavelength, wBeam, dOpt):
return dTO, dTR, dTE, T, R
-def getCoatThickCorr(f, materials, wavelength, dOpt, dTE, dTR):
+def getCoatThickCorr(f, mirror, wavelength, dTE, dTR):
"""Finite coating thickness correction
:f: frequency array in Hz
- :materials: gwinc optic materials structure
+ :mirror: gwinc optic mirror structure
:wavelength: laser wavelength
- :wBeam: beam radius (at 1 / e^2 power)
- :dOpt: coating layer thickness array (Nlayer x 1)
+ :wBeam: beam radius (at 1 / e**2 power)
Uses correction factor from LIGO-T080101, "Thick Coating
Correction" (Evans).
@@ -221,18 +355,18 @@ def getCoatThickCorr(f, materials, wavelength, dOpt, dTE, dTR):
# For comparison in the bTR = 0 limit, the
# equation from Fejer (PRD70, 2004)
# modified so that gFC -> 1 as xi -> 0
- # gTC = (2 ./ (R * xi.^2)) .* (sh - s + R .* (ch - c)) ./ ...
- # (ch + c + 2 * R * sh + R^2 * (ch - c));
+ # gTC = (2 ./ (R * xi.**2)) .* (sh - s + R .* (ch - c)) ./ ...
+ # (ch + c + 2 * R * sh + R**2 * (ch - c));
# which in the limit of xi << 1 becomes
# gTC = 1 - xi * (R - 1 / (3 * R));
# parameter extraction
- pS = materials.Substrate
+ pS = mirror.Substrate
Cs = pS.MassCM * pS.MassDensity
Ks = pS.MassKappa
# compute coating average parameters
- dc, Cc, Kc, junk = getCoatAvg(materials, wavelength, dOpt)
+ dc, Cc, Kc, junk = getCoatAvg(mirror, wavelength)
# R and xi (from T080101, Thick Coating Correction)
w = 2 * pi * f
@@ -261,19 +395,19 @@ def getCoatThickCorr(f, materials, wavelength, dOpt, dTE, dTR):
return gTC
-def getCoatThermal(f, materials, wBeam):
+def getCoatThermal(f, mirror, wBeam):
"""Thermal noise spectra for a surface layer
:f: frequency array in Hz
- :materials: gwinc optic material structure
- :wBeam: beam radius (at 1 / e^2 power)
+ :mirror: mirror parameter Struct
+ :wBeam: beam radius (at 1 / e**2 power)
:returns: tuple of:
- SsurfT = power spectra of thermal fluctuations in K^2 / Hz
+ SsurfT = power spectra of thermal fluctuations in K**2 / Hz
rdel = thermal diffusion length at each frequency in m
"""
- pS = materials.Substrate
+ pS = mirror.Substrate
C_S = pS.MassCM * pS.MassDensity
K_S = pS.MassKappa
kBT2 = const.kB * pS.Temp**2
@@ -290,12 +424,11 @@ def getCoatThermal(f, materials, wBeam):
return SsurfT, rdel
-def getCoatLayers(materials, wavelength, dOpt):
+def getCoatLayers(mirror, wavelength):
"""Layer vectors for refractive index, effective alpha and beta and geometrical thickness
- :materials: gwinc optic materials structure
+ :mirror: mirror parameter Struct
:wavelength: laser wavelength
- :dOpt: coating layer thickness array (Nlayer x 1)
:returns: tuple of:
nLayer = refractive index of each layer, ordered input to output (N x 1)
@@ -309,8 +442,9 @@ def getCoatLayers(materials, wavelength, dOpt):
"""
# coating parameters
- pS = materials.Substrate
- pC = materials.Coating
+ pS = mirror.Substrate
+ pC = mirror.Coating
+ dOpt = mirror.Coating.dOpt
Y_S = pS.MirrorY
sigS = pS.MirrorSigma
@@ -364,12 +498,11 @@ def getCoatLayers(materials, wavelength, dOpt):
return nLayer, aLayer, bLayer, dLayer, sLayer
-def getCoatAvg(materials, wavelength, dOpt):
+def getCoatAvg(mirror, wavelength):
"""Coating average properties
- :materials: gwinc optic materials structure
+ :mirror: gwinc optic mirror structure
:wavelength: laser wavelength
- :dOpt: coating layer thickness array (Nlayer x 1)
:returns: tuple of:
dc = total thickness (meters)
@@ -379,8 +512,9 @@ def getCoatAvg(materials, wavelength, dOpt):
"""
# coating parameters
- pS = materials.Substrate
- pC = materials.Coating
+ pS = mirror.Substrate
+ pC = mirror.Coating
+ dOpt = mirror.Coating.dOpt
alphaS = pS.MassAlpha
C_S = pS.MassCM * pS.MassDensity
@@ -393,7 +527,7 @@ def getCoatAvg(materials, wavelength, dOpt):
K_H = pC.ThermalDiffusivityhighn
# compute refractive index, effective alpha and beta
- junk1, junk2, junk3, dLayer, junk4 = getCoatLayers(materials, wavelength, dOpt)
+ junk1, junk2, junk3, dLayer, junk4 = getCoatLayers(mirror, wavelength)
# heat capacity
dc = np.sum(dLayer)
@@ -423,7 +557,7 @@ def getCoatTOPhase(nIn, nOut, nLayer, dOpt, aLayer, bLayer, sLayer):
:aLayer: change in geometrical thickness with temperature
= the effective thermal expansion coeffient of the coating layer
:bLayer: change in refractive index with temperature
- = dn/dT
+ = dn/dT
= dd/dT - n * a
:returns: tuple of:
@@ -463,13 +597,12 @@ def getCoatTOPhase(nIn, nOut, nLayer, dOpt, aLayer, bLayer, sLayer):
return dphi_dT, dphi_TE, dphi_TR, rCoat
-def getCoatFiniteCorr(materials, wavelength, wBeam, dOpt):
+def getCoatFiniteCorr(mirror, wavelength, wBeam):
"""Finite mirror size correction
- :materials: gwinc optic materials structure
+ :mirror: mirror parameter Struct
:wavelength: laser wavelength
- :wBeam: beam radius (at 1 / e^2 power)
- :dOpt: coating layer thickness array (Nlayer x 1)
+ :wBeam: beam radius (at 1 / e**2 power)
Uses correction factor from PLA 2003 vol 312 pg 244-255
"Thermodynamical fluctuations in optical mirror coatings"
@@ -482,25 +615,26 @@ def getCoatFiniteCorr(materials, wavelength, wBeam, dOpt):
"""
# parameter extraction
- R = materials.MassRadius
- H = materials.MassThickness
-
- alphaS = materials.Substrate.MassAlpha
- C_S = materials.Substrate.MassCM * materials.Substrate.MassDensity
- Y_S = materials.Substrate.MirrorY
- sigS = materials.Substrate.MirrorSigma
-
- alphaL = materials.Coating.Alphalown
- C_L = materials.Coating.CVlown
- Y_L = materials.Coating.Ylown
- sigL = materials.Coating.Sigmalown
- nL = materials.Coating.Indexlown
-
- alphaH = materials.Coating.Alphahighn
- C_H = materials.Coating.CVhighn
- Y_H = materials.Coating.Yhighn
- sigH = materials.Coating.Sigmahighn
- nH = materials.Coating.Indexhighn
+ R = mirror.MassRadius
+ H = mirror.MassThickness
+ dOpt = mirror.Coating.dOpt
+
+ alphaS = mirror.Substrate.MassAlpha
+ C_S = mirror.Substrate.MassCM * mirror.Substrate.MassDensity
+ Y_S = mirror.Substrate.MirrorY
+ sigS = mirror.Substrate.MirrorSigma
+
+ alphaL = mirror.Coating.Alphalown
+ C_L = mirror.Coating.CVlown
+ Y_L = mirror.Coating.Ylown
+ sigL = mirror.Coating.Sigmalown
+ nL = mirror.Coating.Indexlown
+
+ alphaH = mirror.Coating.Alphahighn
+ C_H = mirror.Coating.CVhighn
+ Y_H = mirror.Coating.Yhighn
+ sigH = mirror.Coating.Sigmahighn
+ nH = mirror.Coating.Indexhighn
# coating sums
dL = wavelength * np.sum(dOpt[::2]) / nL
@@ -536,7 +670,7 @@ def getCoatFiniteCorr(materials, wavelength, wBeam, dOpt):
# between eq 77 and 78
km = zeta / R
Qm = exp(-2 * km * H)
- pm = exp(-km**2 * r0**2 / 4) / j0m # left out factor of pi * R^2 in denominator
+ pm = exp(-km**2 * r0**2 / 4) / j0m # left out factor of pi * R**2 in denominator
# eq 88
Lm = Xr - Zf * (1 + sigS) + (Yr * (1 - 2 * sigS) + Zf - 2 * Cr) * \
@@ -754,3 +888,141 @@ def getCoatRefl2(nIn, nOut, nLayer, dOpt):
dcdp = -imag(dr_dphi / rCoat) ### Where did this minus come from???
return rCoat, dcdp, rbar, r
+
+
+def getCoatReflAndDer(nN, nsub, dOpt):
+ '''
+ Helper function for coating_brownian_hong().
+ Follows Hong et al . PRD 87, 082001 (2013) Sec V.A.
+ This function calculates derivatives of complex reflectivity of Coating
+ with respect to phase shifts through each layer and reflectivities of
+ each interface
+ Input:
+
+ nN = Refractive indices of coatings layers
+ nsub = Refractive Index of Substrate
+ dOpt = optical thickness / lambda of each layer,
+ geometrical thickness * refractive index / lambda
+
+ Returns:
+ delLogRho_delPhik = Partial derivative of log of total effective
+ reflectivity of coating with respect to phase shifts
+ in each layer.
+ delLogRho_delReflk = Partial derivative of log of total effective
+ reflectivity of coating with respect to reflectivity
+ of each interface.
+ '''
+ nol = len(dOpt) # Number of layers in coating
+ # Reflectivities and transmitivities
+ # r[j] is reflectivity from (j-1)th and (j)th layer interface
+ # Here r[0] is reflectivity from air and 0th layer
+ # and r[-1] is reflectivity between last layer and substrate
+ Refl = np.zeros(nol+1)
+ Refl[0] = (1 - nN[0]) / (1 + nN[0])
+ Refl[1:-1] = (nN[:-1] - nN[1:]) / (nN[:-1] + nN[1:])
+ Refl[-1] = (nN[-1] - nsub) / (nN[-1] + nsub)
+ # Note the shift from nomenclature
+ # Phi is reserved for denoting one-way phase shift suffered by light
+ # during propagation through a layer
+ Phi = np.asarray(dOpt) * 2 * pi
+
+ # Define rho_k as reflectivity of
+ # k layers starting from (N-k-1)th lyer to (N-1)th layer
+ # So rhoN[-1] is reflectivity for no layers but interface from
+ # last layer to substrate
+ # rhoN[0] is total complex reflectivity of the coating stack.
+ rhoN = np.zeros_like(Refl, np.complex128)
+
+ phiNmkm1 = np.flip(Phi) # phi_{N-k-1}
+ rNmkm1 = np.flip(Refl[:-1]) # r_{N-k-1}
+ exp2iphiNmkm1 = np.exp(2j*phiNmkm1) # exp(2i phi_{N-k-1})
+
+ # Recursion relation for complex reflectivity
+ # See https://dcc.ligo.org/T2000552 for derivation
+ rhoN[0] = Refl[-1]
+ for k in range(len(Refl)-1):
+ rhoN[k+1] = ((rNmkm1[k] + exp2iphiNmkm1[k] * rhoN[k])
+ / (1 + exp2iphiNmkm1[k] * rNmkm1[k] * rhoN[k]))
+
+ denTerm = (1 + exp2iphiNmkm1 * rNmkm1 * rhoN[:-1])**2
+
+ # Derivatives of rho_{k+1} wrt to rho_{k}, r_{N-k-1} and phi_{N-k-1}
+ delRhokp1_delRhok = exp2iphiNmkm1 * (1 - rNmkm1**2) / denTerm
+ delRhokp1_delRNmkm1 = np.append(1, ((1 - (exp2iphiNmkm1*rhoN[:-1])**2)
+ / denTerm))
+ delRhokp1_delPhiNmkm1 = np.append(0, -2j * rhoN[:-1] * delRhokp1_delRhok)
+
+ # Derivative of rho_{N} wrt to rho_{N-j}
+ delRhoN_delRhoNmj = np.append(1, np.cumprod(np.flip(delRhokp1_delRhok)))
+
+ # Derivative of rho_{N} wrt to r_k and phi_k
+ delRho_delRk = - delRhoN_delRhoNmj * np.flip(delRhokp1_delRNmkm1)
+ delRho_delPhik = - delRhoN_delRhoNmj * np.flip(delRhokp1_delPhiNmkm1)
+ delLogRho_delReflk = delRho_delRk / rhoN[-1]
+ delLogRho_delPhik = delRho_delPhik / rhoN[-1]
+ delLogRho_delPhik[-1] = 0 # Define this as per Eq (26)
+
+ return rhoN[-1], delLogRho_delPhik, delLogRho_delReflk, Refl
+
+
+def interpretLossAngles(coat):
+ '''
+ Helper function for coating_brownian().
+
+ Creates function from 100 Hz value of loss angle and is logarithmic
+ slope.
+
+ if separate bulk and shear loss angles are not provided as
+ lossBhighn, lossShighn, lossBlown and lossSlown
+ then earlier version names of Phihighn and Philown are searched for and
+ used as Bulk loss angle while setting Shear loss angles to zero.
+
+ Input argument:
+ coat = Coating object containing loss angle values or expressions.
+ Returns:
+ lossBlown = Coating Bulk Loss Angle of Low Refractive Index layer
+ lossSlown = Coating Shear Loss Angle of Low Refractive Index layer
+ lossBhighn = Coating Bulk Loss Angle of High Refractive Index layer
+ lossShighn = Coating Shear Loss Angle of High Refractive Index layer
+ '''
+ if 'lossBhighn' in coat and 'lossShighn' in coat:
+ if 'lossBhighn_slope' in coat:
+ def lossBhighn(f):
+ return coat.lossBhighn * (f / 100)**coat.lossBhighn_slope
+ else:
+ def lossBhighn(f): return coat.lossBhighn
+ if 'lossShighn_slope' in coat:
+ def lossShighn(f):
+ return coat.lossShighn * (f / 100)**coat.lossShighn_slope
+ else:
+ def lossShighn(f): return coat.lossShighn
+ else:
+ # Use Phihighn if specific Bulk & Shear loss angles not provided
+ if 'Phihighn_slope' in coat:
+ def Phihighn(f):
+ return coat.Phihighn * (f / 100)**coat.Phihighn_slope
+ lossBhighn = lossShighn = Phihighn
+ else:
+ lossBhighn = lossShighn = lambda f: coat.Phihighn
+
+ if 'lossBlown' in coat and 'lossSlown' in coat:
+ if 'lossBlown_slope' in coat:
+ def lossBlown(f):
+ return coat.lossBlown * (f / 100)**coat.lossBlown_slope
+ else:
+ def lossBlown(f): return coat.lossBlown
+ if 'lossSlown_slope' in coat:
+ def lossSlown(f):
+ return coat.lossSlown * (f / 100)**coat.lossSlown_slope
+ else:
+ def lossSlown(f): return coat.lossSlown
+ else:
+ # Use Philown if specific Bulk & Shear loss angles not provided
+ if 'Philown_slope' in coat:
+ def Philown(f):
+ return coat.Philown * (f / 100)**coat.Philown_slope
+ lossBlown = lossSlown = Philown
+ else:
+ lossBlown = lossSlown = lambda f: coat.Philown
+
+ return lossBhighn, lossShighn, lossBlown, lossSlown
diff --git a/gwinc/noise/quantum.py b/gwinc/noise/quantum.py
index e474ff30eeb8c08ea0556a01aad8dc2863dc953a..df5f54f8bf77e318dc84969816e8d41543f86867 100644
--- a/gwinc/noise/quantum.py
+++ b/gwinc/noise/quantum.py
@@ -93,7 +93,8 @@ def shotrad(f, ifo):
"""Quantum noise noise strain spectrum
:f: frequency array in Hz
- :ifo: gwinc IFO structure
+ :ifo: gwinc IFO Struct
+ :power: gwinc power Struct
:returns: displacement noise power spectrum at :f:
@@ -327,7 +328,7 @@ def shotradSignalRecycled(f, ifo):
L = ifo.Infrastructure.Length # Length of arm cavities [m]
l = ifo.Optics.SRM.CavityLength # SRC Length [m]
T = ifo.Optics.ITM.Transmittance # ITM Transmittance [Power]
- m = ifo.Materials.MirrorMass # Mirror mass [kg]
+ m = mirror_mass # Mirror mass [kg]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bsloss = ifo.Optics.BSLoss # BS Loss [Power]
@@ -349,7 +350,7 @@ def shotradSignalRecycled(f, ifo):
R = 1 - T - ifo.Optics.Loss # ITM Reflectivity [Power]
- P = ifo.gwinc.parm # use precomputed value
+ P = power.parm # use precomputed value
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -476,7 +477,7 @@ def shotradSignalRecycled(f, ifo):
return coeff, Mifo, Msig, Mnoise
-def shotradSignalRecycledBnC(f, ifo):
+def shotradSignalRecycledBnC(f, ifo, mirror_mass, power):
"""Quantum noise model for signal recycled IFO
See shotrad for more info.
@@ -508,7 +509,7 @@ def shotradSignalRecycledBnC(f, ifo):
L = ifo.Infrastructure.Length # Length of arm cavities [m]
l = ifo.Optics.SRM.CavityLength # SRC Length [m]
T = ifo.Optics.ITM.Transmittance # ITM Transmittance [Power]
- m = ifo.Materials.MirrorMass # Mirror mass [kg]
+ m = mirror_mass # Mirror mass [kg]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
bsloss = ifo.Optics.BSLoss # BS Loss [Power]
@@ -530,7 +531,7 @@ def shotradSignalRecycledBnC(f, ifo):
gamma_ac = T*c/(4*L) # [KLMTV-PRD2001] Arm cavity half bandwidth [1/s]
epsilon = lambda_arm/(2*gamma_ac*L/c) # [BnC, after 5.2] Loss coefficent for arm cavity
- I_0 = ifo.gwinc.pbs # [BnC, Table 1] Power at BS (Power*prfactor) [W]
+ I_0 = power.pbs # [BnC, Table 1] Power at BS (Power*prfactor) [W]
I_SQL = (m*L**2*gamma_ac**4)/(4*omega_0) # [BnC, 2.14] Power to reach free mass SQL
Kappa = 2*((I_0/I_SQL)*gamma_ac**4)/ \
(Omega**2*(gamma_ac**2+Omega**2)) # [BnC 2.13] Effective Radiation Pressure Coupling
diff --git a/gwinc/noise/seismic.py b/gwinc/noise/seismic.py
index ac856f713b1f183a1a83e1f9082443779e61f913..fe8f0f720d78a9cdf059f10357c3ce74dd9564ec 100644
--- a/gwinc/noise/seismic.py
+++ b/gwinc/noise/seismic.py
@@ -6,20 +6,21 @@ import numpy as np
from scipy.interpolate import PchipInterpolator as interp1d
-def seismic_suspension_fitered(sus, in_trans):
+def seismic_suspension_fitered(sus, sustf, in_trans):
"""Seismic displacement noise for single suspended test mass.
- :sus: gwinc suspension structure
+ :sus: suspension Struct
+ :sustf: sus transfer function Struct
:in_trans: input translational displacement spectrum
:returns: tuple of displacement noise power spectrum at :f:, and
horizontal and vertical components.
"""
- hTable = sus.hTable
- vTable = sus.vTable
+ hTable = sustf.hTable
+ vTable = sustf.vTable
- theta = sus.VHCoupling.theta
+ theta = sustf.VHCoupling.theta
# horizontal noise total
nh = (abs(hTable)**2) * in_trans**2
diff --git a/gwinc/noise/substratethermal.py b/gwinc/noise/substratethermal.py
index 0b9e86ebaf8931da1e8888a511457bc89a027b6f..a263eb7ae3c783127a4130b80316568cb9a77b08 100644
--- a/gwinc/noise/substratethermal.py
+++ b/gwinc/noise/substratethermal.py
@@ -12,52 +12,6 @@ from ..const import BESSEL_ZEROS as zeta
from ..const import J0M as j0m
-def substrate_carrierdensity(f, materials, wBeam, exact=False):
- """Substrate thermal displacement noise spectrum from charge carrier density fluctuations
-
- For semiconductor substrates.
-
- :f: frequency array in Hz
- :materials: gwinc optic materials structure
- :wBeam: beam radius (at 1 / e^2 power)
- :exact: whether to use adiabatic approximation or exact calculation (False)
-
- :returns: displacement noise power spectrum at :f:, in meters
-
- """
- H = materials.MassThickness
- diffElec = materials.Substrate.ElectronDiffusion
- diffHole = materials.Substrate.HoleDiffusion
- mElec = materials.Substrate.ElectronEffMass
- mHole = materials.Substrate.HoleEffMass
- cdDens = materials.Substrate.CarrierDensity
- gammaElec = materials.Substrate.ElectronIndexGamma
- gammaHole = materials.Substrate.HoleIndexGamma
- r0 = wBeam/np.sqrt(2)
- omega = 2*pi*f
-
- if exact:
- def integrand(k, om, D):
- return D * k**3 * exp(-k**2 * wBeam**2/4) / (D**2 * k**4 + om**2)
-
- integralElec = np.array([scipy.integrate.quad(lambda k: integrand(k, om, diffElec), 0, inf)[0] for om in omega])
- integralHole = np.array([scipy.integrate.quad(lambda k: integrand(k, om, diffHole), 0, inf)[0] for om in omega])
-
- # From P1400084 Heinert et al. Eq. 15
- #psdCD = @(gamma,m,int) 2*(3/pi^7)^(1/3)*kBT*H*gamma^2*m/hbar^2*cdDens^(1/3)*int; %units are meters
- # FIXME: why the unused argument here?
- def psdCD(gamma, m, int_):
- return 2/pi * H * gamma**2 * cdDens * int_
-
- psdElec = psdCD(gammaElec, mElec, integralElec)
- psdHole = psdCD(gammaHole, mHole, integralHole)
- else:
- psdElec = 4*H*gammaElec**2*cdDens*diffElec/(pi*r0**4*omega**2)
- psdHole = 4*H*gammaHole**2*cdDens*diffHole/(pi*r0**4*omega**2)
-
- return psdElec + psdHole
-
-
def substrate_thermorefractive(f, materials, wBeam, exact=False):
"""Substrate thermal displacement noise spectrum from thermorefractive fluctuations
@@ -87,7 +41,7 @@ def substrate_thermorefractive(f, materials, wBeam, exact=False):
# From P1400084 Heinert et al. Eq. 15
#psdCD = @(gamma,m,int) 2*(3/pi^7)^(1/3)*kBT*H*gamma^2*m/hbar^2*cdDens^(1/3)*int; %units are meters
- psdTR = lambda int_: 2/pi * H * beta**2 * kBT * Temp / (rho*C) * int_;
+ psdTR = lambda int_: 2/pi * H * beta**2 * kBT * Temp / (rho*C) * int_
psd = psdTR(inte)
psd = 2/pi * H * beta**2 * kBT * Temp / (rho*C) * inte
diff --git a/gwinc/noise/suspensionthermal.py b/gwinc/noise/suspensionthermal.py
index 984c6bdf5bdc46ad4a0d441832cfb10a5562c572..f144a97311e09e1a57fe2da5807bfc8143ff8a91 100644
--- a/gwinc/noise/suspensionthermal.py
+++ b/gwinc/noise/suspensionthermal.py
@@ -9,11 +9,12 @@ from ..const import kB
from ..suspension import getJointParams
-def suspension_thermal(f, sus):
+def suspension_thermal(f, sus, sustf):
"""Suspension thermal noise.
:f: frequency array in Hz
:sus: gwinc suspension structure
+ :sustf: sus transfer function Struct
:returns: displacement noise power spectrum at :f:
@@ -25,23 +26,14 @@ def suspension_thermal(f, sus):
pre-calculated and populated into the relevant `sus` fields.
"""
-
+ # suspension TFs
+ hForce = sustf.hForce
+ vForce = sustf.vForce
# and vertical to beamline coupling angle
- theta = sus.VHCoupling.theta
+ theta = sustf.VHCoupling.theta
noise = np.zeros_like(f)
- ##########################################################
- # Suspension TFs
- ##########################################################
-
- hForce = sus.hForce
- vForce = sus.vForce
-
- ##########################################################
- # Thermal Noise Calculation
- ##########################################################
-
# Here the suspension stage list is reversed so that the last stage
# in the list is the test mass
stages = sus.Stage[::-1]
diff --git a/gwinc/plot.py b/gwinc/plot.py
index 3f296ab550bede1268378f22ef74276d5f4952d7..1753108d36102db91ee83d897b03853754e27fa3 100644
--- a/gwinc/plot.py
+++ b/gwinc/plot.py
@@ -1,14 +1,9 @@
-from numpy import sqrt
-from collections.abc import Mapping
-
-
-def plot_noise(
- freq,
- traces,
+def plot_budget(
+ budget,
ax=None,
**kwargs
):
- """Plot a GWINC noise budget from calculated noises
+ """Plot a GWINC BudgetTrace noise budget from calculated noises
If an axis handle is provided it will be used for the plot.
@@ -22,36 +17,29 @@ def plot_noise(
else:
fig = ax.figure
- ylim = kwargs.get('ylim')
-
- for name, trace in traces.items():
- if isinstance(trace, Mapping):
- trace = trace['Total']
+ total = budget.asd
+ ylim = [min(total)/10, max(total)]
+ style = dict(
+ color='#000000',
+ alpha=0.6,
+ lw=4,
+ )
+ style.update(getattr(budget, 'style', {}))
+ if 'label' in style:
+ style['label'] = 'Total ' + style['label']
+ else:
+ style['label'] = 'Total'
+ ax.loglog(budget.freq, total, **style)
- try:
- data = trace[0]
- style = dict(**trace[1])
- except TypeError:
- data = trace
- style = {}
- # assuming all data is PSD
- data = sqrt(data)
- if name == 'Total' and not style:
- style = dict(
- color='#000000',
- alpha=0.6,
- lw=4,
- )
- if ylim is None:
- ylim = [min(data)/10, max(data)]
+ for name, trace in budget.items():
+ style = trace.style
if 'label' not in style:
- style['label'] = name
+ style['label'] = budget.name
if 'linewidth' in style:
style['lw'] = style['linewidth']
elif 'lw' not in style:
style['lw'] = 3
-
- ax.loglog(freq, data, **style)
+ ax.loglog(budget.freq, trace.asd, **style)
ax.grid(
True,
@@ -67,14 +55,19 @@ def plot_noise(
)
ax.autoscale(enable=True, axis='y', tight=True)
- if ylim:
- ax.set_ylim(ylim)
- ax.set_xlim(freq[0], freq[-1])
-
+ ax.set_ylim(kwargs.get('ylim', ylim))
+ ax.set_xlim(budget.freq[0], budget.freq[-1])
ax.set_xlabel('Frequency [Hz]')
if 'ylabel' in kwargs:
ax.set_ylabel(kwargs['ylabel'])
if 'title' in kwargs:
- ax.set_title(kwargs['title'])
+ title = kwargs['title']
+ if 'subtitle' in kwargs:
+ title += '\n' + kwargs['subtitle']
+ ax.set_title(title)
return fig
+
+
+# FIXME: deprecate
+plot_noise = plot_budget
diff --git a/gwinc/struct.py b/gwinc/struct.py
index 85e63921616ebb4c95e56b1de8c1005cca043913..2de5efef366cee6ad04ac379eea92500c9db0a36 100644
--- a/gwinc/struct.py
+++ b/gwinc/struct.py
@@ -173,7 +173,9 @@ class Struct(object):
"""
d = {}
- for k,v in self.__dict__.items():
+ for k, v in self.__dict__.items():
+ if k[0] == '_':
+ continue
if isinstance(v, Struct):
d[k] = v.to_dict(array=array)
else:
@@ -222,6 +224,8 @@ class Struct(object):
"""
for k,v in self.__dict__.items():
+ if k[0] == '_':
+ continue
if isinstance(v, type(self)):
for sk,sv in v.walk():
yield k+'.'+sk, sv
@@ -233,6 +237,25 @@ class Struct(object):
except (AttributeError, TypeError):
yield k, v
+ def hash(self, keys=None):
+ """Hash of Struct.walk() data (recursive)
+
+ """
+ def filter_keys(kv):
+ k, v = kv
+ if keys:
+ return k in keys
+ else:
+ return True
+ def map_tuple(kv):
+ k, v = kv
+ if isinstance(v, list):
+ return k, tuple(v)
+ else:
+ return k, v
+ return hash(tuple(sorted(
+ map(map_tuple, filter(filter_keys, self.walk()))
+ )))
def diff(self, other):
"""Return tuple of differences between target IFO.
diff --git a/gwinc/suspension.py b/gwinc/suspension.py
index 984f868b94cdb078d102f04633779ffe3ce94a95..b71cf7b4f5018d80e55e74ef7c4a54a975690d91 100644
--- a/gwinc/suspension.py
+++ b/gwinc/suspension.py
@@ -420,10 +420,8 @@ def suspQuad(f, sus):
# Complex spring constants
khr = np.zeros([len(stages), len(w)])
khi = np.zeros([len(stages), 2, len(w)])
- kh = np.zeros([len(stages), len(w)], dtype=np.complex_)
kvr = np.zeros([len(stages), len(w)])
kvi = np.zeros([len(stages), 2, len(w)])
- kv = np.zeros([len(stages), len(w)], dtype=np.complex_)
for n, stage in enumerate(stages):
@@ -570,25 +568,20 @@ def suspQuad(f, sus):
for m in range(2*len(stages)):
# turn on only the loss of the current joint
- lossy_region_lower = np.zeros((len(stages), 1))
- lossy_region_upper = np.zeros((len(stages), 1))
n = int(m/2) # stage number
- if m % 2:
- lossy_region_upper[n] = 1
- else:
- lossy_region_lower[n] = 1
+ isLower = m % 2
+ stage_selection = np.zeros((len(stages), 1))
+ stage_selection[n] = 1
# horizontal
- # only the imaginary part due to the specified region is used.
- k = khr + 1j*(khi[:,0,:]*lossy_region_upper +
- khi[:,1,:]*lossy_region_lower)
+ # only the imaginary part due to the specified joint is used.
+ k = khr + 1j*khi[:,isLower,:]*stage_selection
# calculate TFs
hForce[m,:] = tst_force_to_tst_displ(k, masses, f)
# vertical
- # only the imaginary part due to the specified stage is used
- k = kvr + 1j*(kvi[:,0,:]*lossy_region_upper +
- kvi[:,1,:]*lossy_region_lower)
+ # only the imaginary part due to the specified joint is used
+ k = kvr + 1j*kvi[:,isLower,:]*stage_selection
# calculate TFs
vForce[m,:] = tst_force_to_tst_displ(k, masses, f)
diff --git a/gwinc/test/__main__.py b/gwinc/test/__main__.py
index 19a702ece823e0b94db3b69e20f86c091c32c123..24a546792c604d54a939484d115c8d607514804a 100644
--- a/gwinc/test/__main__.py
+++ b/gwinc/test/__main__.py
@@ -6,10 +6,8 @@ import logging
import tempfile
import argparse
import subprocess
-import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
-from collections.abc import Mapping
from PyPDF2 import PdfFileReader, PdfFileWriter
from .. import IFOS, load_budget
@@ -120,81 +118,45 @@ def load_cache(path):
return cache
-def walk_traces(traces, root=()):
- """recursively walk through traces
-
- yields (key_tuple, value), where key_tuple is a tuple of nested
- dict keys of the traces dict locating the given value.
-
- """
- for key, val in traces.items():
- r = root+(key,)
- if isinstance(val, Mapping):
- for k, v in walk_traces(val, root=r):
- yield k, v
- continue
- else:
- yield r, val
-
-
-def zip_noises(tracesA, tracesB, skip):
+def zip_noises(traceA, traceB, skip):
"""zip matching noises from traces"""
- for keys, (noiseA, _) in walk_traces(tracesA):
- # keys is a tuple of nested keys for noise
- name = keys[-1]
- # skip keys we'll add at the end
- if name in ['Total', 'int73']:
- continue
+ B_dict = dict(traceB.walk())
+ for name, tA in traceA.walk():
if skip and name in skip:
logging.warning("SKIPPING TEST: '{}'".format(name))
continue
-
try:
- t = tracesB
- for key in keys:
- t = t[key]
- noiseB, style = t
+ tB = B_dict[name]
except (KeyError, TypeError):
logging.warning("MISSING B TRACE: '{}'".format(name))
continue
+ yield name, tA, tB
- yield name, noiseA, noiseB
-
- yield 'Total', tracesA['Total'][0], tracesB['Total'][0]
- if 'int73' in tracesA:
- yield 'int73', tracesA['int73'][0], tracesB['int73'][0]
-
-
-def compare_traces(tracesA, tracesB, tolerance=TOLERANCE, skip=None):
- """Compare two gwinc trace dictionaries
+def compare_traces(traceA, traceB, tolerance=TOLERANCE, skip=None):
+ """Compare two gwinc traces
Noises listed in `skip` will not be compared.
Returns a dictionary of noises that differ fractionally (on a
- point-by-point basis) by more that `tolerance` between `tracesA`
- and `tracesB`.
+ point-by-point basis) by more than `tolerance` between the traces
+ in `traceA` and `traceB`.
"""
#name_width = max([len(n[0][-1]) for n in walk_traces(tracesA)])
name_width = 15
diffs = OrderedDict()
- for name, noiseA, noiseB in zip_noises(tracesA, tracesB, skip):
+ for name, tA, tB in zip_noises(traceA, traceB, skip):
logging.debug("comparing {}...".format(name))
-
- ampA = np.sqrt(noiseA)
- ampB = np.sqrt(noiseB)
+ ampA = tA.asd
+ ampB = tB.asd
diff = ampB - ampA
frac = abs(diff / ampA)
-
if max(frac) < tolerance:
continue
-
logging.warning("EXCESSIVE DIFFERENCE: {:{w}} {:6.1f} ppm".format(
name, max(frac)*1e6, w=name_width))
-
- diffs[name] = (noiseA, noiseB, frac)
-
+ diffs[name] = (ampA, ampB, frac)
return diffs
@@ -202,11 +164,11 @@ def plot_diffs(freq, diffs, styleA, styleB):
spec = (len(diffs)+1, 2)
sharex = None
for i, nname in enumerate(diffs):
- noiseA, noiseB, frac = diffs[nname]
+ ampA, ampB, frac = diffs[nname]
axl = plt.subplot2grid(spec, (i, 0), sharex=None)
- axl.loglog(freq, np.sqrt(noiseA), **styleA)
- axl.loglog(freq, np.sqrt(noiseB), **styleB)
+ axl.loglog(freq, ampA, **styleA)
+ axl.loglog(freq, ampB, **styleB)
axl.grid()
axl.legend(loc='upper right')
axl.set_ylabel(nname)
@@ -334,14 +296,17 @@ gwinc/test/cache/<SHA1>. Old caches are automatically pruned.""",
fail |= True
continue
- freq, traces_ref, attrs = load_hdf5(path)
-
- Budget = load_budget(name)
- traces_cur = Budget(freq).run()
-
- if inspiral_range:
- total_ref = traces_ref['Total'][0]
- total_cur = traces_cur['Total'][0]
+ traces_ref = load_hdf5(path)
+ traces_ref.name = 'Total'
+ freq = traces_ref.freq
+ budget = load_budget(name, freq)
+ traces_cur = budget.run()
+ traces_cur.name = 'Total'
+
+ # FIXME: add this back
+ if False: #inspiral_range:
+ total_ref = traces_ref.psd
+ total_cur = traces_cur.psd
range_func = inspiral_range.range
H = inspiral_range.waveform.CBCWaveform(freq)
fom_ref = range_func(freq, total_ref, H=H)
diff --git a/gwinc/trace.py b/gwinc/trace.py
new file mode 100644
index 0000000000000000000000000000000000000000..175d5e4628d89de2f2379dd948588cff7d3c5946
--- /dev/null
+++ b/gwinc/trace.py
@@ -0,0 +1,79 @@
+import numpy as np
+
+
+class BudgetTrace:
+ """Budget trace data calculated from a Noise or Budget
+
+
+ """
+ def __init__(self, name=None, style=None, freq=None, psd=None, budget=None):
+ self.name = name
+ self.style = style or {}
+ self._freq = freq
+ self._psd = psd
+ self.budget = budget or []
+
+ def __repr__(self):
+ if self.budget:
+ bs = ' [{}]'.format(', '.join([str(b.name) for b in self]))
+ else:
+ bs = ''
+ return '<{} {}{}>'.format(
+ self.__class__.__name__,
+ self.name,
+ bs,
+ )
+
+ @property
+ def freq(self):
+ """trace frequency array, in Hz"""
+ return self._freq
+
+ @property
+ def psd(self):
+ """trace power spectral density array"""
+ return self._psd
+
+ @property
+ def asd(self):
+ """trace amplitude spectral density array"""
+ return np.sqrt(self._psd)
+
+ def len(self):
+ """Length of data"""
+ return len(self._freq)
+
+ def __iter__(self):
+ """iterator of budget traces"""
+ return iter(self.budget)
+
+ @property
+ def _bdict(self):
+ bdict = {trace.name: trace for trace in self.budget}
+ return bdict
+
+ def __getattr__(self, name):
+ return self._bdict[name]
+
+ def __getitem__(self, name):
+ """get budget trace by name"""
+ return self._bdict[name]
+
+ def items(self):
+ """iterator of budget (name, trace) tuples"""
+ return self._bdict.items()
+
+ def keys(self):
+ """iterator of budget trace names"""
+ return self._bdict.keys()
+
+ def values(self):
+ """iterator of budget traces"""
+ return self._bdict.values()
+
+ def walk(self):
+ """walk recursively through all traces"""
+ yield self.name, self
+ for trace in self:
+ for t in trace.walk():
+ yield t