Add pure parallactic correction option

Author: AlanLoh

docs/source/_autosummary/nenupy.astro.beam_correction.compute_projection_corrections.rst

@@ -0,0 +1,6 @@
+nenupy.astro.beam\_correction.compute\_projection\_corrections
+==============================================================
+
+.. currentmodule:: nenupy.astro.beam_correction
+
+.. autofunction:: compute_projection_corrections

docs/source/_autosummary/nenupy.astro.beam_correction.rst

@@ -1,4 +1,4 @@
-nenupy.astro.beam\_correction
+nenupy.astro.beam\_correction
 =============================
 
 .. automodule:: nenupy.astro.beam_correction
@@ -15,6 +15,7 @@ nenupy.astro.beam\_correction
       :toctree:
 
       compute_jones_matrices
+      compute_projection_corrections
       convert_to_mueller
       matrices_to_hdf5
 

docs/source/_autosummary/nenupy.schedule.open_time.NenuEvent.rst

@@ -0,0 +1,29 @@
+nenupy.schedule.open\_time.NenuEvent
+====================================
+
+.. currentmodule:: nenupy.schedule.open_time
+
+.. autoclass:: NenuEvent
+   :members:
+   :show-inheritance:
+   :inherited-members:
+
+   
+   
+   .. automethod:: __init__
+
+   
+   .. rubric:: Methods
+
+   .. autosummary::
+   
+      ~NenuEvent.__init__
+      ~NenuEvent.connect_to_plot_events
+      ~NenuEvent.contains
+      ~NenuEvent.draw
+   
+   
+
+   
+   
+   

docs/source/_autosummary/nenupy.schedule.rst

@@ -33,6 +33,7 @@
    nenupy.schedule.geneticalgo
    nenupy.schedule.obsblocks
    nenupy.schedule.open_time
+   nenupy.schedule.satellite_position
    nenupy.schedule.schedule
    nenupy.schedule.targets
 

nenupy/__init__.py

@@ -5,7 +5,7 @@
 __copyright__ = "Copyright 2023, nenupy"
 __credits__ = ["Alan Loh"]
 __license__ = "MIT"
-__version__ = "2.7.15"
+__version__ = "2.7.16"
 __maintainer__ = "Alan Loh"
 __email__ = "alan.loh@obspm.fr"
 

nenupy/astro/beam_correction.py

@@ -3,7 +3,7 @@
 
 __all__ = [
    "compute_jones_matrices",
-   "compute_projection_correction",
+   "compute_projection_corrections",
    "convert_to_mueller",
    "matrices_to_hdf5"
 ]

nenupy/io/tf.py

@@ -361,6 +361,55 @@ def wrapper_task(
             repeatable=False,
         )
 
+    @classmethod
+    def correct_parallactic_rotation(cls):
+        """:class:`~nenupy.io.tf.TFTask` calling :func:`~nenupy.io.tf_utils.correct_parallactic`.
+        This allows for parallactic angle correction by aplying the inverse rotation.
+
+        Warning
+        -------
+        This task must be computed early in the pipeline process as it involves full Jones operations.
+
+        """
+
+        def wrapper_task(
+            time_unix,
+            frequency_hz,
+            data,
+            dt,
+            channels,
+            dreambeam_dt,
+            dreambeam_skycoord,
+        ):
+            # DreamBeam correction (beam gain + parallactic angle)
+            if (
+                (dreambeam_dt is None)
+                or (dreambeam_skycoord is None)
+            ):
+                return time_unix, frequency_hz, data
+            data = utils.correct_parallactic(
+                time_unix=time_unix,
+                frequency_hz=frequency_hz,
+                data=data,
+                dt_sec=dt.to_value(u.s),
+                time_step_sec=dreambeam_dt.to_value(u.s),
+                n_channels=channels,
+                skycoord=dreambeam_skycoord
+            )
+            return time_unix, frequency_hz, data
+
+        return cls(
+            "Parallactic corection",
+            wrapper_task,
+            [
+                "channels",
+                "dt",
+                "dreambeam_skycoord",
+                "dreambeam_dt",
+            ],
+            repeatable=False,
+        )
+
     @classmethod
     def correct_faraday_rotation(cls):
         """:class:`~nenupy.io.tf.TFTask` calling :func:`~nenupy.io.tf_utils.de_faraday_data` to correct for Faraday rotation for a given ``'rotation_measure'`` set in :attr:`~nenupy.io.tf.TFPipeline.parameters`.

nenupy/io/tf_utils.py

@@ -32,15 +32,17 @@
 
 log = logging.getLogger(__name__)
 
-from nenupy.astro import dispersion_delay, faraday_angle
+from nenupy.astro import dispersion_delay, faraday_angle, parallactic_angle
 from nenupy.astro.beam_correction import compute_jones_matrices, compute_projection_corrections
 
 __all__ = [
+    "apply_dreambeam_corrections",
     "blocks_to_tf_data",
     "compute_spectra_frequencies",
     "compute_spectra_time",
     "compute_stokes_parameters",
     "correct_bandpass",
+    "correct_parallactic",
     "crop_subband_edges",
     "de_disperse_array",
     "de_faraday_data",
@@ -580,6 +582,90 @@ def correct_bandpass(data: np.ndarray, n_channels: int) -> np.ndarray:
     return data.reshape((n_times, n_freqs, 2, 2))
 
 
+# ============================================================= #
+# -------------------- correct_parallactic -------------------- #
+def correct_parallactic(
+        time_unix: np.ndarray,
+        frequency_hz: np.ndarray,
+        data: np.ndarray,
+        dt_sec: float,
+        time_step_sec: float,
+        n_channels: int,
+        skycoord: SkyCoord
+    ) -> np.ndarray:
+
+    # Basic checks to make sure the dimensions are correct
+    freq_size = frequency_hz.size
+    time_size = time_unix.size
+    if time_size != data.shape[0]:
+        raise ValueError("There is a problem in the time dimension!")
+    if (freq_size != data.shape[1]) or (freq_size % n_channels != 0):
+        raise ValueError("There is a problem in the frequency dimension!")
+    n_subbands = int(freq_size / n_channels)
+
+    # Compute the number of time samples that will be corrected together
+    time_group_size = int(np.round(time_step_sec / dt_sec))
+    log.debug(
+        f"\tGroups of {time_group_size} time blocks will be corrected altogether ({dt_sec*time_group_size} sec resolution)."
+    )
+    n_time_groups = time_size // time_group_size
+    leftover_time_samples = time_size % time_group_size
+
+    # Build the time array for the Jones solutions
+    start_time = Time(time_unix[0], format="unix", precision=7)
+    jones_times = start_time + np.arange(time_group_size) * TimeDelta(time_group_size * dt_sec, format="sec")
+    jones_unix = jones_times.unix
+
+    # Reshape the data at the time and frequency resolutions
+    # Take into account leftover times
+    data_leftover = data[-leftover_time_samples:, ...].reshape(
+        (leftover_time_samples, n_subbands, n_channels, 2, 2)
+    )
+    data = data[: time_size - leftover_time_samples, ...].reshape(
+        (n_time_groups, time_group_size, n_subbands, n_channels, 2, 2)
+    )
+
+    # Do the same with the time
+    group_start_time = time_unix[: time_size - leftover_time_samples].reshape(
+        (n_time_groups, time_group_size)
+    )[:, 0]
+    time_start_idx = np.argmax(jones_unix >= group_start_time[0])
+    time_stop_idx = time_start_idx + n_time_groups - 1
+
+    par_angle = parallactic_angle(coordinates=skycoord, time=jones_times).rad
+    jones_parallactic = np.array([
+        [np.cos(par_angle), - np.sin(par_angle)],
+        [np.sin(par_angle), np.cos(par_angle)]
+    ])
+    jones_parallactic = np.swapaxes(jones_parallactic, 2, 0)
+    jones_parallactic = np.linalg.inv(jones_parallactic)
+
+    jones = jones_parallactic[
+        time_start_idx : time_stop_idx + 1, :, :
+    ][:, None, None, None, :, :]
+    jones_leftover = jones_parallactic[
+        -1, :, :
+    ][None, None, None, :, :]
+
+    # Compute the Hermitian matrices
+    jones_transpose = np.swapaxes(jones, -2, -1)
+    jones_leftover_transpose = np.swapaxes(jones_leftover, -2, -1)
+    jones_hermitian = np.conjugate(jones_transpose)
+    jones_leftover_hermitian = np.conjugate(jones_leftover_transpose)
+
+    return np.concatenate(
+        (
+            np.matmul(jones, np.matmul(data, jones_hermitian)).reshape(
+                (time_size - leftover_time_samples, freq_size, 2, 2)
+            ),
+            np.matmul(
+                jones_leftover, np.matmul(data_leftover, jones_leftover_hermitian)
+            ).reshape((leftover_time_samples, freq_size, 2, 2)),
+        ),
+        axis=0
+    )
+
+
 # ============================================================= #
 # -------------------- crop_subband_edges --------------------- #
 def crop_subband_edges(