time_delays.py

Go to the documentation of this file.

r"""!
@file  time_delays.py
@brief `main()` shows how to use
[make_time_delay_skymap()](#time_delays.make_time_delay_skymap)
 
* Running the `main` function through `python3 time_delays.py` will print out a DataFrame
  containing one pair of antennas (`Ant 203` and `Ant 104`)
 
* The dataframe will also have a column `time delays [sec]`, which is the "time delay skymap"
  (for this pair of antennas).
 
* [_plot_one_random_row_as_example()](#time_delays._plot_one_random_row_as_example) will plot the
  time delay skymap for a randomly selected pair of antennas.
    * \f$\phi\f$ denotes the azimuthal angle of the incoming signal's direction, and
    * \f$\theta\f$ the zenith angle (in payload coordinates).
    * time (the colorbar) in seconds.
 
![Time Delays of Signals from Different Directions](time_delay_map_in_seconds.svg)
"""
 
import polars as pl
import numpy as np
 

 
 
def _make_grid(sparse=False) -> tuple[np.ndarray, np.ndarray]:
    r"""!
    @ingroup delay
    @brief Creates a 2D grid \f$(\phi, \theta)\f$ of azimuthal and zenith angles.
    @param[in]     sparse see [NumPy's documentation](https://numpy.org/doc/2.2/reference/generated/numpy.meshgrid.html#numpy-meshgrid)
    @retval phi    azimuthal angle coordinates
    @retval theta  zenith angle coordinates
 
    * This grid represents the "bottom hemisphere":
        * Azimuthal angle range:\f$0 < \phi < 2 \pi\f$
        * **Zenith** (not elevation!) angle range: \f$ \pi/2 < \theta < \pi\f$
    * The figure below shows a sparse grid as an example.
    @note The real grid currently in use is 180-by-360
          (i.e. 360 points for \f$\phi\f$, 180 points for \f$\theta\f$)
 
    ![Example 30-by-3 grid](southern_hemisphere.svg)
    """
    azimuthal_angle_mesh = np.linspace(0, 2 * np.pi, 360)
    zenith_angle_mesh = np.linspace(np.pi / 2, np.pi, 180)
    phi, theta = np.meshgrid(azimuthal_angle_mesh, zenith_angle_mesh, sparse=sparse)
 
    return phi, theta
 
 
def _make_unit_vectors_on_bottom_hemisphere() -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    r"""!
    @ingroup delay
    @brief Returns unit vectors pointing to the origin from the surface of a unit bottom hemisphere.
    @retval incmoing_x
    @retval incmoing_y
    @retval incmoing_z
 
    * Calls #_make_grid() to generate a 2D grid of \f$\phi\f$s and \f$\theta\f$s.
    * Converts each point on the grid into Cartesian coordinates; assuming `r=1`:
        * \f$ x = \sin\theta \cdot \cos\phi \f$
        * \f$ y = \sin\theta \cdot \sin\phi \f$
        * \f$ z = \cos\theta                \f$
    * Each point is a unit vector pointing from the origin to somewhere on the hemisphere.
    * Thus, the unit vector \f$\hat{\bf{u}}\equiv(-x,-y,-z)\f$ represents a possible (incoming) signal
        direction, assuming that the payload is at the center of the sphere.
    * The unit bottom hemisphere:
        ```python
        incoming_x, incoming_y, incoming_z = _make_unit_vectors_on_bottom_hemisphere()
        fig = plt.figure()
        ax = fig.add_subplot(projection='3d')
        ax.plot_surface(-incoming_x, -incoming_y, -incoming_z)
        ax.set_aspect('equal')
        ```
 
    ![A Unit "Southern" Hemisphere](southern_hemisphere_3D.svg)
    """
    phi, theta = _make_grid()
    incoming_x = -np.sin(theta) * np.cos(phi)
    incoming_y = -np.sin(theta) * np.sin(phi)
    incoming_z = -np.cos(theta)
 
    return incoming_x, incoming_y, incoming_z
 
 
def make_time_delay_skymap(antenna_pairs: pl.DataFrame) -> pl.DataFrame:
    r"""! Computes the expected incoming-signal time delays between
    [antenna pairs](#antenna_pairs.generate_MI_antenna_pairs)
    based on pairwise displacement vectors and signal directions.
    @ingroup delay
    @param[in]  antenna_pairs The output of #antenna_pairs.generate_MI_antenna_pairs
    @retval time_delays  See sample output below
 
* The following columns are **required** in `antenna_pairs`
        * `A1_X[m]`
        * `A1_Y[m]`
        * `A1_Z[m]`
        * `A2_X[m]`
        * `A2_Y[m]`
        * `A2_Z[m]`
 
* Sample schema of `time_delays`, assuming only the essential columns have been supplied in `antenna_pairs`
```
$ A1_X[m]                              <f64>
$ A1_Y[m]                              <f64>
$ A1_Z[m]                              <f64>
$ A2_X[m]                              <f64>
$ A2_Y[m]                              <f64>
$ A2_Z[m]                              <f64>
$ time delays [sec] <array[f64, (180, 360)]>
```
* Explanation:
    -#   The time delay between any antenna pair is computed via
         \f$ (\bf{X_1} - \bf{X_2}) \cdot \hat{\bf{u}} \f$ / \f$c\f$
         * \f$\bf{X}\f$ denotes the antenna position,
         * \f$\hat{\bf{u}}\f$ is a unit vector that represents the
         [signal direction](#_make_unit_vectors_on_bottom_hemisphere), and
         * \f$c\f$ is the speed of light.
 
    -#  Each entry in the column `time delays [sec]` would be a 2D matrix with the same dimensions
        as the 2D grid generated by #_make_grid; let's call them **time delay skymaps**.
    -#  Each point on a time delay skymap represents the expected time-delay based on where the
        signal is coming from and where the antennas are on the payload.
    -#  Time delay occurs since the signal has to travel an extra distance to reach the second antenna,
        after it hits the first antenna:
        \anchor timedelaybydotproduct
        \image html time_delay_by_dot_product.png width=50%
 
    -#  As an example, the figure below shows the expected time delays (in seconds) between antenna
        203 and 104.
 
    ![Time Delays of Signals from Different Directions](time_delay_map_in_seconds.svg)
    """
    from scipy import constants
    incoming_x, incoming_y, incoming_z = _make_unit_vectors_on_bottom_hemisphere()
 
    # Compute the displacements between antenna pairs
    time_delays = (
        antenna_pairs
        .with_columns(
            pl.concat_arr("A1_X[m]", "A1_Y[m]", "A1_Z[m]").alias("A1_position[m]"),
            pl.concat_arr("A2_X[m]", "A2_Y[m]", "A2_Z[m]").alias("A2_position[m]")
        )
        .with_columns(  # dot product between the antenna displacement vector and signal direction
            (pl.col("A1_position[m]") - pl.col("A2_position[m]")).map_batches(
                lambda displacement_vector: np.array([
                    (dx * incoming_x + dy * incoming_y + dz * incoming_z) / constants.c
                    for (dx, dy, dz) in displacement_vector
                ])
            ).alias("time delays [sec]")
        )
        .select(pl.all().exclude(r"^A[12]_position\[m\]$"))
    )
 
    return time_delays
 
 
def _plot_one_random_row_as_example(one_row: pl.DataFrame, plot_name: str) -> None:
    """!Makes an example plot"""
    import matplotlib.pyplot as plt
    plt.rcParams.update({"font.size": 16})
 
    pl.Config.set_tbl_cols(one_row.width)
    print(one_row)
    phi, theta = _make_grid()
    plt.contourf(phi, theta, one_row["time delays [sec]"][0])
    plt.axis('scaled')
    plt.xlabel(r"$\phi$")
    plt.ylabel(r"$\theta$")
    plt.xticks([0, np.pi, 2 * np.pi], labels=["0", r"$\pi$", r"$2\pi$"])
    plt.yticks([np.pi / 2, np.pi], labels=[r"$\pi / 2$", r"$\pi$"])
    plt.colorbar(label="time delay [sec]")
    plt.savefig(plot_name)
 
 
if __name__ == "__main__":
 
    import os
    from pathlib import Path
    from antenna_attributes import read_MI_antenna_geometry
    from antenna_pairs import generate_MI_antenna_pairs
 
    uncalibrated_antenna_positions: Path = (
        os.environ.get("PUEO_UTIL_INSTALL_DIR") / Path("share/pueo/geometry/jun25/qrh.dat")
    )
    antennas: pl.DataFrame = read_MI_antenna_geometry(uncalibrated_antenna_positions)
 
    antenna_pairs = generate_MI_antenna_pairs(antennas)
 
    try:
        bad = antenna_pairs.select(pl.all().exclude(r"^A[12]_[XYZ]\[m\]$"))
        time_delay = make_time_delay_skymap(bad)
    except pl.exceptions.ColumnNotFoundError:
        print("The above should fail due to missing columns!")
 
    print("But below is okay and the output is:")
    okay = antenna_pairs.select(r"^A[12]_[XYZ]\[m\]$")
    time_delay = make_time_delay_skymap(okay)
 
    _plot_one_random_row_as_example(time_delay.sample(1), plot_name="foobar")