BJD isn't being computed

[rrubenza]

Currently, the pipeline is not computing barycentric julian dates for any of the timestamps saved in FITS headers. The current flow, as I understand it, is:

1. DATE-MID is the geometric mean of shutter stop/start, in UTC
2. GEOMID_UTC is set to DATE-MID for every order
3. GEOMID_BJD = Time(GEOMID_UTC).jd, which is simply making a copy of GEOMID_UTC since that quantity is already in the format JD_UTC
4. PHOTON_BJD is the flux-weighted GEOMID_BJD per order (could @cwang2016 or @bjfultn comment on the code used to do the flux weighting?)
5. CCD1JD is set to PHOTON_BJD at order 0 in green
6. CCD2JD is set to PHOTON_BJD at order 0 in red

Here's what the above keywords look like for a cloudy SoCal exposure

What immediately stands out from the plot above is GEOMID_BJD is not ~8 minutes different than GEOMID_UTC, because there is no barycentric conversion happening in step 3 above. What this should be changed to is

1. Same
2. Same
3. Using barycorrpy, compute GEOMID_BJD according to:

   if star == 'sun':
   GEOMID_BJD, warn, stat = barycorrpy.utc_tdb.JDUTC_to_HJDTDB(GEOMID_UTC,
                                             lat=obs_config[BarycentricCorrectionAlg.LAT],
                                             longi=obs_config[BarycentricCorrectionAlg.LON],
                                             alt=obs_config[BarycentricCorrectionAlg.ALT],
                                             ra=None,
                                             dec=None,
                                             epoch=None,
                                             pmra=None,
                                             pmdec=None,
                                             px=None,
   )
   else:
   GEOMID_BJD, warn, stat = barycorrpy.utc_tdb.JDUTC_to_BJDTDB(GEOMID_UTC,
                                             ra=obs_config[BarycentricCorrectionAlg.RA],
                                             dec=obs_config[BarycentricCorrectionAlg.DEC],
                                             epoch=obs_config[BarycentricCorrectionAlg.EPOCH],
                                             pmra=obs_config[BarycentricCorrectionAlg.PMRA],
                                             pmdec=obs_config[BarycentricCorrectionAlg.PMDEC],
                                             px=obs_config[BarycentricCorrectionAlg.PX],
                                             lat=obs_config[BarycentricCorrectionAlg.LAT],
                                             longi=obs_config[BarycentricCorrectionAlg.LON],
                                             alt=obs_config[BarycentricCorrectionAlg.ALT],
                                             rv=obs_config[BarycentricCorrectionAlg.RV]
   )

   To be precise, HJD is not exactly the same as BJD, so perhaps it should get its own name GEOMID_HJD -- but it is the correct quantity to compute for the Sun, and for simplicity in the pipeline the keyword could be left as GEOMID_BJD. BJD and HJD will differ by +/- 4 seconds at most.

4. To compute the flux-weighted midpoint, I can't find a dedicated function in barycorrpy to do just this calculation. Instead, there is a function that computes the flux-weighted barycentric RV AND the midpoint at the same time. This step would therefore happen in modules/barycentric_correction/src/alg_barycentric_corr.py to replace get_BC_vel like so:

       if star == 'sun':
           bcvel, PHOTONMID_UTC, warn, stat = barycorrpy.exposure_meter_BC_vel(expmeter_jdutc, expmeter_flux,
                               ra=None,
                               dec=None,
                               epoch=None,
                               pmra=None,
                               pmdec=None,
                               px=None,
                               lat=obs_config[BarycentricCorrectionAlg.LAT],
                               longi=obs_config[BarycentricCorrectionAlg.LON],
                               alt=obs_config[BarycentricCorrectionAlg.ALT],
                               SolSystemTarget='Sun',
                               predictive=True, zmeas=0,
                               rv=None,
                               #rv=obs_config[BarycentricCorrectionAlg.RV]
                               )
           return -bc_obj[0][0]
       else:
           bcvel, PHOTONMID_UTC, warn, stat = barycorrpy.exposure_meter_BC_vel(expmeter_jdutc, expmeter_flux,
                           ra=obs_config[BarycentricCorrectionAlg.RA],
                           dec=obs_config[BarycentricCorrectionAlg.DEC],
                           epoch=obs_config[BarycentricCorrectionAlg.EPOCH],
                           pmra=obs_config[BarycentricCorrectionAlg.PMRA],
                           pmdec=obs_config[BarycentricCorrectionAlg.PMDEC],
                           px=obs_config[BarycentricCorrectionAlg.PX],
                           lat=obs_config[BarycentricCorrectionAlg.LAT],
                           longi=obs_config[BarycentricCorrectionAlg.LON],
                           alt=obs_config[BarycentricCorrectionAlg.ALT],
                           rv=obs_config[BarycentricCorrectionAlg.RV])

   ****Note the above function takes in JD_UTC from the exposure meter timeseries and returns the flux weighted midpoint in JD_UTC as well. So, to get PHOTONMID_BJD, we have to do the same calculation as in step 3, but simply replacing GEOMID_UTC with this new PHOTONMID_UTC.

5. CCD1JD should be set to the mean PHOTON_BJD across all orders in green
6. CCD2JD should be set to the mean PHOTON_BJD across all orders in red

Q: Do steps 5 and 6 here need to be weighted using the same order-weights as are used to compute CCD1RV and CCD2RV, or is it correct to take a simple average across all the orders since the relative weighting will be accounted for in the final RV?

When doing the above for the same example solar frame, I get a correct value for HJD that is ~8 minutes different than the times in the headers (light travel time from Sun to Earth)

[rrubenza]

I'll also note that the time format that should be passed into get_BC_vel is JD_UTC. So, the barycentric corrections as they are currently being computed are correct, because they are receiving (not sure which of GEOMID_UTC/GEOMID_BJD/PHOTON_BJD, PHOTON_BJD would be "correct", but all are actually in JD_UTC as described above). In the fixed scheme, this would mean passing PHOTON_UTC to get_BC_vel, but since we can simply change to use exposure_meter_BC_vel we will get both the correct flux-weighted midpoint and flux-weighted BC velocity in one step, so we can just use exposure_meter_BC_vel instead of get_BC_vel

[cwang2016]

regarding

4. PHOTON_BJD is the flux-weighted GEOMID_BJD per order (could @cwang2016 or @bjfultn comment on the code used to do the flux weighting?)

PHOTON_BJD is computed by midpoint_photo_arrival.orderedMidPoint at radial_velocity/src based on the following input,

• date begin and date end which are either from DATE-BEG and ELAPSED in the header or 'Date-Beg' and 'Date-end' columns in the EXPMETER_SCI table of level 0.
• minimum and maximum wavelength of each order
• EXPMETER_SCI table

I think the data returned from above Midpoint photon arrival should be converted from JDUTC to BJDTDB and replace the column of PHOTON_BJD in current BARY_CORR table.