Make refpol the one that doesn't change in cross_pol_phase_cal, as the docstring describes

hera_cal/abscal.py

@@ -992,6 +992,7 @@ def RFI_delay_slope_cal(reds, antpos, red_data, freqs, rfi_chans, rfi_headings,
         gains[ant] = np.exp(2j * np.pi * np.outer(dlys, freqs))
     return gains
 
+
 def cross_pol_phase_cal(model, data, model_bls, data_bls, wgts={}, return_gains=False, gain_ants=[], refpol="Jee"):
     """
     Solve for the relative phase degeneracy between two polarizations of redundantly calibrated data.
@@ -1009,7 +1010,7 @@ def cross_pol_phase_cal(model, data, model_bls, data_bls, wgts={}, return_gains=
     wgts : dict
         Dictionary mapping baseline tuples to weights. Default is no weights.
     return_gains : bool
-        If True, return gains instead of the relative phase difference between the two polarizations. 
+        If True, return gains instead of the relative phase difference between the two polarizations.
         Default is False.
     gain_ants : list of tuples
         List of antenna-pol tuples to return gains for. Default is None.
@@ -1021,7 +1022,7 @@ def cross_pol_phase_cal(model, data, model_bls, data_bls, wgts={}, return_gains=
     --------
     If return_gains is False:
         delta : np.ndarray
-            Array of relative phase differences between the two polarizations in units of radians. 
+            Array of relative phase differences between the two polarizations in units of radians.
     If return_gains is True:
         delta_gains : dict
             Dictionary mapping antenna keys like (0, 'Jee') to gains of the same shape of the data
@@ -1050,26 +1051,27 @@ def cross_pol_phase_cal(model, data, model_bls, data_bls, wgts={}, return_gains=
 
         if refpol == data_pol1:
             summation += (
-                data[data_bl] * np.conj(model[model_bl]) * wgts.get(data_bl, 1.0)
+                np.conj(data[data_bl]) * model[model_bl] * wgts.get(data_bl, 1.0)
             )
         else:
             summation += (
-                np.conj(data[data_bl]) * model[model_bl] * wgts.get(data_bl, 1.0)
+                data[data_bl] * np.conj(model[model_bl]) * wgts.get(data_bl, 1.0)
             )
 
     if return_gains:
         delta_gains = {}
         for (ant, pol) in gain_ants:
             if pol == refpol:
-                delta_gains[(ant, pol)] = np.exp(1j * np.angle(summation))
-            else:
                 delta_gains[(ant, pol)] = np.ones_like(summation)
+            else:
+                delta_gains[(ant, pol)] = np.exp(1j * np.angle(summation))
 
         return delta_gains
-    
+
     else:
         return np.angle(summation)
 
+
 def dft_phase_slope_solver(xs, ys, data, flags=None):
     '''Solve for spatial phase slopes across an array by looking for the peak in the DFT.
     This is analogous to the method in utils.fft_dly(), except its in 2D and does not

hera_cal/lst_stack/calibration.py

@@ -109,7 +109,7 @@ def _expand_degeneracies_to_ant_gains(
     # Map antenna-polarizations to visibility indices
     antpol_to_idx = {
         utils.split_pol(pol)[0]: polidx
-        for polidx, pol in enumerate(stack.pols) 
+        for polidx, pol in enumerate(stack.pols)
         if utils.split_pol(pol)[0] == utils.split_pol(pol)[1]
     }
 
@@ -157,8 +157,6 @@ def _expand_degeneracies_to_ant_gains(
                         per_day_flags, 1.0 + 0.0j, smooth_ant_gain * rephasor.conj()
                     )
 
-
-
     return gains
 
 
@@ -580,7 +578,7 @@ def lstbin_absolute_calibration(
             )
 
     if run_cross_pol_phase_cal and not all_nights_flagged:
-        cross_pols_in_data  = any(utils.split_pol(pol)[0] != utils.split_pol(pol)[1] for pol in stack.pols)
+        cross_pols_in_data = any(utils.split_pol(pol)[0] != utils.split_pol(pol)[1] for pol in stack.pols)
 
         if cross_pols_in_data:
             delta = _lstbin_cross_pol_phase_calibration(stack=stack, model=model)
@@ -595,7 +593,7 @@ def lstbin_absolute_calibration(
 
             for ant in gain_ants:
                 for pol in unique_pols:
-                    if pol != refpol:
+                    if pol == refpol:
                         phase_gains[(ant, pol)] = phase_gains.get(
                             (ant, pol), np.ones_like(delta)
                         )

hera_cal/lst_stack/tests/test_calibration.py

@@ -306,8 +306,8 @@ def test_relative_phase_calibration(self):
             [
                 np.ones((20, 10, 100)),
                 np.ones((20, 10, 100)),
-                np.exp(1j * delta),
                 np.exp(-1j * delta),
+                np.exp(1j * delta),
             ]
         )
         gains = np.transpose(gains, (1, 2, 3, 0))

hera_cal/tests/test_abscal.py

@@ -755,7 +755,7 @@ def test_complex_phase_abscal(self):
     def test_cross_pol_phase_cal(self):
         rng = np.random.default_rng(42)
         antpos = hex_array(3, split_core=True, outriggers=0)
-        
+
         model, data = {}, {}
         ntimes, nfreqs = 1, 2
 
@@ -775,21 +775,22 @@ def test_cross_pol_phase_cal(self):
             data[bls] = model[bls] * gain
 
         # Solve for the phase degeneracy
-        solved_delta = abscal.cross_pol_phase_cal(model, data, model_bls, data_bls)
-    
+        solved_delta = abscal.cross_pol_phase_cal(model, data, model_bls, data_bls, refpol='Jnn')
+
         # Check that the phase degeneracy was solved for correctly
         np.testing.assert_array_almost_equal(solved_delta, delta, decimal=5)
         gain_ants = set()
         for bl in data_bls:
             gain_ants.update(set(utils.split_bl(bl)))
         gain_ants = list(gain_ants)
-        gains = abscal.cross_pol_phase_cal(model, data, model_bls, data_bls, return_gains=True, gain_ants=gain_ants)
+        gains = abscal.cross_pol_phase_cal(model, data, model_bls, data_bls, return_gains=True, gain_ants=gain_ants, refpol='Jnn')
 
         for k in gains:
             # Check that the gains are correct
             if k[-1] == 'Jnn':
                 np.testing.assert_array_almost_equal(gains[k], 1.0 + 0.0j, decimal=5)
 
+
 @pytest.mark.filterwarnings("ignore:The default for the `center` keyword has changed")
 @pytest.mark.filterwarnings("ignore:invalid value encountered in true_divide")
 @pytest.mark.filterwarnings("ignore:divide by zero encountered in true_divide")