add additional waveform and detector library examples (gwastro#4456)

* add waveform and detector examples

* add coursel to highlight pycbc features

* fixes suggested by Ian

docs/catalog.rst

@@ -1,3 +1,5 @@
+.. _catalog:
+
 ###################################################
 Catalog of Observed Gravitational-wave Mergers
 ###################################################

docs/conf.py

@@ -38,6 +38,7 @@
               'sphinx.ext.napoleon', 'sphinx.ext.mathjax',
               'matplotlib.sphinxext.plot_directive', 'sphinx.ext.autosummary',
               'sphinx.ext.inheritance_diagram', 'sphinx_design',
+              'sphinx_carousel.carousel',
               "sphinxcontrib.jquery",
               ]
 

docs/detector.rst

@@ -33,4 +33,26 @@ Antenna Patterns and Projecting a Signal into the Detector Frame
 
 .. literalinclude:: ../examples/detector/ant.py
 .. command-output:: python ../examples/detector/ant.py
+
+==============================================================
+Adding a custom detector / overriding existing ones
+==============================================================
+PyCBC supports observatories with arbitrary locations. For the study
+of possible new observatories you can add them explicitly within a script
+or by means of a config file to make the detectors visible to all codes
+that use the PyCBC detector interfaces.
 
+An example of adding a detector directly within a script.
+
+.. plot:: ../examples/detector/custom.py
+   :include-source:
+
+
+The following demonstrates a config file which similarly can provide
+custom observatory information. The options are the same as for the
+direct function calls. To tell PyCBC the location of the config file, 
+set the PYCBC_DETECTOR_CONFIG variable to the location of the file e.g.
+PYCBC_DETECTOR_CONFIG=/some/path/to/detectors.ini. The following would
+provide new detectors 'f1' and 'f2'.
+
+.. literalinclude:: ../examples/detector/custom.ini

docs/index.rst

@@ -16,6 +16,38 @@ gravitational-wave astronomy, please consider `contributing <https://github.com/
 providing feedback, examples, documentation or helping to improve the core
 library and application suite.
 
+.. carousel::
+    :show_controls:
+    :show_indicators:
+    :show_dark:
+    :show_shadows:
+    :show_captions_below:
+
+    .. figure:: https://pycbc.org/pycbc/latest/html/_images/data.png
+        :height: 400px
+        :target: catalog.html
+
+        Working with gravitational wave data
+
+    .. figure:: https://pycbc.org/pycbc/latest/html/_images/plot_detwaveform.png
+        :height: 400px
+        :target: waveform.html
+
+        Your interface to generating gravitational wave signals
+
+    .. figure:: https://collincapano.com/wp-content/uploads/2020/02/posterior3d-gw150914_masses-e1581860222179.png
+        :height: 400px
+        :target: inference.html
+
+        Flexible, easy-to-use, parameter estimation for GW Astronomy
+
+    .. figure:: https://pycbc.org/pycbc/latest/html/_images/demarg_150914.png
+        :height: 400px
+        :target: inference/examples/margtime.html
+
+        Fast Parameter Estimation with Advanced Marginalizations
+
+
 .. toctree::
    :hidden:
    :maxdepth: 1

docs/waveform.rst

@@ -24,6 +24,21 @@ Generating one waveform in multiple detectors
 .. plot:: ../examples/waveform/plot_detwaveform.py
    :include-source:
 
+
+===============================================
+Selecting which modes to include
+===============================================
+Gravitational waves can be decomposed into a set
+of modes. Some approximants only calculate the dominant
+l=2, m=2 mode, while others included higher-order modes. These
+often, but not always, include 'HM' in the name. The modes
+present in the output polarizations can be selected for these waveforms
+as demonstrated below. By default, all modes that the waveform model
+supports are typically generated.
+
+.. plot:: ../examples/waveform/higher_modes.py
+   :include-source:
+
 =======================================
 Calculating the match between waveforms
 =======================================
@@ -51,3 +66,15 @@ Plotting frequency evolution of TD waveform
 
 .. plot:: ../examples/waveform/plot_freq.py
    :include-source:
+
+=====================================
+Adding a custom waveform
+=====================================
+
+You can also add your own custom waveform and make it available through 
+the waveform interface standard. You can directly call the code as below
+or if you include it in a python package, :ref:`PyCBC can directly detect it! <waveform_plugin>`
+
+.. plot:: ../examples/waveform/add_waveform.py
+   :include-source:
+

examples/detector/ant.py

@@ -31,4 +31,4 @@
 # The projection process can also take into account the rotation of the
 # earth using the project wave function.
 hp.start_time = hc.start_time = time
-ht = d.project_wave(hp, hc, ra, dec, pol)
+ht2 = d.project_wave(hp, hc, ra, dec, pol)

examples/detector/custom.ini

@@ -0,0 +1,11 @@
+[detector-f1]
+method = earth_normal
+latitude = 0.7120943348136864
+longitude = -1.9861846887695471
+yangle = 2.0943951023931953
+
+[detector-f2]
+method = earth_normal
+latitude = -0.5497787143782138
+longitude = -2.0594885173533086
+yangle = 2.5830872929516078

examples/detector/custom.py

@@ -0,0 +1,45 @@
+import matplotlib.pyplot as plt
+from pycbc.detector import add_detector_on_earth, Detector
+import pycbc.psd
+import numpy as np
+
+# Set up potential Cosmic Explorer detector locations
+
+# 40 km detector
+lon = -125 / 180.0 * np.pi
+lat = 46 / 180.0 * np.pi
+yangle = 100.0 / 180.0 * np.pi 
+# yangle is the rotation clockwise from pointing north at 0
+# xangle can also be specified and allows for detectors that don't have
+# 90 degree opening between arms. By default we assume xangle is yangle + pi/2
+add_detector_on_earth("C4", lon, lat, yangle=yangle,
+                      xlength=40000, ylength=40000)
+
+# 20 km detector
+# Arm length is optional, but if provided, you can accurately calcuale
+# high-frequency corrects if you provide a frequency argument to the
+# antenna pattern method
+lon = -94 / 180.0 * np.pi
+lat = 29 / 180.0 * np.pi
+yangle = 160.0 / 180.0 * np.pi
+add_detector_on_earth("C2", lon, lat, yangle=yangle,
+                      xlength=20000, ylength=20000)
+
+ra, dec = np.meshgrid(np.arange(0, np.pi*2.0, .1), 
+                      np.arange(-np.pi / 2.0, np.pi / 2.0, .1))
+ra = ra.flatten()
+dec = dec.flatten()
+
+pol = 0
+time = 1e10 + 8000 # A time when ra ~ lines up with lat/lon coordinates
+
+for d in [Detector("C4"), Detector("C2")]:
+    fp, fc = d.antenna_pattern(ra, dec, pol, time)
+
+    plt.figure()
+    plt.subplot(111, projection="mollweide")
+    ra[ra>np.pi] -= np.pi * 2.0
+    plt.scatter(ra, dec, c=fp**2.0 + fc**2.0)
+    plt.title("Mollweide")
+    plt.grid(True)
+    plt.show()

examples/waveform/higher_modes.py

@@ -0,0 +1,42 @@
+import matplotlib.pyplot as pp
+from pycbc.waveform import get_td_waveform
+
+# Let's plot what our new waveform looks like
+pp.figure()
+
+# You can select sets of modes or individual modes using the 'mode_array'
+# The standard format is to provide a list of (l, m) modes, however
+# a string format is also provided to aid use in population from config files.
+# e.g. "22 33" is also acceptable to select these two modes.
+# "None" will result in the waveform return its default which is usually
+# to return all implemented modes.
+for mode_select in [None,  
+                    [(2, 2), (3, 3)], # Select two modes at once
+                    [(2, 2)],
+                    [(2, 1)], 
+                    [(3, 2)],
+                    [(4, 4)],
+                   ]: 
+    hp, hc = get_td_waveform(approximant="IMRPhenomXPHM",
+                         mass1=7,
+                         mass2=40,
+                         f_lower=20.0,
+                         mode_array=mode_select,
+                         inclination = 1.0,
+                         delta_t=1.0/4096)
+
+
+
+    if mode_select is None:
+        label = 'Full Waveform'
+        a = hp.max()
+    else:
+        label = "l, m = " + '  '.join([f"{l}, {m}" for l, m in mode_select])
+
+    (hp / a).plot(label=label)
+
+pp.xlim(-1, 0.05)
+pp.legend()
+pp.xlabel('Time [s]')
+pp.ylabel('Relative Strain')
+pp.show()

requirements.txt

@@ -34,6 +34,7 @@ dynesty
 
 # For building documentation
 Sphinx>=4.2.0
+sphinx-carousel
 sphinx-rtd-theme>=1.0.0
 sphinxcontrib-programoutput>=0.11
 sphinx_design