Merge pull request #1609 from LuluAgazie/wavex_model

Addition of WaveX model

AUTHORS.rst

@@ -16,6 +16,7 @@ Contributors
 
 Active developers are indicated by (*). Authors of the PINT paper are indicated by (#).
 
+* Gabriella Agazie (*)
 * Akash Anumarlapudi (*)
 * Anne Archibald (#*)
 * Matteo Bachetti (#)

CHANGELOG-unreleased.md

@@ -9,10 +9,13 @@ the released changes.
 
 ## Unreleased
 ### Changed
+- WAVE parameters can be added to a Wave model with `add_wave_component()` in wave.py 
 - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function.
 - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`).
 ### Added
+- Added WaveX model as DelayComponent with wave amplitudes as fitted parameters
 ### Fixed
+- Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter
 - Fixed RTD by specifying theme explicitly.
 - `.value()` now works for pairParameters
 - Setting `model.PARAM1 = model.PARAM2` no longer overrides the name of `PARAM1`

src/pint/models/__init__.py

@@ -42,6 +42,7 @@
 from pint.models.timing_model import DEFAULT_ORDER, TimingModel
 from pint.models.troposphere_delay import TroposphereDelay
 from pint.models.wave import Wave
+from pint.models.wavex import WaveX
 
 # Define a standard basic model
 StandardTimingModel = TimingModel(

src/pint/models/timing_model.py

@@ -117,6 +117,7 @@
     "spindown",
     "phase_jump",
     "wave",
+    "wavex",
 ]
 
 

src/pint/models/wave.py

@@ -28,7 +28,13 @@ class Wave(PhaseComponent):
 
     def __init__(self):
         super().__init__()
-
+        self.add_param(
+            MJDParameter(
+                name="WAVEEPOCH",
+                description="Reference epoch for wave solution",
+                time_scale="tdb",
+            )
+        )
         self.add_param(
             floatParameter(
                 name="WAVE_OM",
@@ -46,13 +52,6 @@ def __init__(self):
                 parameter_type="pair",
             )
         )
-        self.add_param(
-            MJDParameter(
-                name="WAVEEPOCH",
-                description="Reference epoch for wave solution",
-                time_scale="tdb",
-            )
-        )
         self.phase_funcs_component += [self.wave_phase]
 
     def setup(self):
@@ -64,14 +63,14 @@ def validate(self):
         super().validate()
         self.setup()
         if self.WAVEEPOCH.quantity is None:
-            if self.PEPOCH.quantity is None:
+            if self._parent.PEPOCH.quantity is None:
                 raise MissingParameter(
                     "Wave",
                     "WAVEEPOCH",
                     "WAVEEPOCH or PEPOCH are required if " "WAVE_OM is set.",
                 )
             else:
-                self.WAVEEPOCH = self.PEPOCH
+                self.WAVEEPOCH.quantity = self._parent.PEPOCH.quantity
 
         if (not hasattr(self._parent, "F0")) or (self._parent.F0.quantity is None):
             raise MissingParameter(
@@ -95,6 +94,52 @@ def print_par(self, format="pint"):
 
         return result
 
+    def add_wave_component(self, amps, index=None):
+        """Add Wave Component
+
+        Parameters
+        ----------
+
+        index : int
+            Interger label for Wave components.
+        amps : tuple of float or astropy.quantity.Quantity
+            Sine and cosine amplitudes
+
+        Returns
+        -------
+
+        index :
+            Index that has been assigned to new Wave component
+        """
+        #### If index is None, increment the current max Wave index by 1. Increment using WAVE
+        if index is None:
+            dct = self.get_prefix_mapping_component("WAVE")
+            index = np.max(list(dct.keys())) + 1
+        i = f"{int(index):04d}"
+
+        if int(index) in self.get_prefix_mapping_component("WAVE"):
+            raise ValueError(
+                f"Index '{index}' is already in use in this model. Please choose another"
+            )
+
+        for amp in amps:
+            if isinstance(amp, u.quantity.Quantity):
+                amp = amp.to_value(u.s)
+        self.add_param(
+            prefixParameter(
+                name=f"WAVE{index}",
+                value=amps,
+                units="s",
+                description="Wave components",
+                type_match="pair",
+                long_double=True,
+                parameter_type="pair",
+            )
+        )
+        self.setup()
+        self.validate()
+        return f"{index}"
+
     def wave_phase(self, toas, delays):
         times = 0
         wave_names = ["WAVE%d" % ii for ii in range(1, self.num_wave_terms + 1)]