diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index c4ec3337..309bc491 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -2,6 +2,13 @@
 Changelog
 =========
 
+dev
+===
+
+Fixed
+-----
+- API calls for pyuvdata v2.4.0.
+
 v4.1.0 [2023.06.26]
 ===================
 This release heavily focuses on performance of the visibility simulators, and in
diff --git a/docs/conf.py b/docs/conf.py
index 0190d785..bffc1fd0 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -94,7 +94,7 @@
 ]
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = "sphinx"
+pygments_style = "trac"
 
 # -- Options for HTML output -------------------------------------------------
 
diff --git a/hera_sim/adjustment.py b/hera_sim/adjustment.py
index ec8c3644..eab6438b 100644
--- a/hera_sim/adjustment.py
+++ b/hera_sim/adjustment.py
@@ -1,6 +1,5 @@
 """Module providing tools for adjusting simulation data/metadata to a reference."""
 
-import copy
 import logging
 import numpy as np
 import os
@@ -274,7 +273,7 @@ def match_antennas(
     target = _to_uvdata(target)
     if not target.future_array_shapes:  # pragma: nocover
         target.use_future_array_shapes()
-    target_copy = copy.deepcopy(target)
+    target_copy = target.copy()
     reference = _to_uvdata(reference)
     if not reference.future_array_shapes:  # pragma: nocover
         reference.use_future_array_shapes()
@@ -517,6 +516,12 @@ def iswrapped(lsts):
 
     # Ensure reference parameters are a subset of target parameters.
     if axis in ("time", "both"):
+        # Raise an error if the phasing isn't trivial
+        if not target._check_for_cat_type("unprojected").all():
+            raise ValueError(
+                "Time interpolation only supported for unprojected telescopes."
+            )
+
         # Unwrap the LST axis if we have a phase wrap.
         if iswrapped(target_lsts):
             target_lsts[target_lsts < target_lsts[0]] += 2 * np.pi
@@ -546,6 +551,12 @@ def iswrapped(lsts):
         new_baseline_array = np.empty(new_Nblts, dtype=int)
         new_uvw_array = np.empty((new_Nblts, 3), dtype=float)
         new_integration_times = np.empty(new_Nblts, dtype=float)
+        new_phase_center_id_array = np.zeros(new_Nblts, dtype=int)
+        new_phase_center_app_ra = np.empty(new_Nblts, dtype=float)
+        new_phase_center_app_dec = (
+            np.ones(new_Nblts, dtype=float) * target.phase_center_app_dec[0]
+        )
+        new_phase_center_frame_pa = np.zeros(new_Nblts, dtype=float)
         if axis == "both":
             new_data_shape = (new_Nblts, ref_freqs.size, target.Npols)
         else:
@@ -556,6 +567,7 @@ def iswrapped(lsts):
 
     # Actually update metadata and interpolate the data.
     new_data = np.empty(new_data_shape, dtype=complex)
+    history_update = "" if target.history.endswith("\n") else "\n"
     for i, antpair in enumerate(target.get_antpairs()):
         if axis == "freq":
             for pol_ind, pol in enumerate(target.polarization_array):
@@ -571,10 +583,10 @@ def iswrapped(lsts):
         # Preparation for updating metadata.
         ant1, ant2 = antpair
         this_slice = slice(i, None, target.Nbls)
-        old_blt = target._key2inds(antpair)[0][0]  # As a reference
-        this_uvw = target.uvw_array[old_blt]
-        this_baseline = target.baseline_array[old_blt]
-        this_integration_time = target.integration_time[old_blt]
+        old_blts = target._key2inds(antpair)[0]  # As a reference
+        this_uvw = target.uvw_array[old_blts[0]]
+        this_baseline = target.baseline_array[old_blts[0]]
+        this_integration_time = target.integration_time[old_blts[0]]
 
         # Now actually update the metadata.
         new_ant_1_array[this_slice] = ant1
@@ -584,6 +596,10 @@ def iswrapped(lsts):
         new_time_array[this_slice] = ref_times
         new_lst_array[this_slice] = ref_lsts
         new_integration_times[this_slice] = this_integration_time
+        phase_center_interp = interp1d(
+            target_lsts, target.phase_center_app_ra[old_blts], kind="linear"
+        )
+        new_phase_center_app_dec[this_slice] = phase_center_interp(ref_lsts)
 
         # Update the data.
         for pol_ind, pol in enumerate(target.polarization_array):
@@ -617,6 +633,7 @@ def iswrapped(lsts):
     if axis in ("freq", "both"):
         target.Nfreqs = ref_freqs.size
         target.freq_array = ref_freqs
+        history_update += "Data interpolated in frequency with hera_sim.\n"
     if axis in ("time", "both"):
         target.Nblts = ref_times.size * target.Nbls
         target.Ntimes = ref_times.size
@@ -627,9 +644,14 @@ def iswrapped(lsts):
         target.baseline_array = new_baseline_array
         target.uvw_array = new_uvw_array
         target.integration_time = new_integration_times
+        target.phase_center_app_dec = new_phase_center_app_dec
+        target.phase_center_app_ra = new_phase_center_app_ra
+        target.phase_center_id_array = new_phase_center_id_array
+        target.phase_center_frame_pa = new_phase_center_frame_pa
         target.blt_order = None
+        history_update += "Data interpolated in time with hera_sim.\n"
 
-    # Now update the data-like attributes
+    # Now update the data-like attributes; assumes input is unflagged, unavged
     target.flag_array = np.zeros(new_data.shape, dtype=bool)
     target.nsample_array = np.ones(new_data.shape, dtype=float)
     target.data_array = new_data
@@ -789,7 +811,7 @@ def rephase_to_reference(
             vis = data[antpairpol]
             bl = bls[antpairpol]
             new_data[this_slice, :, pol_ind] = lst_rephase(
-                vis, bl, data.freqs, dlst, lat=lat, array=True
+                vis, bl, data.freqs, dlst, lat=lat, inplace=False, array=True
             )
 
     # Convert from HERAData object to UVData object
diff --git a/hera_sim/io.py b/hera_sim/io.py
index 620d491d..79cbf1c0 100644
--- a/hera_sim/io.py
+++ b/hera_sim/io.py
@@ -106,7 +106,7 @@ def empty_uvdata(
     # This is a bit of a hack, but this seems like the only way?
     if pyuvdata.__version__ < "2.2.0":
         uvd.set_drift()
-    else:
+    elif next(iter(uvd.phase_center_catalog.values()))["cat_type"] != "unprojected":
         uvd.fix_phase()
 
     # TODO: the following is a hack patch for pyuvsim which should be fixed there.
diff --git a/hera_sim/sigchain.py b/hera_sim/sigchain.py
index b825e189..42c11098 100644
--- a/hera_sim/sigchain.py
+++ b/hera_sim/sigchain.py
@@ -1086,12 +1086,12 @@ def build_coupling_matrix(
             else:
                 power_beam = uvbeam.copy()
                 power_beam.efield_to_power()
-                if power_beam.interpolation_function is None:
-                    power_beam.interpolation_function = pixel_interp
                 if power_beam.freq_interp_kind is None:
                     power_beam.freq_interp_kind = freq_interp
                 power_beam = power_beam.interp(
-                    freq_array=freqs * units.GHz.to("Hz"), new_object=True
+                    freq_array=freqs * units.GHz.to("Hz"),
+                    new_object=True,
+                    interpolation_function=pixel_interp,
                 )  # Interpolate to the desired frequencies
                 power_beam.to_healpix()
                 power_beam.peak_normalize()
diff --git a/hera_sim/tests/test_adjustment.py b/hera_sim/tests/test_adjustment.py
index b4aef102..3a219d73 100644
--- a/hera_sim/tests/test_adjustment.py
+++ b/hera_sim/tests/test_adjustment.py
@@ -449,6 +449,27 @@ def test_interpolation_in_frequency_with_simulators(base_config, base_sim):
     assert interpolated_sim.data.Nfreqs == base_config["Nfreqs"]
 
 
+def test_interpolation_phased(base_config, base_sim):
+    base_config["Nfreqs"] = 200  # Increase frequency resolution.
+    base_config["start_freq"] = 105e6  # Ensure reference sim freqs contained in target.
+    base_config["bandwidth"] = 40e6
+    ref_sim = Simulator(**base_config)
+
+    # Simulate foregrounds.
+    base_sim.add("diffuse_foreground", Tsky_mdl=Tsky_mdl, omega_p=omega_p)
+    base_sim.data.phase_center_catalog[0]["cat_type"] = "sidereal"
+
+    # Interpolate base_sim to ref_sim along the frequency axis.
+    with pytest.raises(
+        ValueError, match="Time interpolation only supported for unprojected telescopes"
+    ):
+        adjustment.interpolate_to_reference(
+            target=base_sim,
+            reference=ref_sim,
+            axis="time",
+        )
+
+
 def test_interpolation_in_frequency_with_array(base_config, base_sim):
     # Do the same as above, but this time pass a frequency array.
     ref_freqs = np.linspace(105e6, 145e6, 200)  # Hz
diff --git a/hera_sim/tests/test_beams.py b/hera_sim/tests/test_beams.py
index cf1611ca..0a6cac6b 100644
--- a/hera_sim/tests/test_beams.py
+++ b/hera_sim/tests/test_beams.py
@@ -61,7 +61,7 @@ def evaluate_polybeam(polybeam):
     za = np.pi / 2 - np.arcsin(n)
 
     freqs = np.arange(1e8, 2.01e8, 0.04e8)
-    print(len(freqs))
+
     eval_beam = polybeam.interp(az, za, freqs)
 
     # Check that calling the interp() method with wrongly sized
@@ -137,6 +137,7 @@ def run_sim(
             )
             * units.Jy,
             name=["derp"] * flux.shape[1],
+            frame="icrs",
         ),
     )
     simulation = VisibilitySimulation(
diff --git a/hera_sim/tests/test_compare_pyuvsim.py b/hera_sim/tests/test_compare_pyuvsim.py
index 90bbc25d..cf6d5170 100644
--- a/hera_sim/tests/test_compare_pyuvsim.py
+++ b/hera_sim/tests/test_compare_pyuvsim.py
@@ -92,6 +92,7 @@ def get_sky_model(uvdata, nsource):
         spectral_index=sources[:, 3],
         stokes=stokes * units.Jy,
         reference_frequency=Quantity(reference_frequency, "Hz"),
+        frame="icrs",
     )
 
     # Calculate stokes at all the frequencies.
diff --git a/hera_sim/tests/test_vis.py b/hera_sim/tests/test_vis.py
index 695dee28..0a5c2954 100644
--- a/hera_sim/tests/test_vis.py
+++ b/hera_sim/tests/test_vis.py
@@ -169,6 +169,7 @@ def make_point_sky(uvdata, ra: np.ndarray, dec: np.ndarray, align=True):
         name=np.array(["derp"] * len(ra)),
         spectral_type="full",
         freq_array=freqs * units.Hz,
+        frame="icrs",
     )
 
 
@@ -241,6 +242,7 @@ def create_uniform_sky(freq, nbase=4, scale=1) -> SkyModel:
         spectral_type="full",
         freq_array=freq * units.Hz,
         name=np.array([str(i) for i in range(npix)]),
+        frame="icrs",
     )
 
 
diff --git a/hera_sim/visibilities/cli.py b/hera_sim/visibilities/cli.py
index d4677a4d..b239fe91 100644
--- a/hera_sim/visibilities/cli.py
+++ b/hera_sim/visibilities/cli.py
@@ -228,6 +228,8 @@ def run_vis_sim(args):
                     lst_range=None,
                     polarizations=None,
                     blt_inds=None,
+                    phase_center_ids=None,
+                    catalog_names=None,
                 )[0]
 
                 np.save(args.compress, blt_inds)
diff --git a/setup.cfg b/setup.cfg
index 0b9e4e2f..c1b7fbfc 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -63,7 +63,7 @@ docs =
     nbsphinx
     numpydoc>=0.8
     pyradiosky>=0.1.2
-    sphinx>=1.8
+    sphinx>=1.8,<7.2
     sphinx-autorun
     vis-cpu>=1.1.0
 gpu =