Repository: feder-observatory/stellarphot

Issue:
BJD calculation is per star, not per frame
I noticed yesterday that the BJD calculation, which is done in [`PhotometryData`](https://github.com/feder-observatory/stellarphot/blob/main/stellarphot/core.py#L260), has an issue.

In the past, `BJD_TDB` was calculated using the center of the frame as the RA/Dec. Now it is being calculated for each star individually in [`add_bjd_column`](https://github.com/feder-observatory/stellarphot/blob/main/stellarphot/core.py#L260), which results with each star in a single image having a slightly different BJD. 

While at some level that is maybe the correct behavior (or at least not incorrect behavior), since BJD does depend on RA and Dec, it is not usually what one wants. Instead, having a single BJD for each image is more useful.

I have not yet looked into how large a field of view one would need to have for there to be meaningful difference in BJD across the image, but I'm guessing it is tens of degrees.

@JuanCab -- what do you think about this? Maybe meet up to hash out how to handle it?

Discussion:

Patch:
diff --git a/stellarphot/core.py b/stellarphot/core.py
index af29601..e06ffdf 100644
--- a/stellarphot/core.py
+++ b/stellarphot/core.py
@@ -353,11 +353,10 @@ class PhotometryData(BaseEnhancedTable):
     In addition to these required columns, the following columns are created based
     on the input data during creation.
 
-    bjd    (only if ra and dec are all real numbers, otherwise set to np.nan)
     night
 
     If these computed columns already exist in `data` class the class
-    will throw an error a ValueError UNLESS`ignore_computed=True`
+    will throw an error a `ValueError` UNLESS ``ignore_computed=True``
     is passed to the initializer, in which case the columns will be
     retained and not replaced with the computed values.
     """
@@ -485,8 +484,8 @@ class PhotometryData(BaseEnhancedTable):
                         f"{input_data[this_col].unit})."
                     )
 
-            # Compute additional columns (not done yet)
-            computed_columns = ["bjd", "night"]
+            # Compute additional columns
+            computed_columns = ["night"]
 
             # Check if columns exist already, if they do and retain_user_computed is
             # False,  throw an error.
@@ -502,8 +501,6 @@ class PhotometryData(BaseEnhancedTable):
                     # Compute the columns that need to be computed (match requires
                     # python>=3.10)
                     match this_col:
-                        case "bjd":
-                            self["bjd"] = self.add_bjd_col(self.observatory)
 
                         case "night":
                             # Generate integer counter for nights. This should be
@@ -536,7 +533,7 @@ class PhotometryData(BaseEnhancedTable):
                                 name="night",
                             )
 
-                        case _:
+                        case _:  # pragma: no cover
                             raise ValueError(
                                 f"Trying to compute column ({this_col}). "
                                 "This should never happen."
@@ -547,18 +544,36 @@ class PhotometryData(BaseEnhancedTable):
                 self._passband_map = passband_map.model_copy()
                 self._update_passbands()
 
-    def add_bjd_col(self, observatory):
+    def add_bjd_col(self, observatory=None, bjd_coordinates=None):
         """
         Returns a astropy column of barycentric Julian date times corresponding to
         the input observations.  It modifies that table in place.
+
+        Parameters
+        ----------
+        observatory: `stellarphot.settings.Observatory`
+            Information about the observatory. Defaults to the observatory in
+            the table metadata.
+
+        bjd_coordinates: `astropy.coordinates.SkyCoord`, optional (Default: None)
+            The coordinates to use for computing the BJD. If None, the RA and Dec
+            columns in the table will be used.
         """
+        if observatory is None:
+            if self.observatory is None:
+                raise ValueError(
+                    "You must provide an observatory object to compute BJD."
+                )
+            observatory = self.observatory
 
-        if np.isnan(self["ra"]).any() or np.isnan(self["dec"]).any():
+        if bjd_coordinates is None and (
+            np.isnan(self["ra"]).any() or np.isnan(self["dec"]).any()
+        ):
             print(
-                "WARNING: BJD could not be computed in output PhotometryData object "
+                "WARNING: BJD could not be computed "
                 "because some RA or Dec values are missing."
             )
-            return np.full(len(self), np.nan)
+            self["bjd"] = np.full(len(self), np.nan)
         else:
             # Convert times at start of each observation to TDB (Barycentric Dynamical
             # Time)
@@ -567,14 +582,18 @@ class PhotometryData(BaseEnhancedTable):
             times_tdb.format = "jd"  # Switch to JD format
 
             # Compute light travel time corrections
-            ip_peg = SkyCoord(ra=self["ra"], dec=self["dec"], unit="degree")
+            if bjd_coordinates is None:
+                sky_coords = SkyCoord(ra=self["ra"], dec=self["dec"], unit="degree")
+            else:
+                sky_coords = bjd_coordinates
+
             ltt_bary = times.light_travel_time(
-                ip_peg, location=observatory.earth_location
+                sky_coords, location=observatory.earth_location
             )
             time_barycenter = times_tdb + ltt_bary
 
             # Return BJD at midpoint of exposure at each location
-            return Time(time_barycenter + self["exposure"] / 2, scale="tdb")
+            self["bjd"] = Time(time_barycenter + self["exposure"] / 2, scale="tdb")
 
 
 class CatalogData(BaseEnhancedTable):


Test patch:
diff --git a/stellarphot/tests/test_core.py b/stellarphot/tests/test_core.py
index 80964e5..77257b9 100644
--- a/stellarphot/tests/test_core.py
+++ b/stellarphot/tests/test_core.py
@@ -491,8 +491,11 @@ def test_photometry_blank():
     assert test_base.observatory is None
 
 
-def test_photometry_data(feder_cg_16m, feder_passbands, feder_obs):
+@pytest.mark.parametrize("bjd_coordinates", [None, "custom"])
+def test_photometry_data(feder_cg_16m, feder_passbands, feder_obs, bjd_coordinates):
     # Create photometry data instance
+    custom_coord = SkyCoord(30.0, "22:30:20.77733059", unit="degree")
+
     phot_data = PhotometryData(
         observatory=feder_obs,
         camera=feder_cg_16m,
@@ -517,16 +520,33 @@ def test_photometry_data(feder_cg_16m, feder_passbands, feder_obs):
     assert phot_data.observatory.earth_location.height.unit == u.m
     assert phot_data["night"][0] == 59909
 
-    # Checking the BJD computation against Ohio State online calculator for
-    # UTC 2022 11 27 06 27 29.620
-    # Latitude 46.86678
-    # Longitude -96.45328
-    # Elevation 311
-    # RA 78.17278712191924
-    # Dec 22 30 20.77733059
-    # which returned 2459910.775405664 (Uses custom IDL, astropy is SOFA checked).
-    # Demand a difference of less than 1/20 of a second.
-    assert (phot_data["bjd"][0].value - 2459910.775405664) * 86400 < 0.05
+    if bjd_coordinates is None:
+        # BJD calculation is done based on the locations of individual stars
+        phot_data.add_bjd_col()
+        # Checking the BJD computation against Ohio State online calculator for
+        # UTC 2022 11 27 06 27 29.620
+        # Latitude 46.86678
+        # Longitude -96.45328
+        # Elevation 311
+        # RA 78.17278712191924
+        # Dec 22 30 20.77733059
+        # which returned 2459910.775405664 (Uses custom IDL, astropy is SOFA checked).
+        # Demand a difference of less than 1/20 of a second.
+        assert (phot_data["bjd"][0].value - 2459910.775405664) * 86400 < 0.05
+    else:
+        # Do the BJD calculation based on the provided coordinates, same for
+        # all strs in the image.
+        phot_data.add_bjd_col(bjd_coordinates=custom_coord)
+        # Checking the BJD computation against Ohio State online calculator for
+        # UTC 2022 11 27 06 27 29.620
+        # Latitude 46.86678
+        # Longitude -96.45328
+        # Elevation 311
+        # RA 30
+        # Dec 22 30 20.77733059
+        # which returned 2459910.774772048 (Uses custom IDL, astropy is SOFA checked).
+        # Demand a difference of less than 1/20 of a second.
+        assert (phot_data["bjd"][0].value - 2459910.774772048) * 86400 < 0.05
 
 
 def test_photometry_data_short_filter_name(feder_cg_16m, feder_obs):