diff --git a/kadi/commands/observations.py b/kadi/commands/observations.py index fadf3cc3..a96ecb8e 100644 --- a/kadi/commands/observations.py +++ b/kadi/commands/observations.py @@ -109,7 +109,9 @@ def set_star_ids(aca): from kadi.commands import conf q_att = Quat(aca.att) - stars = get_agasc_cone_fast(q_att.ra, q_att.dec, radius=1.2, date=aca.date) + stars = get_agasc_cone_fast( + q_att.ra, q_att.dec, radius=1.2, date=aca.date, matlab_pm_bug=True + ) yang_stars, zang_stars = radec_to_yagzag( stars["RA_PMCORR"], stars["DEC_PMCORR"], q_att ) @@ -127,8 +129,18 @@ def set_star_ids(aca): idx = np.argmin(stars["MAG_ACA"][ok]) aca[idx_aca]["id"] = stars["AGASC_ID"][ok][idx] aca[idx_aca]["mag"] = stars["MAG_ACA"][ok][idx] + # Some debug code from the proper motion stuff. This can be deleted + # in the future. + # if abs(dys[ok][idx]) > 0.5 or abs(dzs[ok][idx]) > 0.5: + # print( + # f'{aca.date}, {stars["AGASC_ID"][ok][idx]}, ' + # f"{aca.obsid}, " + # f"{dys[ok][idx]:.2f}, {dzs[ok][idx]:.2f}, " + # f'{stars["PM_RA"][ok][idx]}, {stars["PM_DEC"][ok][idx]}, ' + # f'{stars["RA"][ok][idx]}, {stars["DEC"][ok][idx]}' + # ) else: - logger.info( + logger.warning( f"WARNING: star idx {idx_aca + 1} not found in obsid {aca.obsid} at " f"{aca.date}" ) @@ -429,8 +441,8 @@ def get_observations( ): """Get observations corresponding to input parameters. - The ``start``, ``stop``, ``starcat_date`` and ``obsid`` parameters serve as matching filters - on the list of observations that is returned. + The ``start``, ``stop``, ``starcat_date`` and ``obsid`` parameters serve as + matching filters on the list of observations that is returned. Over the mission there are thousands of instances of multiple observations with the same obsid, so this function always returns a list of observation @@ -457,9 +469,8 @@ def get_observations( >>> obs_all = get_observations() # All observations in commands archive - # Might be convenient to handle this as a Table - >>> from astropy.table import Table - >>> obs_all = Table(obs_all) + # Might be convenient to handle this as a Table >>> from astropy.table + import Table >>> obs_all = Table(obs_all) >>> from kadi.commands import get_observations >>> get_observations(starcat_date='2022:001:17:00:58.521') @@ -475,20 +486,15 @@ def get_observations( 'starcat_idx': 171677, 'source': 'DEC3021A'}] - :param start: CxoTime-like, None - Start time (default=beginning of commands) - :param stop: CxoTime-like, None - Stop time (default=end of commands) - :param obsid: int, None - ObsID - :param scenario: str, None - Scenario - :param cmds: CommandTable, None - Use this command table instead of querying the archive. - :param starcat_date: CxoTime-like, None - Date of the observation's star catalog - :returns: list of dict - Observation parameters for matching observations. + :param start: CxoTime-like, None Start time (default=beginning of commands) + :param stop: CxoTime-like, None Stop time (default=end of commands) + :param obsid: int, None ObsID + :param scenario: str, None Scenario + :param cmds: CommandTable, None Use this command table instead of querying + the archive. + :param starcat_date: CxoTime-like, None Date of the observation's star + catalog + :returns: list of dict Observation parameters for matching observations. """ from kadi.commands.commands_v2 import get_cmds @@ -517,9 +523,10 @@ def get_observations( else: cmds_obs = cmds[cmds["tlmsid"] == "OBS"] - # Get observations in date range with padding - # _some_ padding is necessary because start/stop are used to filter observations based on - # obs_start, which can be more than 30 minutes after starcat_date. I was generous with padding. + # Get observations in date range with padding. _some_ padding is necessary + # because start/stop are used to filter observations based on obs_start, + # which can be more than 30 minutes after starcat_date. I was generous with + # padding. i0, i1 = cmds_obs.find_date([(start - 7 * u.day).date, (stop + 7 * u.day).date]) cmds_obs = cmds_obs[i0:i1] @@ -550,7 +557,7 @@ def get_observations( return obss -def get_agasc_cone_fast(ra, dec, radius=1.5, date=None): +def get_agasc_cone_fast(ra, dec, radius=1.5, date=None, matlab_pm_bug=False): """ Get AGASC catalog entries within ``radius`` degrees of ``ra``, ``dec``. @@ -560,6 +567,8 @@ def get_agasc_cone_fast(ra, dec, radius=1.5, date=None): :param dec: Declination (deg) :param radius: Cone search radius (deg) :param date: Date for proper motion (default=Now) + :param matlab_pm_bug: bool Apply MATLAB proper motion bug prior to + the MAY2118A loads (default=False) :returns: astropy Table of AGASC entries """ @@ -592,6 +601,19 @@ def get_agasc_cone_fast(ra, dec, radius=1.5, date=None): ok = dists <= radius stars = STARS_AGASC[idx0:idx1][ok] + # Account for a bug in MATLAB proper motion correction that was fixed + # starting with the MAY2118A loads (MATLAB Tools 2018115). The bug was not + # dividing the RA proper motion by cos(dec), so here we premultiply by that + # factor so that add_pmcorr_columns() will match MATLAB. This is purely for + # use in set_star_ids() to match flight catalogs created with MATLAB. + if matlab_pm_bug and CxoTime(date).date < "2018:141:03:35:03.000": + ok = stars["PM_RA"] != -9999 + # Note this is an int16 field so there is some rounding error, but for + # the purpose of star identification this is fine. + stars["PM_RA"][ok] = np.round( + stars["PM_RA"][ok] * np.cos(np.deg2rad(stars["DEC"][ok])) + ) + add_pmcorr_columns(stars, date) return stars