Merge pull request nanograv#1827 from dlakaplan/freezeswxp

Made SWXP frozen by default, even for first segment

CHANGELOG-unreleased.md

@@ -24,6 +24,7 @@ the released changes.
 - Documentation: Fixed empty descriptions in the timing model components table
 - BIC implementation
 - `event_optimize`: Fixed a bug that was causing the results.txt file to be written without the median values. 
+- SWX model now has SWXP_0001 frozen by default, and new segments should also have SWXP frozen
 ### Removed
 - Removed the argument `--usepickle` in `event_optimize` as the `load_events_weights` function checks the events file type to see if the 
 file is a pickle file.

src/pint/models/solar_wind_dispersion.py

@@ -526,7 +526,8 @@ class SolarWindDispersionX(SolarWindDispersionBase):
     """This class provides a SWX model - multiple Solar Wind segments.
 
     This model lets the user specify time ranges and fit for a different
-    SWXDM (max solar wind DM) value and SWXP (radial power-law index) in each time range.
+    SWXDM (max solar wind DM) value with a different SWXP (radial power-law index)
+    in each time range.
 
     The default radial power-law index value of 2 corresponds to the Edwards et al. model.
     Other values are for the You et al./Hazboun et al. model.
@@ -537,6 +538,8 @@ class SolarWindDispersionX(SolarWindDispersionBase):
     to the standard model: the standard model goes from opposition (min) to conjunction (max),
     while this model goes from 0 to conjunction (max), so the scaling is ``((conjunction - opposition)/conjuction)``.
 
+    Also note that fitting for SWXP is possible but not advised.
+
     See `Solar Wind Examples <examples/solar_wind.html>`_.
 
     To compare against a standard model:
@@ -568,7 +571,7 @@ class SolarWindDispersionX(SolarWindDispersionBase):
     def __init__(self):
         super().__init__()
 
-        self.add_swx_range(None, None, swxdm=0, swxp=2, frozen=False, index=1)
+        self.add_swx_range(None, None, swxdm=0, swxp=2, index=1)
 
         self.set_special_params(["SWXDM_0001", "SWXP_0001", "SWXR1_0001", "SWXR2_0001"])
         self.dm_value_funcs += [self.swx_dm]
@@ -698,7 +701,7 @@ def add_swx_range(
         swxp : float or astropy.quantity.Quantity
             Solar wind power-law index
         frozen : bool
-            Indicates whether SWXDM and SWXP will be fit.
+            Indicates whether SWXDM will be fit.
 
         Returns
         -------
@@ -752,7 +755,7 @@ def add_swx_range(
                 value=swxp,
                 description="Solar wind power-law index",
                 parameter_type="float",
-                frozen=frozen,
+                frozen=True,
                 tcb2tdb_scale_factor=u.Quantity(1),
             )
         )

tests/test_solar_wind.py

@@ -118,6 +118,33 @@ def test_solar_wind_generalmodel_p1():
         toas = make_fake_toas_uniform(54000, 54000 + year, 13, model=model, obs="gbt")
 
 
+def test_swx_frozen():
+    # SWX model with a single segment to match the default model
+    model = get_model(
+        StringIO(
+            "\n".join(
+                [par2, "SWXDM_0001 1\nSWXP_0001 2\nSWXR1_0001 53999\nSWXR2_0001 55000"]
+            )
+        )
+    )
+    assert model.SWXP_0001.frozen is True
+    assert model.SWXDM_0001.frozen is True
+    model = get_model(
+        StringIO(
+            "\n".join(
+                [
+                    par2,
+                    "SWXDM_0001 1\nSWXP_0001 2 1\nSWXR1_0001 53999\nSWXR2_0001 55000",
+                ]
+            )
+        )
+    )
+    assert model.SWXP_0001.frozen is False
+    model.add_swx_range(54000, 54180, swxdm=10, swxp=1.5, frozen=False)
+    assert model.SWXP_0002.frozen is True
+    assert model.SWXDM_0002.frozen is False
+
+
 def test_swx_minmax():
     # default model
     model = get_model(StringIO("\n".join([par2, "NE_SW 1"])))
@@ -266,8 +293,8 @@ def test_swfit_swx():
     model3.SWXR1_0001.value = t.get_mjds().min().value
     model3.SWXR2_0001.value = t.get_mjds().max().value
     model3.SWXP_0001.value = 2
-    model3.SWXP_0001.frozen = True
     model3.SWXDM_0001.value = 1e-3
+    model3.SWXDM_0001.frozen = False
 
     f1 = Fitter.auto(t, model1)
     f3 = Fitter.auto(t, model3)