New drembeam correction, stokes V fix

docs/source/_autosummary/nenupy.schedule.constraints.AzimuthCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.AzimuthCnst
+nenupy.schedule.constraints.AzimuthCnst
 =======================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.Constraint.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.Constraint
+nenupy.schedule.constraints.Constraint
 ======================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.Constraints.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.Constraints
+nenupy.schedule.constraints.Constraints
 =======================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.ElevationCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.ElevationCnst
+nenupy.schedule.constraints.ElevationCnst
 =========================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.LocalSiderealTimeCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.LocalSiderealTimeCnst
+nenupy.schedule.constraints.LocalSiderealTimeCnst
 =================================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.LocalTimeCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.LocalTimeCnst
+nenupy.schedule.constraints.LocalTimeCnst
 =========================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.MeridianTransitCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.MeridianTransitCnst
+nenupy.schedule.constraints.MeridianTransitCnst
 ===============================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.NightTimeCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.NightTimeCnst
+nenupy.schedule.constraints.NightTimeCnst
 =========================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.ScheduleConstraint.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.ScheduleConstraint
+nenupy.schedule.constraints.ScheduleConstraint
 ==============================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.SolarProximityCnst.rst

@@ -0,0 +1,35 @@
+nenupy.schedule.constraints.SolarProximityCnst
+==============================================
+
+.. currentmodule:: nenupy.schedule.constraints
+
+.. autoclass:: SolarProximityCnst
+   :members:
+   :show-inheritance:
+   :inherited-members:
+
+
+
+   .. automethod:: __call__
+   .. automethod:: __init__
+
+
+   .. rubric:: Methods
+
+   .. autosummary::
+
+      ~SolarProximityCnst.__init__
+      ~SolarProximityCnst.get_score
+      ~SolarProximityCnst.plot
+
+
+
+
+
+   .. rubric:: Attributes
+
+   .. autosummary::
+
+      ~SolarProximityCnst.weight
+
+

docs/source/_autosummary/nenupy.schedule.constraints.TargetConstraint.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.TargetConstraint
+nenupy.schedule.constraints.TargetConstraint
 ============================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.TimeRangeCnst.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints.TimeRangeCnst
+nenupy.schedule.constraints.TimeRangeCnst
 =========================================
 
 .. currentmodule:: nenupy.schedule.constraints

docs/source/_autosummary/nenupy.schedule.constraints.rst

@@ -1,4 +1,4 @@
-nenupy.schedule.constraints
+nenupy.schedule.constraints
 ===========================
 
 .. automodule:: nenupy.schedule.constraints
@@ -28,6 +28,7 @@ nenupy.schedule.constraints
       MeridianTransitCnst
       NightTimeCnst
       ScheduleConstraint
+      SolarProximityCnst
       TargetConstraint
       TimeRangeCnst
 

docs/source/_autosummary/nenupy.schedule.schedule.Schedule.rst

@@ -23,6 +23,7 @@
       ~Schedule.export
       ~Schedule.extend_scheduled_observations
       ~Schedule.fine_tune
+      ~Schedule.from_ics
       ~Schedule.insert
       ~Schedule.match_booking
       ~Schedule.plot

nenupy/__init__.py

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

nenupy/astro/beam_correction.py

@@ -3,6 +3,7 @@
 
 __all__ = [
    "compute_jones_matrices",
+   "compute_projection_correction",
    "convert_to_mueller",
    "matrices_to_hdf5"
 ]
@@ -12,6 +13,7 @@
 from astropy.coordinates import SkyCoord
 from astropy.time import Time, TimeDelta
 import astropy.units as u
+import numpy as np
 import h5py
 from typing import Union, Tuple
 import logging
@@ -20,6 +22,8 @@
 
 try:
     from dreambeam.rime.scenarios import on_pointing_axis_tracking
+    from dreambeam.rime.jones import DualPolFieldPointSrc, PJones
+    from dreambeam.telescopes.rt import load_mountedfeed
 except ModuleNotFoundError:
     # This will raise an error eventually with an appropriate message
     pass
@@ -65,6 +69,78 @@ def compute_jones_matrices(
     return Time(times, format="datetime"), frequencies*u.Hz, JonesMatrix(Jn)
 # ============================================================= #
 
+# ============================================================= #
+# -------------- compute_projection_corrections --------------- #
+def compute_projection_corrections(
+        start_time: Time,
+        time_step: TimeDelta,
+        duration: TimeDelta,
+        skycoord: SkyCoord,
+        parallactic: bool = True
+    ) -> Tuple[Time, u.Quantity, JonesMatrix]:
+    """_summary_
+    Took apart on_pointing_axis_tracking method from Dreambeam to only take into account
+    parallactic angle and projection effects.
+
+    Parameters
+    ----------
+    start_time : Time
+        _description_
+    time_step : TimeDelta
+        _description_
+    duration : TimeDelta
+        _description_
+    skycoord : SkyCoord
+        _description_
+    parallactic : bool, optional
+        _description_, by default True
+
+    Returns
+    -------
+    Tuple[Time, u.Quantity, JonesMatrix]
+        _description_
+
+    Raises
+    ------
+    ValueError
+        _description_
+    """
+    log.info("\tComputing Jones projection matrices using DreamBeam...")
+
+    if time_step.sec <= 1.:
+        raise ValueError("DreamBeam does not allow for time intervals lesser than 1 sec.")
+
+    try:
+        stnfeed = load_mountedfeed(
+            tscopename="NenuFAR",
+            station="NenuFAR",
+            band="LBA",
+            modelname="Hamaker-NEC4_Charrier_v1r1"
+        )
+        stnrot = stnfeed.stnRot
+        freqs = stnfeed.getfreqs() # list
+
+        pointingdir = (skycoord.ra.rad, skycoord.dec.rad, "J2000")
+        srcfld = DualPolFieldPointSrc(pointingdir)
+
+        timespy = []
+        nrTimSamps = int((duration.datetime.total_seconds() / time_step.datetime.seconds)) + 1
+        for ti in range(0, nrTimSamps):
+            timespy.append(start_time.datetime + ti * time_step.datetime)
+        pjones = PJones(timespy, np.transpose(stnrot), do_parallactic_rot=parallactic)
+        pjonesOfSrc = pjones.op(srcfld)
+        jones = pjonesOfSrc.getValue()
+        jones = np.repeat(jones[None, ...], len(freqs), axis=0)
+
+    except NameError:
+        log.error(
+            "DreamBeam is not installed. "
+            "See installation instructions https://github.com/2baOrNot2ba/dreamBeam"
+        )
+        raise
+
+    return Time(timespy, format="datetime"), freqs * u.Hz, JonesMatrix(jones)
+# ============================================================= #
 
 # ============================================================= #
 # -------------------- convert_to_mueller --------------------- #

nenupy/io/tf.py

@@ -797,7 +797,8 @@ def info(self, return_str: bool = False) -> str:
             print(message)
 
     def set_default(self) -> None:
-        """Set the default pipeline.
+        """Set the default pipeline. This will reset any existing list of :class:`~nenupy.io.tf.TFTask`.
+        This method also reset the parameters configuration (listed in :attr:`~nenupy.io.tf.TFPipeline.parameters`).
         The list of tasks is:
 
         - Bandpass correction (:meth:`~nenupy.io.tf.TFTask.correct_bandpass`)
@@ -808,6 +809,7 @@ def set_default(self) -> None:
 
         """
 
+        self.parameters.reset()
         self.tasks = [
             TFTask.correct_bandpass(),
             TFTask.remove_channels(),

nenupy/io/tf_utils.py

@@ -33,7 +33,7 @@
 log = logging.getLogger(__name__)
 
 from nenupy.astro import dispersion_delay, faraday_angle
-from nenupy.astro.beam_correction import compute_jones_matrices
+from nenupy.astro.beam_correction import compute_jones_matrices, compute_projection_corrections
 
 __all__ = [
     "blocks_to_tf_data",
@@ -132,7 +132,14 @@ def apply_dreambeam_corrections(
     leftover_time_samples = time_size % time_group_size
 
     # Computing DreamBeam matrices
-    db_time, db_frequency, db_jones = compute_jones_matrices(
+    # db_time, db_frequency, db_jones = compute_jones_matrices(
+    #     start_time=Time(time_unix[0], format="unix", precision=7),
+    #     time_step=TimeDelta(time_group_size * dt_sec, format="sec"),
+    #     duration=TimeDelta(time_unix[-1] - time_unix[0], format="sec"),
+    #     skycoord=skycoord,
+    #     parallactic=parallactic,
+    # )
+    db_time, db_frequency, db_jones = compute_projection_corrections(
         start_time=Time(time_unix[0], format="unix", precision=7),
         time_step=TimeDelta(time_group_size * dt_sec, format="sec"),
         duration=TimeDelta(time_unix[-1] - time_unix[0], format="sec"),
@@ -141,7 +148,7 @@ def apply_dreambeam_corrections(
     )
     db_time = db_time.unix
     db_frequency = db_frequency.to_value(u.Hz)
-    db_jones = np.swapaxes(db_jones, 0, 1)
+    db_jones = np.swapaxes(db_jones, 0, 1) # swap frequency and time axes
 
     # Reshape the data at the time and frequency resolutions
     # Take into account leftover times
@@ -427,7 +434,7 @@ def compute_stokes_parameters(
         I &= \Re(X\overline{X}) + \Re(Y\overline{Y})\\
         Q &= \Re(X\overline{X}) - \Re(Y\overline{Y})\\
         U &= 2\Re(X\overline{Y})\\
-        V &= 2\Im(X\overline{Y})
+        V &= -2\Im(X\overline{Y})
         \end{align}
         
     Parameters
@@ -496,7 +503,7 @@ def compute_stokes_parameters(
         elif stokes_i == "U":
             data_i = data_array[..., 0, 1].real * 2
         elif stokes_i == "V":
-            data_i = data_array[..., 0, 1].imag * 2
+            data_i = - data_array[..., 0, 1].imag * 2 # no negative sign? because data_array[..., 0, 1] = [XrYr + XiYi ; XrYi - XiYr] which is the opposite of XY*=YrXr+YiXi + i(XiYr - XrYi)
         elif stokes_i == "Q/I":
             data_i = (data_array[..., 0, 0].real - data_array[..., 1, 1].real) / (
                 data_array[..., 0, 0].real + data_array[..., 1, 1].real
@@ -509,7 +516,7 @@ def compute_stokes_parameters(
             )
         elif stokes_i == "V/I":
             data_i = (
-                data_array[..., 0, 1].imag
+                - data_array[..., 0, 1].imag
                 * 2
                 / (data_array[..., 0, 0].real + data_array[..., 1, 1].real)
             )
@@ -1845,13 +1852,26 @@ def spectra_data_to_matrix(fft0: da.Array, fft1: da.Array) -> da.Array:
     :class:`~dask.array.Array`
         Reshaped data :math:`\mathbf{d}_{\rm J}`.
 
-    """    
-    row1 = da.stack(
-        [fft0[..., 0], fft1[..., 0] - 1j * fft1[..., 1]], axis=-1  # XX  # XY*
-    )
-    row2 = da.stack(
-        [fft1[..., 0] + 1j * fft1[..., 1], fft0[..., 1]], axis=-1  # YX*  # YY
-    )
+    """
+    # fft0 = [XrXr+XiXi : YrYr+YiYi]
+    # fft1 = [XrYr+XiYi : XrYi-XiYr]
+    xx = fft0[..., 0] # XrXr + XiXi = XX*
+    yy = fft0[..., 1] # YrYr + YiYi = YY*
+    xy = fft1[..., 0] - 1j * fft1[..., 1] # XrYr + XiYi - i(XrYi - XiYr) = XY*
+    yx = fft1[..., 0] + 1j * fft1[..., 1] # XrYr + XiYi + i(XrYi - XiYr) = YX*
+    # row1 = da.stack(
+    #     [fft0[..., 0], fft1[..., 0] - 1j * fft1[..., 1]], axis=-1  # XX  # XY*
+    # )
+    # row2 = da.stack(
+    #     [fft1[..., 0] + 1j * fft1[..., 1], fft0[..., 1]], axis=-1  # YX*  # YY
+    # )
+    # Check using:
+    # xx = np.ones(30).reshape((3, 10))
+    # xy = 2 * np.ones(30).reshape((3, 10))
+    # yx = 3 * np.ones(30).reshape((3, 10))
+    # yy = 4 * np.ones(30).reshape((3, 10))
+    row1 = da.stack([xx, yx], axis=-1)
+    row2 = da.stack([xy, yy], axis=-1)
     return da.stack([row1, row2], axis=-1)