diff --git a/flystar/align.py b/flystar/align.py index 994a3b1..b87f44d 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -1288,6 +1288,8 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot ye=ye_trans_arr[:,boot_idx], me=me_trans_arr[:,boot_idx], t=np.tile(t_boot, (len(ref_table),1))) + if 'motion_model_used' in ref_table.columns: + star_table['motion_model_input'] = ref_table['motion_model_used'] # Now, do proper motion calculation, making sure to fix t0 to the # orig value (so we can get a reasonable error on x0, y0) diff --git a/flystar/analysis.py b/flystar/analysis.py index 81ab3f4..3121458 100644 --- a/flystar/analysis.py +++ b/flystar/analysis.py @@ -438,7 +438,8 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False): combined astrometry + uncombined photometry table. """ # Multiples: ['x', 'y', 'm', 'name_in_list', 'xe', 'ye', 'me', 't', - # 'x_orig', 'y_orig', 'm_orig', 'xe_orig', 'ye_orig', 'me_orig', 'used_in_trans'] + # 'x_orig', 'y_orig', 'm_orig', 'xe_orig', 'ye_orig', 'me_orig', 'used_in_trans', + # 'xe_boot','ye_boot','me_boot'] # Single: ['name', 'm0', 'm0_err', 'use_in_trans', 'ref_orig', 'n_detect', # 'x0', 'vx', 'y0', 'vy', 'x0_err', 'vx_err', 'y0_err', 'vy_err', 't0'] # Don't include n_vfit @@ -446,8 +447,8 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False): new_tab = copy.deepcopy(tab) #new_tab.remove_column('n_fit') new_tab.remove_column('n_detect') - for col in ['x','y','m','xe','ye','me','t','x_orig','y_orig','m_orig', - 'xe_orig','ye_orig','me_orig','used_in_trans']: + for col in ['x','y','m','name_in_list','xe','ye','me','t','x_orig','y_orig','m_orig', + 'xe_orig','ye_orig','me_orig','used_in_trans','xe_boot','ye_boot','me_boot']: new_tab[col] = tab[col][:,idx] new_tab.combine_lists('m', weights_col='me', sigma=3, ismag=True) diff --git a/flystar/motion_model.py b/flystar/motion_model.py index e24067b..0b86d07 100644 --- a/flystar/motion_model.py +++ b/flystar/motion_model.py @@ -26,27 +26,47 @@ class MotionModel(ABC): optional_param_names = [] def __init__(self, *args, **kwargs): - # TODO: do we need this? - '''for param in self.fitter_param_names: - param_var = getattr(self, param) - if not isinstance(param_var, (list, np.ndarray)): - setattr(self, param, np.array([param_var]))''' + """ + Make a motion model object. This object defines the fitter and fixed parameters, + and if needed stores metadata such as RA and Dec for Parallax, + for the given motion model and contains functions to fit these values to data + and apply the values to compute expected positions at given times. Each instance + corresponds to a given motion model, not an individual star, and thus the fit + values are only input/returned in functions and not stored in the object. + """ return def get_pos_at_time(self, params, t): + """ + Position calculator for a single star using a given motion model and input + model parameters and times. + """ #return x, y pass def get_batch_pos_at_time(self, t): + """ + Position calculator for a set of stars using a given motion model and input + model parameters and times. + """ #return x, y, x_err, y_err pass def run_fit(self, t, x, y, xe, ye, t0, weighting='var', use_scipy=True, absolute_sigma=True): + """ + Run a single fit of the data to the motion model and return the best parameters. + This function is used by the overall fit_motion_model function once for a basic fit + or several times for a bootstrap fit. + """ # Run a single fit (used both for overall fit + bootstrap iterations) pass def get_weights(self, xe, ye, weighting='var'): + """ + Get the weights for each data point for fitting. Options are 'var' (default) + and 'std'. + """ if weighting=='std': return 1./xe, 1./ye elif weighting=='var': @@ -56,6 +76,9 @@ def get_weights(self, xe, ye, weighting='var'): return 1./xe**2, 1./ye**2 def scale_errors(self, errs, weighting='var'): + """ + Rescale the fit result errors as needed, according to the weighting scheme used. + """ if weighting=='std': return np.array(errs)**2 elif weighting=='var': @@ -70,6 +93,7 @@ def fit_motion_model(self, t, x, y, xe, ye, t0, bootstrap=0, weighting='var', Fit the input positions on the sky and errors to determine new parameters for this motion model (MM). Best-fit parameters will be returned along with uncertainties. + Optionally, bootstrap error estimation can be performed. """ params, param_errs = self.run_fit(t, x, y, xe, ye, t0=t0, weighting=weighting, use_scipy=use_scipy, absolute_sigma=absolute_sigma) @@ -99,8 +123,7 @@ def fit_motion_model(self, t, x, y, xe, ye, t0, bootstrap=0, weighting='var', def get_chi2(self, fit_params, fixed_params, t, x, y, xe, ye, reduced=False): """ - Get the chi^2 value for the current MM and - the input data. + Get the chi^2 value for the input motion model parameters and data. """ x_pred, y_pred = self.get_pos_at_time(fit_params, fixed_params, t) chi2x = np.sum((x-x_pred)**2 / xe**2) @@ -345,8 +368,7 @@ class Parallax(MotionModel): """ Motion model for linear proper motion + parallax - Requires RA, Dec, and PA parameters (degrees) for parallax calculation. - RA, Dec in J2000 + Requires RA & Dec (J2000) for parallax calculation. Optional PA is counterclockwise offset of the image y-axis from North. Optional obs parameter describes observer location, default is 'earth'. """ @@ -450,11 +472,13 @@ def fit_func(use_t, x0,vx, y0,vy, pi): param_errors = [x0_err, vx_err, y0_err, vy_err, pi_err] return params, param_errors -""" -Check that everything is set up properly for motion models to run and their -required metadata. -""" + def validate_motion_model_dict(motion_model_dict, startable, default_motion_model): + """ + Check that everything is set up properly for motion models to run and their + required metadata. + """ + # Collect names of all motion models that might get used. all_motion_model_names = ['Fixed'] if default_motion_model is not None: @@ -478,11 +502,12 @@ def validate_motion_model_dict(motion_model_dict, startable, default_motion_mode return motion_model_dict -""" -Get all the motion model parameters for a given motion_model_name. -Optionally, include fixed and error parameters (included by default). -""" + def get_one_motion_model_param_names(motion_model_name, with_errors=True, with_fixed=True): + """ + Get all the motion model parameters for a given motion_model_name. + Optionally, include fixed and error parameters (included by default). + """ mod = eval(motion_model_name) list_of_parameters = [] list_of_parameters += getattr(mod, 'fitter_param_names') @@ -492,11 +517,12 @@ def get_one_motion_model_param_names(motion_model_name, with_errors=True, with_f list_of_parameters += [par+'_err' for par in getattr(mod, 'fitter_param_names')] return list_of_parameters -""" -Get all the motion model parameters for all models given in motion_model_list. -Optionally, include fixed and error parameters (included by default). -""" + def get_list_motion_model_param_names(motion_model_list, with_errors=True, with_fixed=True): + """ + Get all the motion model parameters for all models given in motion_model_list. + Optionally, include fixed and error parameters (included by default). + """ list_of_parameters = [] all_motion_models = [eval(mm) for mm in np.unique(motion_model_list).tolist()] for aa in range(len(all_motion_models)): @@ -512,11 +538,12 @@ def get_list_motion_model_param_names(motion_model_list, with_errors=True, with_ return np.unique(list_of_parameters).tolist() -""" -Get all the motion model parameters for all models defined in this module. -Optionally, include fixed and error parameters (included by default). -""" + def get_all_motion_model_param_names(with_errors=True, with_fixed=True): + """ + Get all the motion model parameters for all models defined in this module. + Optionally, include fixed and error parameters (included by default). + """ list_of_parameters = [] all_motion_models = MotionModel.__subclasses__() for aa in range(len(all_motion_models)): diff --git a/flystar/parallax.py b/flystar/parallax.py index b4b3a1b..4792ec6 100755 --- a/flystar/parallax.py +++ b/flystar/parallax.py @@ -64,39 +64,6 @@ def parallax_in_direction(RA, Dec, mjd, obsLocation='earth', PA=0): return pvec -def dparallax_dt_in_direction(RA, Dec, mjd, obsLocation='earth'): - """ - R.A. in degrees. (J2000) - Dec. in degrees. (J2000) - MJD - - Equations following MulensModel. - Time derivative --> units are yr^-1 - - """ - # print('parallax_in_direction: len(t) = ', len(mjd)) - # Munge inputs into astropy format. - times = Time(mjd + 2400000.5, format='jd', scale='tdb') - coord = SkyCoord(RA, Dec, unit=(units.deg, units.deg)) - - direction = coord.cartesian.xyz.value - north = np.array([0., 0., 1.]) - _east_projected = np.cross(north, direction) / np.linalg.norm(np.cross(north, direction)) - _north_projected = np.cross(direction, _east_projected) / np.linalg.norm(np.cross(direction, _east_projected)) - - obs_posvel = get_observer_barycentric(obsLocation, times, velocity=True)[1] - sun_posvel = get_body_barycentric_posvel('Sun', times)[1] - sun_obs_vel = sun_posvel - obs_posvel - vel = sun_obs_vel.xyz.T.to(units.au / units.year) - - e = np.dot(vel, _east_projected) - n = np.dot(vel, _north_projected) - - dpvec_dt = np.array([e.value, n.value]).T - - return dpvec_dt - - def get_observer_barycentric(body, times, min_ephem_step=1, velocity=False): """ Get the barycentric position of a satellite or other Solar System body @@ -180,152 +147,3 @@ def get_observer_barycentric(body, times, min_ephem_step=1, velocity=False): return obs_pos -def sun_position(mjd, radians=False): - """ - - NAME: - SUNPOS - - PURPOSE: - To compute the RA and Dec of the Sun at a given date. - - INPUTS: - mjd - The modified Julian date of the day (and time), scalar or vector - - OUTPUTS: - ra: - | The right ascension of the sun at that date in DEGREES - | double precision, same number of elements as jd - dec: - The declination of the sun at that date in DEGREES - elong: - Ecliptic longitude of the sun at that date in DEGREES. - obliquity: - the obliquity of the ecliptic, in DEGREES - - OPTIONAL INPUT KEYWORD: - RADIAN [def=False] - If this keyword is set to True, then all output variables - are given in Radians rather than Degrees - - NOTES: - Patrick Wallace (Rutherford Appleton Laboratory, UK) has tested the - accuracy of a C adaptation of the sunpos.pro code and found the - following results. From 1900-2100 SUNPOS gave 7.3 arcsec maximum - error, 2.6 arcsec RMS. Over the shorter interval 1950-2050 the figures - were 6.4 arcsec max, 2.2 arcsec RMS. - - The returned RA and Dec are in the given date's equinox. - - Procedure was extensively revised in May 1996, and the new calling - sequence is incompatible with the old one. - METHOD: - Uses a truncated version of Newcomb's Sun. Adapted from the IDL - routine SUN_POS by CD Pike, which was adapted from a FORTRAN routine - by B. Emerson (RGO). - EXAMPLE: - (1) Find the apparent RA and Dec of the Sun on May 1, 1982 - - | IDL> jdcnv, 1982, 5, 1,0 ,jd ;Find Julian date jd = 2445090.5 - | IDL> sunpos, jd, ra, dec - | IDL> print,adstring(ra,dec,2) - | 02 31 32.61 +14 54 34.9 - - The Astronomical Almanac gives 02 31 32.58 +14 54 34.9 so the error - in SUNPOS for this case is < 0.5". - - (2) Find the apparent RA and Dec of the Sun for every day in 1997 - - | IDL> jdcnv, 1997,1,1,0, jd ;Julian date on Jan 1, 1997 - | IDL> sunpos, jd+ dindgen(365), ra, dec ;RA and Dec for each day - - MODIFICATION HISTORY: - - * Written by Michael R. Greason, STX, 28 October 1988. - * Accept vector arguments, W. Landsman - April,1989 - * Eliminated negative right ascensions - MRG, Hughes STX, 6 May 1992. - * Rewritten using the 1993 Almanac. Keywords added. MRG, HSTX, 10 February 1994. - * Major rewrite, improved accuracy, always return values in degrees - W. Landsman May, 1996 - * Added /RADIAN keyword; W. Landsman; August, 1997 - * Converted to IDL V5.0; W. Landsman; September 1997 - * Converted to python; J. R. Lu; August 2016 - """ - # form time in Julian centuries from 1900.0 - t_obj = Time(mjd, format='mjd') - t = (t_obj.jd - 2415020.0) / 36525.0 - - # form sun's mean longitude - l = (279.696678 + ((36000.768925 * t) % 360.0)) * 3600.0 - - # allow for ellipticity of the orbit (equation of centre) - # using the Earth's mean anomaly ME - me = 358.475844 + ((35999.049750 * t) % 360.0) - ellcor = (6910.1 - 17.2 * t) * np.sin(np.radians(me)) + 72.3 * np.sin( - np.radians(2.0 * me)) - l = l + ellcor - - # allow for the Venus perturbations using the mean anomaly of Venus MV - mv = 212.603219 + ((58517.803875 * t) % 360.0) - vencorr = 4.8 * np.cos(np.radians(299.1017 + mv - me)) + \ - 5.5 * np.cos(np.radians(148.3133 + 2.0 * mv - 2.0 * me)) + \ - 2.5 * np.cos(np.radians(315.9433 + 2.0 * mv - 3.0 * me)) + \ - 1.6 * np.cos(np.radians(345.2533 + 3.0 * mv - 4.0 * me)) + \ - 1.0 * np.cos(np.radians(318.1500 + 3.0 * mv - 5.0 * me)) - l += vencorr - - # Allow for the Mars perturbations using the mean anomaly of Mars MM - mm = 319.529425 + ((19139.858500 * t) % 360.0) - marscorr = 2.0 * np.cos(np.radians(343.8883 - 2.0 * mm + 2.0 * me)) + \ - 1.8 * np.cos(np.radians(200.4017 - 2.0 * mm + me)) - l += marscorr - - # Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ - mj = 225.328328 + ((3034.6920239 * t) % 360.0) - jupcorr = 7.2 * np.cos(np.radians(179.5317 - mj + me)) + \ - 2.6 * np.cos(np.radians(263.2167 - mj)) + \ - 2.7 * np.cos(np.radians(87.1450 - 2.0 * mj + 2.0 * me)) + \ - 1.6 * np.cos(np.radians(109.4933 - 2.0 * mj + me)) - l += jupcorr - - # Allow for the Moons perturbations using the mean elongation of - # the Moon from the Sun D - d = 350.7376814 + ((445267.11422 * t) % 360.0) - mooncorr = 6.5 * np.sin(np.radians(d)) - l += mooncorr - - # Allow for long period terms - longterm = + 6.4 * np.sin(np.radians(231.19 + 20.20 * t)) - l += longterm - l = (l + 2592000.0) % 1296000.0 - longmed = l / 3600.0 - - # Allow for Aberration - l -= 20.5 - - # Allow for Nutation using the longitude of the Moons mean node OMEGA - omega = 259.183275 - ((1934.142008 * t) % 360.0) - l -= 17.2 * np.sin(np.radians(omega)) - - # Form the True Obliquity - oblt = 23.452294 - 0.0130125 * t + ( - 9.2 * np.cos(np.radians(omega))) / 3600.0 - - # Form Right Ascension and Declination - l = l / 3600.0 - l_rad = np.radians(l) - oblt_rad = np.radians(oblt) - ra = np.arctan2(np.sin(l_rad) * np.cos(oblt_rad), np.cos(l_rad)) - - if (len(ra) > 1): - neg = np.where(ra < 0.0)[0] - ra[neg] = ra[neg] + 2.0 * math.pi - - dec = np.arcsin(np.sin(l_rad) * np.sin(oblt_rad)) - - if radians: - oblt = oblt_rad - longmed = np.radians(longmed) - else: - ra = np.degrees(ra) - dec = np.degrees(dec) - - return ra, dec, longmed, oblt diff --git a/flystar/plots.py b/flystar/plots.py index 707d6e6..d077571 100755 --- a/flystar/plots.py +++ b/flystar/plots.py @@ -2212,7 +2212,12 @@ def plot_chi2_dist(tab, Ndetect, motion_model_dict={}, xlim=40, n_bins=50, boot_ # Fitting position and velocity... so subtract 2 to get Ndof n_params = np.nanmean(tab['n_params'][idx]) Ndof = Ndetect - n_params - print(f"Ndof={Ndof}, Ndetect={Ndetect}, Nparams={n_params}") + if len(np.unique(tab['n_params'][idx]))>1: + print("** Warning: using average Ndof for multiple motion models. **") + print("** Consider using plot_chi2_reduced_dist. **") + print(f"Ndof={Ndof:.2f}, Ndetect={Ndetect}, Nparams={n_params:.2f}") + else: + print(f"Ndof={Ndof}, Ndetect={Ndetect}, Nparams={n_params}") chi2_xaxis = np.linspace(0, xlim, xlim*3) chi2_bins = np.linspace(0, xlim, n_bins) @@ -2221,8 +2226,8 @@ def plot_chi2_dist(tab, Ndetect, motion_model_dict={}, xlim=40, n_bins=50, boot_ plt.hist(x[idx], bins=chi2_bins, histtype='step', label='X', density=True) plt.hist(y[idx], bins=chi2_bins, histtype='step', label='Y', density=True) plt.plot(chi2_xaxis, chi2.pdf(chi2_xaxis, Ndof), 'r-', alpha=0.6, - label='$\chi^2$ ' + str(Ndof) + ' dof') - plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) + label='$\chi^2$ ' + str(round(Ndof,2)) + ' dof') + plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(round(Ndof,2))) plt.xlim(0, xlim) plt.legend() @@ -2244,6 +2249,85 @@ def plot_chi2_dist(tab, Ndetect, motion_model_dict={}, xlim=40, n_bins=50, boot_ return +def plot_chi2_reduced_dist(tab, Ndetect, motion_model_dict={}, xlim=8, n_bins=50, boot_err=False): + """ + tab = flystar table + Ndetect = Number of epochs star detected in + """ + chi2_x_list = [] + chi2_y_list = [] + fnd_list = [] # Number of non-NaN error measurements + + motion_model_dict = motion_model.validate_motion_model_dict(motion_model_dict, tab, None) + i_all_detected = np.where(~np.any(np.isnan(tab['t']),axis=1))[0][0] + xt_mod_all, yt_mod_all, xt_mod_err, yt_mod_err = tab.get_star_positions_at_time(tab['t'][i_all_detected], motion_model_dict, allow_alt_models=True) + + for ii in range(len(tab)): + # Ignore the NaNs + fnd = np.argwhere(~np.isnan(tab['xe'][ii,:])) + fnd_list.append(len(fnd)) + + x = tab['x'][ii, fnd] + y = tab['y'][ii, fnd] + if boot_err: + xerr = np.hypot(tab['xe_boot'][ii, fnd], tab['xe'][ii, fnd]) + yerr = np.hypot(tab['ye_boot'][ii, fnd], tab['ye'][ii, fnd]) + else: + xerr = tab['xe'][ii, fnd] + yerr = tab['ye'][ii, fnd] + + fitLineX = xt_mod_all[ii, fnd] + fitLineY = yt_mod_all[ii,fnd] + + diffX = x - fitLineX + diffY = y - fitLineY + sigX = diffX / xerr + sigY = diffY / yerr + + chi2_x = np.sum(sigX**2) + chi2_y = np.sum(sigY**2) + chi2_x_list.append(chi2_x) + chi2_y_list.append(chi2_y) + + x = np.array(chi2_x_list) + y = np.array(chi2_y_list) + fnd = np.array(fnd_list) + + idx = np.where(fnd == Ndetect)[0] + n_params = tab['n_params'] + Ndof = Ndetect - n_params + print("Reduced chi2 for Ndetect="+str(Ndetect)) + chi2_bins = np.linspace(0, xlim, n_bins) + + plt.figure(figsize=(6,4)) + plt.clf() + plt.hist(x[idx]/Ndof[idx], bins=chi2_bins, histtype='step', label='X', density=True, color='tab:blue') + plt.hist(y[idx]/Ndof[idx], bins=chi2_bins, histtype='step', label='Y', density=True, color='tab:orange') + plt.axvline(np.median(x[idx]/Ndof[idx]), color='tab:blue', linestyle='--', label='X median') + plt.axvline(np.median(y[idx]/Ndof[idx]), color='tab:orange', linestyle='--', label='Y median') + plt.title('Reduced chi2, $N_{epoch} = $' + str(Ndetect)) + plt.xlim(0, xlim) + plt.legend() + + chi2red_x = x / Ndof + chi2red_y = y / Ndof + chi2red_t = (x + y) / (2.0 * Ndof + 1*(tab['motion_model_used']=='Parallax')) + + print('Mean reduced chi^2: (Ndetect = {0:d} of {1:d})'.format(len(idx), len(tab))) + fmt = ' {0:s} = {1:.1f} for N_detect and {2:.1f} for all' + med_chi2red_x_f = np.median(chi2red_x[idx]) + med_chi2red_x_a = np.median(chi2red_x) + med_chi2red_y_f = np.median(chi2red_y[idx]) + med_chi2red_y_a = np.median(chi2red_y) + med_chi2red_t_f = np.median(chi2red_t[idx]) + med_chi2red_t_a = np.median(chi2red_t) + print(fmt.format(' X', med_chi2red_x_f, med_chi2red_x_a)) + print(fmt.format(' Y', med_chi2red_y_f, med_chi2red_y_a)) + print(fmt.format('Tot', med_chi2red_t_f, med_chi2red_t_a)) + + return + + def plot_chi2_dist_per_filter(tab, Ndetect, motion_model_dict={}, xlim=40, n_bins=50, filter=None, boot_err=False): """ tab = flystar table @@ -2756,7 +2840,7 @@ def rs(x): chi2_y = np.sum(sigY**2) chi2_m = np.sum(sigM**2) - dof = len(x) - 2 + dof = int(tab['n_fit'][ii]-tab['n_params'][ii]) dofM = len(m) - 1 chi2_red_x = chi2_x / dof @@ -3112,6 +3196,11 @@ def rs(x): x = tab['x0'] y = tab['y0'] r = np.hypot(x, y) + motion_model_dict = motion_model.validate_motion_model_dict(motion_model_dict, tab, None) + i_all_detected = np.where(~np.any(np.isnan(tab['t']),axis=1))[0][0] + cont_times = np.arange(np.min(tab['t'][i_all_detected]), np.max(tab['t'][i_all_detected]), 0.01) + xt_mod_all, yt_mod_all, xt_mod_err, yt_mod_err = tab.get_star_positions_at_time(tab['t'][i_all_detected], motion_model_dict, allow_alt_models=True) + xt_cont_all, yt_cont_all, xt_cont_err, yt_cont_err = tab.get_star_positions_at_time(cont_times, motion_model_dict, allow_alt_models=True) for i in range(Nstars): for ea, epoch_array in enumerate(epoch_array_list): @@ -3242,9 +3331,9 @@ def rs(x): ind = int((row-1)*Ncols + col) paxes = plt.subplot(Nrows, Ncols, ind) - plt.plot(time, fitLineX, 'b-') - plt.plot(time, fitLineX + fitSigX, 'b--') - plt.plot(time, fitLineX - fitSigX, 'b--') + plt.plot(cont_times, xt_cont_all[ii], 'b-') + plt.plot(cont_times, xt_cont_all[ii] + xt_cont_err[ii], 'b--') + plt.plot(cont_times, xt_cont_all[ii] - xt_cont_err[ii], 'b--') print(np.shape(xerr.reshape(len(xerr),))) if not color_time: plt.errorbar(rs(time), rs(x), yerr=rs(xerr), marker='.', color=color, ls='none') @@ -3275,9 +3364,9 @@ def rs(x): ind = int((row-1)*Ncols + col) paxes = plt.subplot(Nrows, Ncols, ind) - plt.plot(time, fitLineY, 'b-') - plt.plot(time, fitLineY + fitSigY, 'b--') - plt.plot(time, fitLineY - fitSigY, 'b--') + plt.plot(cont_times, yt_cont_all[ii], 'b-') + plt.plot(cont_times, yt_cont_all[ii] + yt_cont_err[ii], 'b--') + plt.plot(cont_times, yt_cont_all[ii] - yt_cont_err[ii], 'b--') if not color_time: plt.errorbar(rs(time), rs(y), yerr=rs(yerr), marker='.', color=color, ls='none') else: @@ -3338,8 +3427,8 @@ def rs(x): paxes = plt.subplot(Nrows, Ncols, ind) plt.plot(time, np.zeros(len(time)), 'b-') - plt.plot(time, fitSigX*1e3, 'b--') - plt.plot(time, -fitSigX*1e3, 'b--') + plt.plot(cont_times, xt_cont_err[ii]*1e3, 'b--') + plt.plot(cont_times, -xt_cont_err[ii]*1e3, 'b--') if not color_time: plt.errorbar(rs(time), rs(x - fitLineX)*1e3, yerr=rs(xerr)*1e3, marker='.', color=color, ls='none') else: @@ -3366,8 +3455,8 @@ def rs(x): paxes = plt.subplot(Nrows, Ncols, ind) plt.plot(time, np.zeros(len(time)), 'b-') - plt.plot(time, fitSigY*1e3, 'b--') - plt.plot(time, -fitSigY*1e3, 'b--') + plt.plot(cont_times, yt_cont_err[ii]*1e3, 'b--') + plt.plot(cont_times, -yt_cont_err[ii]*1e3, 'b--') if not color_time: plt.errorbar(rs(time), rs(y - fitLineY)*1e3, yerr=rs(yerr)*1e3, marker='.', color=color, ls='none') else: diff --git a/flystar/startables.py b/flystar/startables.py index a0bf3e3..d75fca9 100644 --- a/flystar/startables.py +++ b/flystar/startables.py @@ -129,7 +129,7 @@ def __init__(self, *args, ref_list=0, **kwargs): # We have to have special handling of meta-data (i.e. info that has # dimensions of n_lists). - meta_tab = ('LIST_TIMES', 'LIST_NAMES') + meta_tab = ('list_times', 'list_names') meta_type = ((float, int), str) for mm in range(len(meta_tab)): meta_test = meta_tab[mm] @@ -156,6 +156,9 @@ def __init__(self, *args, ref_list=0, **kwargs): if meta_arg in kwargs: self.meta[meta_arg] = kwargs[meta_arg] del kwargs[meta_arg] + elif meta_arg.upper() in kwargs: + self.meta[meta_arg] = kwargs[meta_arg.upper()] + del kwargs[meta_arg] for arg in kwargs: if arg in ['name', 'x', 'y', 'm']: @@ -568,8 +571,8 @@ def fit_velocities(self, weighting='var', use_scipy=True, absolute_sigma=True, b if weighting not in ['var', 'std']: raise ValueError(f"fit_velocities: Weighting must either be 'var' or 'std', not {weighting}!") - if ('t' not in self.colnames) and ('LIST_TIMES' not in self.meta): - raise KeyError("fit_velocities: Failed to access time values. No 't' column in table, no 'LIST_TIMES' in meta.") + if ('t' not in self.colnames) and ('list_times' not in self.meta): + raise KeyError("fit_velocities: Failed to access time values. No 't' column in table, no 'list_times' in meta.") # Check if we have the required columns if not all([_ in self.colnames for _ in ['x', 'y']]): @@ -624,7 +627,7 @@ def fit_velocities(self, weighting='var', use_scipy=True, absolute_sigma=True, b if 't' in self.colnames: self['t0'] = self['t'] else: - self['t0'] = self.meta['LIST_TIMES'][0] + self['t0'] = self.meta['list_times'][0] if 'xe' in self.colnames: self['x0_err'] = self['xe'] self['y0_err'] = self['ye'] @@ -639,7 +642,7 @@ def fit_velocities(self, weighting='var', use_scipy=True, absolute_sigma=True, b if 't' in self.colnames: self['t0'] = self['t'][:, 0] else: - self['t0'] = self.meta['LIST_TIMES'][0] + self['t0'] = self.meta['list_times'][0] if 'xe' in self.colnames: self['x0_err'] = self['xe'][:,0] self['y0_err'] = self['ye'][:,0] @@ -756,7 +759,7 @@ def fit_velocity_for_star(self, ss, motion_model_dict, weighting='var', use_scip if 't' in self.colnames: t = np.ma.masked_invalid(self['t'][ss, :].data) else: - t = np.ma.masked_invalid(self.meta['LIST_TIMES']) + t = np.ma.masked_invalid(self.meta['list_times']) if mask_val: t = np.ma.masked_values(t, mask_val) @@ -1086,3 +1089,38 @@ def shift_reference_frame(self, delta_vx=0.0, delta_vy=0.0, delta_pi=0.0, self['x'] += delta_pi*pvec[0] self['y'] += delta_pi*pvec[1] return + +def shift_reference_frame(table, delta_vx=0.0, delta_vy=0.0, delta_pi=0.0, + motion_model_dict={}): + """ + After completing an alignment, shift from your relative reference frame to + the absolute frame using either Gaia or a Galactic model. This modified the + motion model fit parameters as well as the time series astrometry, assuming + zero error on the shift values. + + Parameters + ---------- + delta_vx : float, optional + velocity shift in x-direction (as/yr) + delta_vy : float, optional + velocity shift in y-direction (as/yr) + delta_pi : float, optional + parallax shift (as) + """ + motion_model_dict = motion_model.validate_motion_model_dict(motion_model_dict, table, None) + if delta_vx==0.0 and delta_vy==0.0 and delta_pi==0.0: + print("No shifts input, reference frame unchanged.") + print("Specify delta_vx, delta_vy, and/or delta_pi to perform a reference frame shift.") + return + table['vx'] += delta_vx + table['x'] += delta_vx*(table['t']-table['t0'][:, np.newaxis]) + table['vy'] += delta_vy + table['y'] += delta_vy*(table['t']-table['t0'][:, np.newaxis]) + if delta_pi!=0.0: + t_all = table['t'][np.where(~np.any(np.isnan(table['t']), axis=1))[0][0]] + t_mjd = Time(t_all, format='decimalyear', scale='utc').mjd + pvec = motion_model_dict['Parallax'].get_parallax_vector(t_mjd) + table['pi'] += delta_pi + table['x'] += delta_pi*pvec[0] + table['y'] += delta_pi*pvec[1] + return table diff --git a/flystar/tests/test_align.py b/flystar/tests/test_align.py index 9b65eb6..2d6b0dc 100644 --- a/flystar/tests/test_align.py +++ b/flystar/tests/test_align.py @@ -326,8 +326,8 @@ def test_MosaicToRef_p0_vel(): # The velocities should be almost the same (but not as close as before) # as the input velocities since update_ref == True. assert (msc.ref_table['name']==ref_list['name']).all() - np.testing.assert_allclose(msc.ref_table['vx'], ref_list['vx'], rtol=1e-1) - np.testing.assert_allclose(msc.ref_table['vy'], ref_list['vy'], rtol=1e-1) + assert np.max(np.abs(msc.ref_table['vx']-ref_list['vx']))<3e-4 + assert np.max(np.abs(msc.ref_table['vy']-ref_list['vy']))<3e-4 # Also double check that they aren't exactly the same for the reference stars. #assert np.any(np.not_equal(msc.ref_table['vx'], ref_list['vx'])) @@ -392,8 +392,8 @@ def test_MosaicToRef_vel(): # The velocities should be almost the same (but not as close as before) # as the input velocities since update_ref == True. assert (msc.ref_table['name']==ref_list['name']).all() - np.testing.assert_allclose(msc.ref_table['vx'], ref_list['vx'], rtol=1e-1) - np.testing.assert_allclose(msc.ref_table['vy'], ref_list['vy'], rtol=1e-1) + np.testing.assert_allclose(msc.ref_table['vx'], ref_list['vx'], rtol=1e-1, atol=3e-4) + np.testing.assert_allclose(msc.ref_table['vy'], ref_list['vy'], rtol=1e-1, atol=3e-4) # Also double check that they aren't exactly the same for the reference stars. #assert np.any(np.not_equal(msc.ref_table['vx'], ref_list['vx'])) @@ -476,11 +476,11 @@ def test_MosaicToRef_acc(): if ~np.isnan(msc.ref_table['ax'][ix_fit]): i_orig.append(i) i_fit.append(ix_fit) - np.testing.assert_allclose(msc.ref_table['ax'][i_fit], ref_list['ax'][i_orig], rtol=1e-1) - np.testing.assert_allclose(msc.ref_table['ay'][i_fit], ref_list['ay'][i_orig], rtol=1e-1) + np.testing.assert_allclose(msc.ref_table['ax'][i_fit], ref_list['ax'][i_orig], rtol=1e-1, atol=3e-4) + np.testing.assert_allclose(msc.ref_table['ay'][i_fit], ref_list['ay'][i_orig], rtol=1e-1, atol=3e-4) # Also double check that they aren't exactly the same for the reference stars. - assert np.any(np.not_equal(msc.ref_table['ax'], ref_list['ax'])) + assert np.any(np.not_equal(msc.ref_table['ax'][:200], ref_list['ax'][:200])) return msc