Skip to content

Commit

Permalink
Moved radecdists to utils.py. Updated modules to handle change.
Browse files Browse the repository at this point in the history
  • Loading branch information
@gotten
gotten committed Jan 8, 2024
1 parent fb5e40f commit 4384216
Show file tree
Hide file tree
Showing 3 changed files with 82 additions and 78 deletions.
74 changes: 3 additions & 71 deletions backtracks/backtracks.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@

from schwimmbad import MPIPool

from backtracks.utils import pol2car, transform_gengamm, transform_normal, transform_uniform, utc2tt, HostStarPriors
from backtracks.utils import pol2car, transform_gengamm, transform_normal, transform_uniform, utc2tt, HostStarPriors, radecdists
from backtracks.plotting import diagnostic, neighborhood, plx_prior, posterior, trackplot, stationtrackplot


Expand Down Expand Up @@ -271,7 +271,7 @@ def dummy_loglike(radec):

ra, dec = radec
dummy_param = ra, dec, 0, 0, 0, self.rao, self.deco, self.pmrao, self.pmdeco, self.paro, self.radvelo
xs, ys = self.radecdists([utc2tt(self.ref_epoch)], dummy_param)
xs, ys = radecdists(self, [utc2tt(self.ref_epoch)], dummy_param)

distance = np.linalg.norm(np.array([self.ras[self.ref_epoch_idx],self.decs[self.ref_epoch_idx]])-np.array([xs[0],ys[0]]))

Expand Down Expand Up @@ -373,74 +373,6 @@ def set_prior_attr(self):

print(f'[BACKTRACK INFO]: Queried distance prior parameters, L={self.L:.2f}, alpha={self.alpha:.2f}, beta={self.beta:.2f}')

def radecdists(self, days, param): # for multiple epochs
"""
Function that calculates the offset between companion and host star at a certain set of Epochs assuming background star tracks.
Args:
days (np.array of float): Array of Julian days (Terrestrial Time) at which to calculate the offsets.
param (np.array of float): Array of host star and background star parameters.
Returns:
tuple of arrays: RA and DEC offsets from host star position at Epochs.
References:
* Bangert, J., Puatua, W., Kaplan, G., Bartlett, J., Harris, W., Fredericks, A., & Monet, A. (2011) User's Guide to NOVAS Version C3.1 (Washington, DC: USNO).
* Barron, E. G., Kaplan, G. H., Bangert, J., Bartlett, J. L., Puatua, W., Harris, W., & Barrett, P. (2011) Bull. AAS, 43, 2011.
* `NOVAS website <https://aa.usno.navy.mil/software/novas_info>`__
Notes:
* The NOVAS function `app_star` does the heavy lifting here and assumes a geocentric observer location. \
It takes into account gravitational lensing by solar system bodies, stellar aberration, 3D motion and parallax.
* For Proxima Centauri you may expect a 33 uas offset throughout a day due to an additional small parallax component \
(up to 1 Earth radius offset of the observer w.r.t. geocenter). `topo_star` (not implemented in backtracks) would have to be used to account for this effect.
* The default ephemeris file in `backtracks` is DE405 (valid from 1599 DEC 09 to 2201 FEB 20). To extend time coverage use DE430/DE440.
* Additional JPL Ephemeris files are found `here <https://ssd.jpl.nasa.gov/planets/eph_export.html>`__ \
and `here <https://ssd.jpl.nasa.gov/ftp/eph/planets/>`__
"""

jd_start, jd_end, number = ephem_open() # can't we do this in the System class?

if len(param) == 4:
ra, dec, pmra, pmdec = param
par=0
host_icrs=self.host_icrs
elif len(param) == 5:
ra, dec, pmra, pmdec, par = param
host_icrs=self.host_icrs

else:
ra, dec, pmra, pmdec, par, ra_host, dec_host, pmra_host, pmdec_host, par_host, rv_host = param

host_gaia= novas.make_cat_entry(star_name="HST", catalog="HIP", star_num=1,
ra=ra_host/15., dec=dec_host, pm_ra=pmra_host, pm_dec=pmdec_host,
parallax=par_host, rad_vel=rv_host)

host_icrs = novas.transform_cat(option=1, incat=host_gaia, date_incat=self.gaia_epoch,
date_newcat=2000., newcat_id="HIP")

star2_gaia = novas.make_cat_entry(star_name="BGR", catalog="HIP", star_num=2,
ra=ra/15., dec=dec, pm_ra=pmra, pm_dec=pmdec,
parallax=par, rad_vel=0)

star2_icrs = novas.transform_cat(option=1, incat=star2_gaia, date_incat=self.gaia_epoch,
date_newcat=2000., newcat_id="HIP")

posx=[]
posy=[]
for i, day in enumerate(days):
raa,deca = novas.app_star(day,host_icrs)
rab,decb = novas.app_star(day,star2_icrs)
c_a=SkyCoord(raa,deca,unit=("hourangle","deg"))
c_b=SkyCoord(rab,decb,unit=("hourangle","deg"))
offset=c_a.spherical_offsets_to(c_b)
position_x=offset[0].mas
position_y=offset[1].mas
posx.append(position_x)
posy.append(position_y)

return np.array(posx),np.array(posy)

def fmodel(self, param):
"""
Function that models the offset of the background star with respect to the host star at the observed epochs.
Expand All @@ -451,7 +383,7 @@ def fmodel(self, param):
Returns:
RA and DEC offsets at given observed epochs.
"""
return self.radecdists(self.epochs_tt, param)
return radecdists(self, self.epochs_tt, param)

def loglike(self, param):
"""
Expand Down
14 changes: 7 additions & 7 deletions backtracks/plotting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
from dynesty import plotting as dyplot
from matplotlib.ticker import FuncFormatter

from backtracks.utils import transform_gengamm, radec2seppa, seppa2radec, transform_errors, utc2tt
from backtracks.utils import transform_gengamm, radec2seppa, seppa2radec, transform_errors, utc2tt, radecdists

plt.style.use('default')
plt.rcParams['figure.figsize'] = [9., 6.]
Expand Down Expand Up @@ -200,7 +200,7 @@ def trackplot(
if plot_stationary:
stat_pars = backtracks.stationary_params

ra_stat, dec_stat = backtracks.radecdists(plot_epochs_tt, stat_pars)
ra_stat, dec_stat = radecdists(backtracks, plot_epochs_tt, stat_pars)
axs['A'].plot(ra_stat, dec_stat, color="lightgray", ls='--',
label="Stationary track, $\chi^2_r={}$".format(round(backtracks.stationary_chi2_red,2)))

Expand All @@ -223,7 +223,7 @@ def trackplot(
dec_samples = np.zeros((n_samples, plot_epochs.size))

for i, idx_item in enumerate(random_idx):
ra_samples[i, ], dec_samples[i, ] = backtracks.radecdists(plot_epochs_tt, post_samples[idx_item, ])
ra_samples[i, ], dec_samples[i, ] = radecdists(backtracks, plot_epochs_tt, post_samples[idx_item, ])
axs['A'].plot(ra_samples[i, ], dec_samples[i, ], color='cornflowerblue', lw=0.3, alpha=0.3)

# Create 1 sigma and 3 sigma percentiles for envelopes
Expand Down Expand Up @@ -254,7 +254,7 @@ def trackplot(

# Plot background track with best-fit parameters

ra_bg, dec_bg = backtracks.radecdists(plot_epochs_tt, backtracks.run_median)
ra_bg, dec_bg = radecdists(backtracks, plot_epochs_tt, backtracks.run_median)

axs['A'].plot(ra_bg, dec_bg, color="black", label="Median track, $\chi^2_r={}$".format(round(backtracks.median_chi2_red,2)))

Expand All @@ -269,7 +269,7 @@ def trackplot(

# Connect data points with best-fit model epochs

ra_bg_best, dec_bg_best = backtracks.radecdists(backtracks.epochs_tt, backtracks.run_median)
ra_bg_best, dec_bg_best = radecdists(backtracks, backtracks.epochs_tt, backtracks.run_median)

for i in range(len(backtracks.ras)):
comp_ra = (backtracks.ras[i], ra_bg_best[i])
Expand Down Expand Up @@ -466,7 +466,7 @@ def stationtrackplot(

stat_pars = backtracks.stationary_params

ra_stat, dec_stat = backtracks.radecdists(plot_epochs_tt, stat_pars)
ra_stat, dec_stat = radecdists(backtracks, plot_epochs_tt, stat_pars)
axs['A'].plot(ra_stat, dec_stat, color="lightgray", ls='--',
label="Stationary track, $\chi^2_r={}$".format(round(backtracks.stationary_chi2_red,2)))

Expand All @@ -480,7 +480,7 @@ def stationtrackplot(

# Connect data points with best-fit model epochs

ra_bg_best, dec_bg_best = backtracks.radecdists(backtracks.epochs_tt, stat_pars)
ra_bg_best, dec_bg_best = radecdists(backtracks, backtracks.epochs_tt, stat_pars)

for i in range(len(backtracks.ras)):
comp_ra = (backtracks.ras[i], ra_bg_best[i])
Expand Down
72 changes: 72 additions & 0 deletions backtracks/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,78 @@
from scipy.stats import norm
from scipy.special import gammainc, gammaincc, gammainccinv, gammaincinv
from astropy.time import Time
import novas.compat as novas
from novas.compat.eph_manager import ephem_open
from astropy.coordinates import SkyCoord

def radecdists(backtracks, days, param): # for multiple epochs
"""
Function that calculates the offset between companion and host star at a certain set of Epochs assuming background star tracks.
Args:
backtracks (class): backtracks.System class which carries the needed methods and attributes.
days (np.array of float): Array of Julian days (Terrestrial Time) at which to calculate the offsets.
param (np.array of float): Array of host star and background star parameters.
Returns:
tuple of arrays: RA and DEC offsets from host star position at Epochs.
References:
* Bangert, J., Puatua, W., Kaplan, G., Bartlett, J., Harris, W., Fredericks, A., & Monet, A. (2011) User's Guide to NOVAS Version C3.1 (Washington, DC: USNO).
* Barron, E. G., Kaplan, G. H., Bangert, J., Bartlett, J. L., Puatua, W., Harris, W., & Barrett, P. (2011) Bull. AAS, 43, 2011.
* `NOVAS website <https://aa.usno.navy.mil/software/novas_info>`__
Notes:
* The NOVAS function `app_star` does the heavy lifting here and assumes a geocentric observer location. \
It takes into account gravitational lensing by solar system bodies, stellar aberration, 3D motion and parallax.
* For Proxima Centauri you may expect a 33 uas offset throughout a day due to an additional small parallax component \
(up to 1 Earth radius offset of the observer w.r.t. geocenter). `topo_star` (not implemented in backtracks) would have to be used to account for this effect.
* The default ephemeris file in `backtracks` is DE405 (valid from 1599 DEC 09 to 2201 FEB 20). To extend time coverage use DE430/DE440.
* Additional JPL Ephemeris files are found `here <https://ssd.jpl.nasa.gov/planets/eph_export.html>`__ \
and `here <https://ssd.jpl.nasa.gov/ftp/eph/planets/>`__
"""

jd_start, jd_end, number = ephem_open() # can't we do this in the System class?

if len(param) == 4:
ra, dec, pmra, pmdec = param
par=0
host_icrs=backtracks.host_icrs
elif len(param) == 5:
ra, dec, pmra, pmdec, par = param
host_icrs=backtracks.host_icrs

else:
ra, dec, pmra, pmdec, par, ra_host, dec_host, pmra_host, pmdec_host, par_host, rv_host = param

host_gaia= novas.make_cat_entry(star_name="HST", catalog="HIP", star_num=1,
ra=ra_host/15., dec=dec_host, pm_ra=pmra_host, pm_dec=pmdec_host,
parallax=par_host, rad_vel=rv_host)

host_icrs = novas.transform_cat(option=1, incat=host_gaia, date_incat=backtracks.gaia_epoch,
date_newcat=2000., newcat_id="HIP")

star2_gaia = novas.make_cat_entry(star_name="BGR", catalog="HIP", star_num=2,
ra=ra/15., dec=dec, pm_ra=pmra, pm_dec=pmdec,
parallax=par, rad_vel=0)

star2_icrs = novas.transform_cat(option=1, incat=star2_gaia, date_incat=backtracks.gaia_epoch,
date_newcat=2000., newcat_id="HIP")

posx=[]
posy=[]
for i, day in enumerate(days):
raa,deca = novas.app_star(day,host_icrs)
rab,decb = novas.app_star(day,star2_icrs)
c_a=SkyCoord(raa,deca,unit=("hourangle","deg"))
c_b=SkyCoord(rab,decb,unit=("hourangle","deg"))
offset=c_a.spherical_offsets_to(c_b)
position_x=offset[0].mas
position_y=offset[1].mas
posx.append(position_x)
posy.append(position_y)

return np.array(posx),np.array(posy)

def pol2car(sep, pa, seperr, paerr, corr=np.nan):
"""
Expand Down

0 comments on commit 4384216

Please sign in to comment.