<jupyter_start><jupyter_text>Waveform TutorialMarch 6, 2019--- Finite Radii WaveformsSpEC datafiles with names `*_FiniteRadii_CodeUnits.h5` hold data extracted at finite radii. When a certain waveform quantity (e.g. $h$, $\Psi_4$, etc.) is computed in SpEC, it is decomposed in terms of the spin-weighted spherical harmonics (SWSHs) ${}_{s}Y_{\ell}^{m}$ with appropriate spin-weight $s$ for the quantity. For example,$$ h(r,t,\theta,\phi) = \sum_{\ell} \sum_{m=-\ell}^{\ell} h^{(\ell,m)}(r,t)\; {}_{-2}Y_{\ell}^{m}(\theta,\phi) \hspace{3cm} (1)$$The mode weights of this decomposition are written to disk as an HDF5 file. So in our example, we would write $h^{(\ell,m)}(r,t)$ to disk for all timesteps $t$ in our simulation and for the specific values of the radius $r$ at which we measured this quantity. Since the strain $h$ and the Weyl scalar $\Psi_4$ fall off predominantly as $r^{-1}$, we multiply the values by the extraction radius before writing to disk, which is why the file names are `rh*.h5` and `rPsi4*.h5`. For finite radius waveforms, all the radius values $r$ are the coordinate radius values. For extrapolation we use the areal radius; see the Extrapolated Waveforms section below for more details.The Python module [h5py](https://www.h5py.org/]) provides convenient tools for working with HDF5 files, which contain an internal directory structure that nicely organizes the data. Let's start by opening up one of these `*_FiniteRadii_CodeUnits.h5` files. To follow this notebook, click [here](https://zenodo.org/record/1215624/files/SXS:BBH:1355/Lev1/rh_FiniteRadii_CodeUnits.h5?download=1) to download a short waveform of the gravitational radiation $rh$ from simulation SXS:BBH:1355 in our catalog. This waveform has the lowest resolution data (Lev1) for an eccentric equal-mass non-spinning binary black hole simulation.<jupyter_code>import h5py
import numpy as np
import matplotlib.pyplot as plt

# Please uncomment the line below and fill in the correct path
rh_file = # /path/to/rh_FiniteRadii_CodeUnits.h5"

# Open the HDF5 file in read-only mode:
rh_data = h5py.File(rh_file,'r')<jupyter_output><empty_output><jupyter_text>The second argument `'r'` that we supplied in `h5py.File` sets the file to read-only mode. The following table lists the other arguments that can be supplied instead:| Argument    | Description                                      ||  :----:     | :----:                                           || `r`         | Readonly, file must exist                        || `r+`	      | Read/write, file must exist                      || `w-` or `x` | Create file, fail if exists                      || `a`	      | Read/write if exists, create otherwise (default) |We navigate through the HDF5 directory structure by specifying the path as a key to `rh_data`. HDF5 directories are known as "groups", so the root group is:`rh_data['/']`If you instead just call:`rh_data`this is assumed to be the root group and is equivalent to the `rh_data['/']` call.<jupyter_code># Let's see what groups we have in the root group:
for group in sorted(rh_data): 
    print(group)
    
# Can also just do:
# > print(sorted(rh_data))<jupyter_output><empty_output><jupyter_text>We should find groups named `R.dir` where the four-digit number `` is the radius $r$ (in simulation coordinate units) at which we computed $h$. We want to look at the largest radius that we can in order to avoid near-zone effects. The largest radius is at $545 M$ in group `R0545.dir`. Let's now see what we have in this group.<jupyter_code>for group in sorted(rh_data['R0545.dir']): 
    print(group)<jupyter_output><empty_output><jupyter_text>We now see that this group is full of datasets. The first few datasets are needed for the extrapolation procedure in [GWFrames](https://github.com/moble/GWFrames), and the rest are files labeled `Y_l*_m*.dat`. These are the datasets that hold the mode weights for $rh$. The spin-weight of the spherical harmonic decomposition is assumed! In other words, these are ***not*** the mode weights for $Y_{\ell}^{m}$ but for ${}_{-2}Y_{\ell}^{m}$. The dominant modes for the strain are the $(\ell=2,m=\pm 2)$ modes so let's take a look at one of those datasets.<jupyter_code>rh_data['R0545.dir/Y_l2_m2.dat'].shape

# For the m=-2 mode:
# rh_data['R0545.dir/Y_l2_m-2.dat'].shape<jupyter_output><empty_output><jupyter_text>The first column of the dataset is the time, the second column is the real part of $rh^{(\ell,m)}$, and the third column is the imaginary part of $rh^{(\ell,m)}$. We see from above that we have 12492 timesteps in this dataset. Let's build a 1D array holding the complex waveform instead of having to access different columns of the HDF5 file:<jupyter_code># Let's take the l=m=2 mode

idx = 'R0545.dir/Y_l2_m2.dat'
rh = rh_data[idx][:,1] + 1j*rh_data[idx][:,2] 

time = rh_data['R0545.dir/Y_l2_m2.dat'][:,0]
ret_time = time - 545 # Define retarded time<jupyter_output><empty_output><jupyter_text>Now we can plot the waveform:<jupyter_code># To plot the waveform interactively, uncomment the following line:
# %matplotlib qt

plt.plot(ret_time, np.real(rh), label='Real')
plt.plot(ret_time, np.imag(rh), label='Imag')
plt.plot(ret_time,  np.abs(rh), label='Abs')
plt.title('Finite Radius Waveform, R=545M')
plt.xlabel('$(t - r)/M$')
plt.ylabel('$rh^{(2,2)}$')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Often, the plot above may be colloquially referred to as a plot of $rh$, but remember that it is actually a plot of just the $\ell=m=2$ mode weight of $rh$!It might also be useful to plot $rh^{(\ell,m)}$ for not just the leading order $(2,2)$ mode but for other modes as well.<jupyter_code>for lm in [(2,2),(4,4),(3,2),(6,6),(5,4),(8,8)]:
    l = str(lm[0])
    m = str(lm[1])
    idx = 'R0545.dir/Y_l'+l+'_m'+m+'.dat'
    ret_time = rh_data[idx][:,0] - 545
    rh = rh_data[idx][:,1] + 1j*rh_data[idx][:,2] 
    
    plt.semilogy(ret_time, np.abs(rh), label='('+l+','+m+')')
    
plt.title('Finite Radius Waveform, R=545M')
plt.xlabel('$(t - r)/M$')
plt.ylabel('$rh^{(\ell,m)}$')
plt.legend()
plt.show()
radii = [int(radius_str[1:5]) for radius_str in sorted(rh_data) if radius_str.startswith('R')]
for R in radii[::2]:
    idx = 'R'+str(R).zfill(4)+'.dir/Y_l2_m2.dat'
    ret_time = rh_data[idx][:,0] - R
    rh = rh_data[idx][:,1] + 1j*rh_data[idx][:,2] 
    
    plt.plot(ret_time, np.abs(rh), label='R='+str(R)+' M')

plt.xlim(2600,3300)
plt.title('Finite Radius Waveform')
plt.xlabel('$(t - r)/M$')
plt.ylabel('$rh^{(2,2)}$')
plt.legend()
plt.show()
# Remember to close the file when you are done
rh_data.close()<jupyter_output><empty_output><jupyter_text>--- Extrapolated WaveformsWe expect that the waveform at asymptotic null infinity has the following fall-off behavior:$$ \lim_{r\rightarrow\infty} h \sim \frac{1}{r} \hspace{3cm} (2)$$However, it is evident that our finite radius waveforms have higher order terms (e.g. $r^{-2}$, $r^{-3}$, etc.) that should not be present in a waveform observed at null infinity. We therefore apply an extrapolation procedure and a center-of-mass correction to our finite radius waveforms in order to have a more reasonable final waveform. In our catalog as of March 6, 2019, only the waveforms ending with `*_Asymptotic_GeometricUnits_CoM.h5` should be used for analysis. If more corrections have been added to the waveforms since then, please use the more up-to-date versions and update this tutorial. These waveforms also inlcude a factor of $M$ as specified in the name of the file in order to render the data dimensionless ($G=c=1$):$$ r M^{-1} h \hspace{1cm} r M \Psi_4 $$For a discussion of the extrapolation procedure, please see:[<NAME> and <NAME>, *Phys. Rev.* **D80**, 124045 (2009)](https://arxiv.org/abs/0905.3177)  For a discussion of other corrections applied to the waveforms, please see:[<NAME>, *Phys. Rev.* **D93**, 084031 (2016)](https://arxiv.org/abs/1509.00862)To follow this part of the notebook, click [here](https://zenodo.org/record/1215624/files/SXS:BBH:1355/Lev1/rhOverM_Asymptotic_GeometricUnits_CoM.h5?download=1) to download an extrapolated waveform. This waveform is the extrapolated version of the finite radius waveform we were looking at above.Note that the finite radius waveforms use coordinate radius values and the extrapolated waveforms use areal radius:$$ r = \sqrt{\frac{1}{4\pi} \int_{S^2} |g|\, d\Omega}$$where $|g|$ is the determinant of the metric, and the surface $S^2$ is the surface of constant coordinate radius where the values of the waveform computed in the simulation.<jupyter_code># Please uncomment the line below and fill in the correct path
rh_file = # "/path/to/rhOverM_Asymptotic_GeometricUnits_CoM.h5"

# Open the HDF5 file in read-only mode:
rh_data = h5py.File(rh_file,'r')
# Let's see what groups we have in the root group:
for group in sorted(rh_data): 
    print(group)
for group in sorted(rh_data['Extrapolated_N4.dir']): 
    print(group)
# Let's take the l=m=2 mode again
idx = 'Extrapolated_N4.dir/Y_l2_m2.dat'
ret_time = rh_data[idx][:,0]
rh = rh_data[idx][:,1] + 1j*rh_data[idx][:,2]

plt.plot(ret_time, np.real(rh), label='Real')
plt.plot(ret_time, np.imag(rh), label='Imag')
plt.plot(ret_time,  np.abs(rh), label='Abs')
plt.title('Extrapolated Waveform, N=4')
plt.xlabel('$(t_{corr} - r^{*})/M$')
plt.ylabel('$rh^{(2,2)}$')
plt.legend()
plt.show()
rh_data.close()<jupyter_output><empty_output><jupyter_text>--- Getting a GW Quantity at a Point on the SphereBefore getting started, you will need the following Python modules:* [quaternion](https://github.com/moble/quaternion)* [spherical_functions](https://github.com/moble/spherical_functions)To compute the actual value of $rh$ and not just the mode weights, we need to perform the summation in Eq. (1), reproduced here for convenience:$$ h(r,t,\theta,\phi) = \sum_{\ell} \sum_{m=-\ell}^{\ell} h^{(\ell,m)}(r,t)\; {}_{-2}Y_{\ell}^{m}(\theta,\phi)$$Unlike the mode weights, the value of $rh$ is a function of $(\theta,\phi)$ so we need to specify a point on the sphere at which to define the SWSHs. Let's define two functions to make our life easier first:<jupyter_code>def swsh(s, lm_list, theta, phi, psi=0):
    import spherical_functions as sf
    import quaternion as qt
    """
    Given a list of modes as (l,m) pairs, returns a list of values of those
    spin-weighted spherical harmonic modes at a point.

    s = spin-weight
    lm_list = [(l,m)] or [(l1,m1),(l2,m2),...]
    theta = polar angle
    phi = azimuthal angle 
    psi = frame orientation
    """
    return sf.SWSH(qt.from_spherical_coords(theta, phi), s, lm_list)*np.exp(1j*s*psi)      


def get_at_point(waveform_modes, lm_list, theta, phi, psi=0, s=-2):
    """
    Computes the value of a gravitational wave quantity at a point on the sphere.
    
    waveform_modes = array where each element is a 1D complex waveform of a single (l,m) mode
    lm_list = a list of all the (l,m) pairs in waveform_modes
    theta = polar angle
    phi = azimuthal angle 
    psi = frame orientation
    s = spin-weight   
    """
    waveform_at_point = np.empty(waveform_modes.shape[1], dtype=complex)
    # For each timestep, compute Eq. (1)
    for t in range(waveform_modes.shape[1]):
        waveform_at_point[t] = np.sum(
            waveform_modes[:, t]*swsh(s, lm_list, theta, phi, psi))
    return waveform_at_point<jupyter_output><empty_output><jupyter_text>Now we can load the extrapolated data from before and pre-compute:* a list of the $(l,m)$ pairs that we want to include in our summation (here we are going to use all of them)* an array where each element is a 1D complex mode weight waveform<jupyter_code>rh_file = # "/path/to/rhOverM_Asymptotic_GeometricUnits_CoM.h5"                                                                 
rh_data = h5py.File(rh_file,'r')

# Get a list of all the data's (l,m) pairs
lm_list = []
mode_strs = sorted(rh_data[sorted(rh_data)[0]])
for mode in mode_strs:
    if mode.startswith('Y'):
        # The l value is the 4th char in the string 'Y_l*_m*.dat'
        # but the m value could have a minus sign so we 
        # make sure we include that if it's present
        if mode[6] == '-':
            lm_list.append((int(mode[3]), int(mode[6:8])))
        else:
            lm_list.append((int(mode[3]), int(mode[6])))

# Alternatively, you could also define your own list of lm pairs 
# if you only wanted the l=2 modes in the summation for example:
# lm_list = [(2,-2),(2,-1),(2,0),(2,1),(2,2)]

waveform_modes = []
for lm_pair in lm_list:
    idx = 'Extrapolated_N4.dir/Y_l'+str(lm_pair[0])+'_m'+str(lm_pair[1])+'.dat'
    waveform_modes.append(rh_data[idx][:,1] + 1j*rh_data[idx][:,2])
waveform_modes = np.array(waveform_modes)

ret_time = rh_data[idx][:,0]<jupyter_output><empty_output><jupyter_text>With all this in hand, we can make a plot of the $rh_+$ waveform at a few points on the sphere:<jupyter_code>for theta in [0,0.25,0.4,0.5]:
    plt.plot(ret_time, np.real(get_at_point(waveform_modes, lm_list, theta*np.pi, 0.0*np.pi)), label='$\theta='+str(theta)+'\pi$')
plt.title('$rh_{+}$ Extrapolated Waveform, N=4')
plt.xlabel('$(t_{corr} - r^{*})/M$')
plt.ylabel('$rh$')
plt.legend()
plt.show()

for theta in [0,0.25,0.4,0.5]:
    plt.plot(ret_time, np.imag(get_at_point(waveform_modes, lm_list, theta*np.pi, 0.0*np.pi)), label='$\theta='+str(theta)+'\pi$')
plt.title('$rh_x$ Extrapolated Waveform, N=4')
plt.xlabel('$(t_{corr} - r^{*})/M$')
plt.ylabel('$rh$')
plt.legend()
plt.show()
rh_data.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13798461083160699, 'fraction_numerical': 0.022935779816513763, 'mean_word_length': 4.292482380579483, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 49, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '338131', 'license': 'MIT', 'n_tokens_mistral': 4883, 'n_tokens_neox': 4540, 'n_words': 1521, 'path': 'Examples/waveform_tutorial.ipynb', 'repo_name': 'SergeiOssokine/catalog_tools', 'size': '13525'}