From 71569d6f4cf744d89579b269e7a88a0f8153ecaa Mon Sep 17 00:00:00 2001 From: David Trevascus Date: Fri, 19 Jul 2024 18:29:06 +0200 Subject: [PATCH 1/6] Add method to system object to plot radec astrometry with error bars --- orbitize/system.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/orbitize/system.py b/orbitize/system.py index 1da30d3f..4e1aec99 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -2,6 +2,7 @@ from orbitize import nbody, kepler, basis, hipparcos from astropy import table from orbitize.read_input import read_file +import matplotlib.pyplot as plt class System(object): @@ -717,6 +718,33 @@ def convert_data_table_radec2seppa(self, body_num=1): ) self.seppa[body_num] = np.append(self.seppa[body_num], i) + def plot_astrometry(self): + """ + Plot astrometry to ensure data is correct. + + Written: David Trevascus, 2024 + """ + + # create figure + fig, axs = plt.subplots() + + # get radec astrometry + ra = self.data_table[self.all_radec]['quant1'] + ra_err = self.data_table[self.all_radec]['quant1_err'] + dec = self.data_table[self.all_radec]['quant2'] + dec_err = self.data_table[self.all_radec]['quant2_err'] + + # plot radec astrometry + axs.errorbar( + ra, + dec, + xerr=ra_err, + yerr=dec_err, + linestyle='None', + marker='o', + color='k', + capsize=3 + ) def radec2seppa(ra, dec, mod180=False): """ From 33b499d94154335210c5a389aab4e5f79a2408ce Mon Sep 17 00:00:00 2001 From: DTCupcakes Date: Fri, 19 Jul 2024 22:49:36 +0200 Subject: [PATCH 2/6] Add plotting of seppa data (in radec) with errors --- orbitize/system.py | 37 ++++++++++++++++++++++++++++++++----- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index 4e1aec99..20e20d61 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -728,11 +728,38 @@ def plot_astrometry(self): # create figure fig, axs = plt.subplots() - # get radec astrometry - ra = self.data_table[self.all_radec]['quant1'] - ra_err = self.data_table[self.all_radec]['quant1_err'] - dec = self.data_table[self.all_radec]['quant2'] - dec_err = self.data_table[self.all_radec]['quant2_err'] + n_obs = len(self.data_table) # number of entries in the data table + + ra = np.zeros(n_obs) + ra_err = np.zeros(n_obs) + dec = np.zeros(n_obs) + dec_err = np.zeros(n_obs) + radec_corr = np.zeros(n_obs) + + for i in range(n_obs): + if np.isin(i, self.all_radec): # get radec astrometry + ra[i] = self.data_table[i]['quant1'] + ra_err[i] = self.data_table[i]['quant1_err'] + dec[i] = self.data_table[i]['quant2'] + dec_err[i] = self.data_table[i]['quant2_err'] + radec_corr[i] = self.data_table[i]['quant12_corr'] + + elif np.isin(i, self.all_seppa): # get seppa astrometry, convert to radec + sep = self.data_table[i]['quant1'] + sep_err = self.data_table[i]['quant1_err'] + pa = self.data_table[i]['quant2'] + pa_err = self.data_table[i]['quant2_err'] + seppa_corr = self.data_table[i]['quant12_corr'] + + ra[i], dec[i] = seppa2radec(sep, pa) + + if np.isnan(seppa_corr): + ra_err[i] = np.sqrt((ra[i]/sep * sep_err)**2 + (dec[i] * np.radians(pa_err))**2) + dec_err[i] = np.sqrt((dec[i]/sep * sep_err)**2 + (ra[i] * np.radians(pa_err))**2) + else: + ra_err[i], dec_err[i], radec_corr[i] = transform_errors( + sep, pa, sep_err, pa_err, seppa_corr, seppa2radec + ) # plot radec astrometry axs.errorbar( From efdd2e4db8d4abd3db6d6734a7c2d6cdc52febe2 Mon Sep 17 00:00:00 2001 From: DTCupcakes Date: Fri, 19 Jul 2024 23:49:03 +0200 Subject: [PATCH 3/6] Add axis labels to astrometry plot and set axis limits so star is centred --- orbitize/system.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index 20e20d61..fba72479 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -725,17 +725,16 @@ def plot_astrometry(self): Written: David Trevascus, 2024 """ - # create figure - fig, axs = plt.subplots() - n_obs = len(self.data_table) # number of entries in the data table + # set empty arrays for radec astrometry ra = np.zeros(n_obs) ra_err = np.zeros(n_obs) dec = np.zeros(n_obs) dec_err = np.zeros(n_obs) radec_corr = np.zeros(n_obs) + # get radec astrometry for i in range(n_obs): if np.isin(i, self.all_radec): # get radec astrometry ra[i] = self.data_table[i]['quant1'] @@ -761,6 +760,9 @@ def plot_astrometry(self): sep, pa, sep_err, pa_err, seppa_corr, seppa2radec ) + # create figure + fig, axs = plt.subplots() + # plot radec astrometry axs.errorbar( ra, @@ -773,6 +775,17 @@ def plot_astrometry(self): capsize=3 ) + # set axis labels + axs.set_xlabel('RA offset (mas)') + axs.set_ylabel('Dec offset (mas)') + + # set plot limits so star is centred + max_ra = np.amax(abs(ra)) + max_dec = np.amax(abs(dec)) + max_err = np.amax(np.append(ra_err, dec_err)) + axs.set_xlim(-(max_ra + 2 * max_err), max_ra + 2 * max_err) + axs.set_ylim(-(max_dec + 2 * max_err), max_dec + 2 * max_err) + def radec2seppa(ra, dec, mod180=False): """ Convenience function for converting from From 377804365fa9e4f86252ead1cd977521ea6cc317 Mon Sep 17 00:00:00 2001 From: DTCupcakes Date: Sat, 20 Jul 2024 00:14:20 +0200 Subject: [PATCH 4/6] Add astrometry plotting to init method of system class --- orbitize/system.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/orbitize/system.py b/orbitize/system.py index fba72479..b303a2e8 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -303,6 +303,9 @@ def __init__( self.param_idx = self.basis.param_idx + # plot astrometry + self.plot_astrometry() + def save(self, hf): """ Saves the current object to an hdf5 file From 55ab5941b4408a26dd3e28773310ede26f694f21 Mon Sep 17 00:00:00 2001 From: DTCupcakes Date: Thu, 25 Jul 2024 04:17:33 +0200 Subject: [PATCH 5/6] Label objects separately and set option to centre plot on star --- orbitize/system.py | 122 ++++++++++++++++++++++++++------------------- 1 file changed, 70 insertions(+), 52 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index b303a2e8..b1a9c98d 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -721,73 +721,91 @@ def convert_data_table_radec2seppa(self, body_num=1): ) self.seppa[body_num] = np.append(self.seppa[body_num], i) - def plot_astrometry(self): + def plot_astrometry(self, center_on_star=False): """ Plot astrometry to ensure data is correct. + Args: + center_on_star(bool, optional): center plot on star + Written: David Trevascus, 2024 """ - - n_obs = len(self.data_table) # number of entries in the data table - - # set empty arrays for radec astrometry - ra = np.zeros(n_obs) - ra_err = np.zeros(n_obs) - dec = np.zeros(n_obs) - dec_err = np.zeros(n_obs) - radec_corr = np.zeros(n_obs) - # get radec astrometry - for i in range(n_obs): - if np.isin(i, self.all_radec): # get radec astrometry - ra[i] = self.data_table[i]['quant1'] - ra_err[i] = self.data_table[i]['quant1_err'] - dec[i] = self.data_table[i]['quant2'] - dec_err[i] = self.data_table[i]['quant2_err'] - radec_corr[i] = self.data_table[i]['quant12_corr'] - - elif np.isin(i, self.all_seppa): # get seppa astrometry, convert to radec - sep = self.data_table[i]['quant1'] - sep_err = self.data_table[i]['quant1_err'] - pa = self.data_table[i]['quant2'] - pa_err = self.data_table[i]['quant2_err'] - seppa_corr = self.data_table[i]['quant12_corr'] - - ra[i], dec[i] = seppa2radec(sep, pa) - - if np.isnan(seppa_corr): - ra_err[i] = np.sqrt((ra[i]/sep * sep_err)**2 + (dec[i] * np.radians(pa_err))**2) - dec_err[i] = np.sqrt((dec[i]/sep * sep_err)**2 + (ra[i] * np.radians(pa_err))**2) - else: - ra_err[i], dec_err[i], radec_corr[i] = transform_errors( - sep, pa, sep_err, pa_err, seppa_corr, seppa2radec - ) + n_obj = len(np.unique(self.data_table['object'])) # number of planets in the system # create figure fig, axs = plt.subplots() - # plot radec astrometry - axs.errorbar( - ra, - dec, - xerr=ra_err, - yerr=dec_err, - linestyle='None', - marker='o', - color='k', - capsize=3 - ) + # plot each object separately + for n in range(n_obj): + + # separate out data table entries for this object + mask_obj = self.data_table['object'] == n + 1 + data_table_obj = self.data_table[mask_obj] + + # number of entries for this object + n_obs = len(data_table_obj) # number of entries in the data table + + # set empty arrays for radec astrometry + ra = np.zeros(n_obs) + ra_err = np.zeros(n_obs) + dec = np.zeros(n_obs) + dec_err = np.zeros(n_obs) + radec_corr = np.zeros(n_obs) + + # get radec astrometry + for i in range(n_obs): + if np.isin(i, self.all_radec): # get radec astrometry + ra[i] = data_table_obj[i]['quant1'] + ra_err[i] = data_table_obj[i]['quant1_err'] + dec[i] = data_table_obj[i]['quant2'] + dec_err[i] = data_table_obj[i]['quant2_err'] + radec_corr[i] = data_table_obj[i]['quant12_corr'] + + elif np.isin(i, self.all_seppa): # get seppa astrometry, convert to radec + sep = data_table_obj[i]['quant1'] + sep_err = data_table_obj[i]['quant1_err'] + pa = data_table_obj[i]['quant2'] + pa_err = data_table_obj[i]['quant2_err'] + seppa_corr = data_table_obj[i]['quant12_corr'] + + ra[i], dec[i] = seppa2radec(sep, pa) + + if np.isnan(seppa_corr): + ra_err[i] = np.sqrt((ra[i]/sep * sep_err)**2 + (dec[i] * np.radians(pa_err))**2) + dec_err[i] = np.sqrt((dec[i]/sep * sep_err)**2 + (ra[i] * np.radians(pa_err))**2) + else: + ra_err[i], dec_err[i], radec_corr[i] = transform_errors( + sep, pa, sep_err, pa_err, seppa_corr, seppa2radec + ) + + # plot radec astrometry + axs.errorbar( + ra, + dec, + xerr=ra_err, + yerr=dec_err, + linestyle='None', + marker='o', + #color='k', + capsize=3, + label=str(n + 1) + ) # set axis labels axs.set_xlabel('RA offset (mas)') axs.set_ylabel('Dec offset (mas)') + axs.legend() + + if center_on_star == True: + # plot start at centre + axs.plot(0, 0, marker='*', color='k', markersize=12) - # set plot limits so star is centred - max_ra = np.amax(abs(ra)) - max_dec = np.amax(abs(dec)) - max_err = np.amax(np.append(ra_err, dec_err)) - axs.set_xlim(-(max_ra + 2 * max_err), max_ra + 2 * max_err) - axs.set_ylim(-(max_dec + 2 * max_err), max_dec + 2 * max_err) + # set plot limits so star is centred + xlim = max(axs.get_xlim()) + ylim = max(axs.get_ylim()) + axs.set_xlim(-xlim, xlim) + axs.set_ylim(-ylim, ylim) def radec2seppa(ra, dec, mod180=False): """ From f3ba63fdaa7f50687d3e0a85cbc4215515a61aee Mon Sep 17 00:00:00 2001 From: David Trevascus Date: Mon, 12 Aug 2024 12:19:02 +0200 Subject: [PATCH 6/6] Fix issue where sep/pa astrometry was being plotted as ra/dec --- orbitize/system.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index b1a9c98d..4c28466e 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -755,14 +755,15 @@ def plot_astrometry(self, center_on_star=False): # get radec astrometry for i in range(n_obs): - if np.isin(i, self.all_radec): # get radec astrometry + j = np.where(mask_obj)[0][i] # index obj entries in input table + if np.isin(j, self.all_radec): # get radec astrometry ra[i] = data_table_obj[i]['quant1'] ra_err[i] = data_table_obj[i]['quant1_err'] dec[i] = data_table_obj[i]['quant2'] dec_err[i] = data_table_obj[i]['quant2_err'] radec_corr[i] = data_table_obj[i]['quant12_corr'] - elif np.isin(i, self.all_seppa): # get seppa astrometry, convert to radec + elif np.isin(j, self.all_seppa): # get seppa astrometry, convert to radec sep = data_table_obj[i]['quant1'] sep_err = data_table_obj[i]['quant1_err'] pa = data_table_obj[i]['quant2']