From c6015010afce954282edbbf3da039a7bd83e626f Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Thu, 2 May 2024 08:54:37 -0700 Subject: [PATCH 01/12] add option to exclude hip iad from likelihood --- orbitize/hipparcos.py | 4 ++++ orbitize/sampler.py | 43 ++++++++++++++++++++++--------------------- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 72bc4084..fda8fde5 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -164,6 +164,8 @@ class HipparcosLogProb(object): should be False, but it's helpful for testing. Check out `orbitize.hipparcos.nielsen_iad_refitting_test()` for an example using this renormalization. + include_hip_iad_in_likelihood (bool): if False, then don't add the Hipparcos + log(likelihood) to the overall log(likelihood computed in sampler.py) Written: Sarah Blunt & Rob de Rosa, 2021 """ @@ -175,9 +177,11 @@ def __init__( num_secondary_bodies, alphadec0_epoch=1991.25, renormalize_errors=False, + include_hip_iad_in_likelihood=True ): self.path_to_iad_file = path_to_iad_file self.renormalize_errors = renormalize_errors + self.include_hip_iad_in_likelihood = include_hip_iad_in_likelihood # infer if the IAD file is an older DVD file or a new file with open(path_to_iad_file, "r") as f: diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 5ee50f21..df6c6508 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -100,30 +100,31 @@ def _logl(self, params): lnlikes_sum += self.custom_lnlike(params) if self.system.hipparcos_IAD is not None: - # compute Ra/Dec predictions at the Hipparcos IAD epochs - raoff_model, deoff_model, _ = self.system.compute_all_orbits( - params, epochs=self.system.hipparcos_IAD.epochs_mjd - ) + if self.system.hipparcos_IAD.include_hip_iad_in_likelihood: + # compute Ra/Dec predictions at the Hipparcos IAD epochs + raoff_model, deoff_model, _ = self.system.compute_all_orbits( + params, epochs=self.system.hipparcos_IAD.epochs_mjd + ) - ( - raoff_model_hip_epoch, - deoff_model_hip_epoch, - _, - ) = self.system.compute_all_orbits( - params, epochs=Time([1991.25], format="decimalyear").mjd - ) + ( + raoff_model_hip_epoch, + deoff_model_hip_epoch, + _, + ) = self.system.compute_all_orbits( + params, epochs=Time([1991.25], format="decimalyear").mjd + ) - # subtract off position of star at reference Hipparcos epoch - raoff_model[:, 0, :] -= raoff_model_hip_epoch[:, 0, :] - deoff_model[:, 0, :] -= deoff_model_hip_epoch[:, 0, :] + # subtract off position of star at reference Hipparcos epoch + raoff_model[:, 0, :] -= raoff_model_hip_epoch[:, 0, :] + deoff_model[:, 0, :] -= deoff_model_hip_epoch[:, 0, :] - # select body 0 raoff/deoff predictions & feed into Hip IAD lnlike fn - lnlikes_sum += self.system.hipparcos_IAD.compute_lnlike( - raoff_model[:, 0, :], - deoff_model[:, 0, :], - params, - self.system.param_idx, - ) + # select body 0 raoff/deoff predictions & feed into Hip IAD lnlike fn + lnlikes_sum += self.system.hipparcos_IAD.compute_lnlike( + raoff_model[:, 0, :], + deoff_model[:, 0, :], + params, + self.system.param_idx, + ) if self.system.gaia is not None: gaiahip_epochs = Time( From 31ea1049d86bcab1120e34ed03d7eb96e2ce9197 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 28 May 2024 15:13:26 -0700 Subject: [PATCH 02/12] fix bug with tracking planet perturbations when passing in an array --- orbitize/system.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/orbitize/system.py b/orbitize/system.py index 1da30d3f..cbf5113c 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -454,8 +454,11 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False): within_orbit = np.where(all_smas <= sma) outside_orbit = np.where(all_smas > sma) all_pl_masses = params_arr[self.secondary_mass_indx] - inside_masses = all_pl_masses[within_orbit] - mtot = np.sum(inside_masses) + m0 + inside_masses = all_pl_masses[within_orbit].reshape((-1, n_orbits)) + if n_orbits == 1: + mtot = np.sum(inside_masses) + m0 + else: + mtot = np.sum(inside_masses, axis=0) + m0 else: m_pl = np.zeros(self.num_secondary_bodies) From a75e46182663aed1234c283b56a72935b2e0b41c Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 29 Jul 2024 11:42:37 -0700 Subject: [PATCH 03/12] add ability to fit astrometric jitter --- orbitize/basis.py | 6 ++++++ orbitize/hipparcos.py | 6 +++++- orbitize/sampler.py | 4 ++++ orbitize/system.py | 5 +++++ 4 files changed, 20 insertions(+), 1 deletion(-) diff --git a/orbitize/basis.py b/orbitize/basis.py index 41a582ad..7a629d33 100644 --- a/orbitize/basis.py +++ b/orbitize/basis.py @@ -491,6 +491,7 @@ def __init__( hipparcos_IAD=None, rv=False, rv_instruments=None, + astrometric_jitter=False ): super(Period, self).__init__( @@ -505,6 +506,7 @@ def __init__( rv, rv_instruments, ) + self.astrometric_jitter = astrometric_jitter def construct_priors(self): """ @@ -548,6 +550,10 @@ def construct_priors(self): # Add mass priors self.set_default_mass_priors(basis_priors, basis_labels) + # Add astrometric jitter priors + basis_labels.append("sigma_ast") + basis_priors.append(priors.UniformPrior(0, 10)) + # Define param label dictionary in current basis & standard basis self.param_idx = dict(zip(basis_labels, np.arange(len(basis_labels)))) self.standard_basis_idx = dict(zip(basis_labels, np.arange(len(basis_labels)))) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index fda8fde5..a3934c84 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -429,8 +429,12 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): ) # compute chi2 (Nielsen+ 2020 Eq 7) + if 'sigma_ast' in param_idx: + eps = np.sqrt(self.eps**2 + samples[param_idx["sigma_ast"]]**2) + else: + eps = self.eps chi2 = np.sum( - [(dist[:, i] / self.eps) ** 2 for i in np.arange(n_samples)], + [(dist[:, i] / eps) ** 2 for i in np.arange(n_samples)], axis=1, ) lnlike = -0.5 * chi2 diff --git a/orbitize/sampler.py b/orbitize/sampler.py index df6c6508..0b9eae1c 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -89,6 +89,10 @@ def _logl(self, params): seppa_indices = self.system.all_seppa # compute lnlike + if "sigma_ast" in self.system.param_idx: + errs = np.sqrt(errs**2 + params[self.system.param_idx["sigma_ast"]]) + + lnlikes = self.lnlike( data, errs, corrs, model, jitter, seppa_indices, chi2_type=self.chi2_type ) diff --git a/orbitize/system.py b/orbitize/system.py index cbf5113c..42b02295 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -42,6 +42,9 @@ class System(object): basis to be used. See basis.py for a list of implemented fitting bases. use_rebound (bool): if True, use an n-body backend solver instead of a Keplerian solver. + astrometric_jitter (bool): if True, include a fitted jitter term in the + model (will be applied to Hipparcos data and arbitrary absolute astrometry + data) Priors are initialized as a list of orbitize.priors.Prior objects and stored in the variable ``System.sys_priors``. You should initialize this class, @@ -67,6 +70,7 @@ def __init__( gaia=None, fitting_basis="Standard", use_rebound=False, + astrometric_jitter=False, ): self.num_secondary_bodies = num_secondary_bodies self.data_table = data_table @@ -81,6 +85,7 @@ def __init__( self.gaia = gaia self.fitting_basis = fitting_basis self.use_rebound = use_rebound + self.astrometric_jitter=astrometric_jitter self.best_epochs = [] self.input_table = self.data_table.copy() From 6e36ca71a5cc8a3added77c7b06d319c2067d8c4 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 29 Jul 2024 12:43:48 -0700 Subject: [PATCH 04/12] astrometric jitter fitting working now --- orbitize/basis.py | 5 +++-- orbitize/lnlike.py | 3 +-- orbitize/sampler.py | 4 ++-- orbitize/system.py | 1 + 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/orbitize/basis.py b/orbitize/basis.py index 7a629d33..0238aaf8 100644 --- a/orbitize/basis.py +++ b/orbitize/basis.py @@ -551,8 +551,9 @@ def construct_priors(self): self.set_default_mass_priors(basis_priors, basis_labels) # Add astrometric jitter priors - basis_labels.append("sigma_ast") - basis_priors.append(priors.UniformPrior(0, 10)) + if self.astrometric_jitter: + basis_labels.append("sigma_ast") + basis_priors.append(priors.UniformPrior(0, 10)) # Define param label dictionary in current basis & standard basis self.param_idx = dict(zip(basis_labels, np.arange(len(basis_labels)))) diff --git a/orbitize/lnlike.py b/orbitize/lnlike.py index 8713186c..d1a9a183 100644 --- a/orbitize/lnlike.py +++ b/orbitize/lnlike.py @@ -44,9 +44,8 @@ def chi2_lnlike( in log scale, and defines pa chi-sq using complex phase representation. log sep chisq = (log sep - log sep_true)^2 / (sep_sigma / sep_true)^2 pa chisq = 2 * (1 - cos(pa-pa_true))/pa_sigma^2 -i - """ + """ if np.ndim(model) == 3: # move M dimension to the primary axis, so that numpy knows to iterate over it model = np.rollaxis(model, 2, 0) # now MxNobsx2 in dimensions diff --git a/orbitize/sampler.py b/orbitize/sampler.py index 0b9eae1c..90be012b 100644 --- a/orbitize/sampler.py +++ b/orbitize/sampler.py @@ -89,9 +89,9 @@ def _logl(self, params): seppa_indices = self.system.all_seppa # compute lnlike + # NOTE: this won't work if we're fitting more than just astrometry. This is a hack for betelgeuse. if "sigma_ast" in self.system.param_idx: - errs = np.sqrt(errs**2 + params[self.system.param_idx["sigma_ast"]]) - + jitter = np.ones_like(model) * params[self.system.param_idx["sigma_ast"]] lnlikes = self.lnlike( data, errs, corrs, model, jitter, seppa_indices, chi2_type=self.chi2_type diff --git a/orbitize/system.py b/orbitize/system.py index 42b02295..33e47e97 100644 --- a/orbitize/system.py +++ b/orbitize/system.py @@ -240,6 +240,7 @@ def __init__( hipparcos_IAD=self.hipparcos_IAD, rv=contains_rv, rv_instruments=self.rv_instruments, + astrometric_jitter=self.astrometric_jitter, **self.extra_basis_kwargs ) From 70c3ef21febb88dfd3df8e2bd08ce5d4b2505412 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 29 Jul 2024 13:16:03 -0700 Subject: [PATCH 05/12] fix jitter normalization for hipparcos chi2 calc --- orbitize/hipparcos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index a3934c84..5623c3f7 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -437,7 +437,7 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): [(dist[:, i] / eps) ** 2 for i in np.arange(n_samples)], axis=1, ) - lnlike = -0.5 * chi2 + lnlike = -0.5 * chi2 - np.log(np.sqrt(2 * np.pi * eps)) return lnlike From 103023d1e794d7aef8d5dbbc4bce591e7bc79fdf Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 29 Jul 2024 13:24:57 -0700 Subject: [PATCH 06/12] need to sum jitters to make normalization work --- orbitize/hipparcos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 5623c3f7..187ed395 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -437,7 +437,7 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): [(dist[:, i] / eps) ** 2 for i in np.arange(n_samples)], axis=1, ) - lnlike = -0.5 * chi2 - np.log(np.sqrt(2 * np.pi * eps)) + lnlike = -0.5 * chi2 - np.sum(np.log(np.sqrt(2 * np.pi * eps))) return lnlike From 89a3a58467f2fdeb331ca796997ffabd4c247bf9 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Mon, 26 Aug 2024 16:48:05 -0700 Subject: [PATCH 07/12] plotting tweaks for corner plot --- orbitize/plot.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/orbitize/plot.py b/orbitize/plot.py index a02cfc1e..0b5a1ef3 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -73,25 +73,25 @@ def plot_corner(results, param_list=None, **corner_kwargs): # Define array of default axis labels (overwritten if user specifies list) default_labels = { - "sma": "$a_{0}$ [au]", - "ecc": "$ecc_{0}$", - "inc": "$inc_{0}$ [$^\\circ$]", - "aop": "$\\omega_{0}$ [$^\\circ$]", - "pan": "$\\Omega_{0}$ [$^\\circ$]", - "tau": "$\\tau_{0}$", + "sma": "$a_{{B}}$ [au]", + "ecc": "$ecc_{{B}}$", + "inc": "$inc_{{B}}$ [$^\\circ$]", + "aop": "$\\omega_{{B}}$ [$^\\circ$]", + "pan": "$\\Omega_{{B}}$ [$^\\circ$]", + "tau": "$\\tau_{{B}}$", "tp": "$T_{{\\mathrm{{P}}}}$", "plx": "$\\pi$ [mas]", "gam": "$\\gamma$ [km/s]", "sig": "$\\sigma$ [km/s]", "mtot": "$M_T$ [M$_{{\\odot}}$]", - "m0": "$M_0$ [M$_{{\\odot}}$]", - "m": "$M_{0}$ [M$_{{\\rm Jup}}$]", - "pm_ra": "$\\mu_{{\\alpha}}$ [mas/yr]", + "m0": "$M_{{A}}$ [M$_{{\\odot}}$]", + "m": "$M_{{B}}$ [M$_{{\\odot}}$]", + "pm_ra": "$\\mu_{{\\alpha^{{*}}}}$ [mas/yr]", "pm_dec": "$\\mu_{{\\delta}}$ [mas/yr]", "alpha0": "$\\alpha^{{*}}_{{0}}$ [mas]", "delta0": "$\\delta_0$ [mas]", - "m": "$M_{0}$ [M$_{{\\rm Jup}}$]", - "per": "$P_{0}$ [yr]", + # "m": "$M_{0}$ [M$_{{\\rm Jup}}$]", + "per": "$P_{{B}}$ [yr]", "K": "$K_{0}$ [km/s]", "x": "$X_{0}$ [AU]", "y": "$Y_{0}$ [AU]", @@ -129,9 +129,9 @@ def plot_corner(results, param_list=None, **corner_kwargs): samples[:, angle_indices] = np.degrees( samples[:, angle_indices] ) # convert angles from rad to deg - samples[:, secondary_mass_indices] *= u.solMass.to( - u.jupiterMass - ) # convert to Jupiter masses for companions + # samples[:, secondary_mass_indices] *= u.solMass.to( + # u.jupiterMass + # ) # convert to Jupiter masses for companions if ( "labels" not in corner_kwargs From b3ab0547d3f5824439280647bf318fbd53ca6560 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Wed, 28 Aug 2024 10:41:26 -0700 Subject: [PATCH 08/12] more plotting tweaks --- orbitize/hipparcos.py | 23 +++++++++++++---------- orbitize/plot.py | 14 +++++++------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 187ed395..6c6a4a15 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -521,11 +521,11 @@ def log_prob(model_pars): myHipLogProb.plx0 + 3 * myHipLogProb.plx0_err, 1000, ) - axes[0].hist(sampler.flatchain[:, 0], bins=50, density=True, color="r") + axes[0].hist(sampler.flatchain[:, 0], bins=50, density=True, color="grey") axes[0].plot( xs, norm(myHipLogProb.plx0, myHipLogProb.plx0_err).pdf(xs), color="k" ) - axes[0].set_xlabel("plx [mas]") + axes[0].set_xlabel("$\pi$ [mas]") # PM RA xs = np.linspace( @@ -533,11 +533,11 @@ def log_prob(model_pars): myHipLogProb.pm_ra0 + 3 * myHipLogProb.pm_ra0_err, 1000, ) - axes[1].hist(sampler.flatchain[:, 1], bins=50, density=True, color="r") + axes[1].hist(sampler.flatchain[:, 1], bins=50, density=True, color="grey") axes[1].plot( xs, norm(myHipLogProb.pm_ra0, myHipLogProb.pm_ra0_err).pdf(xs), color="k" ) - axes[1].set_xlabel("PM RA [mas/yr]") + axes[1].set_xlabel("$\mu_{{\\alpha_0^*}}$ [mas/yr]") # PM Dec xs = np.linspace( @@ -545,27 +545,30 @@ def log_prob(model_pars): myHipLogProb.pm_dec0 + 3 * myHipLogProb.pm_dec0_err, 1000, ) - axes[2].hist(sampler.flatchain[:, 2], bins=50, density=True, color="r") + axes[2].hist(sampler.flatchain[:, 2], bins=50, density=True, color="grey") axes[2].plot( xs, norm(myHipLogProb.pm_dec0, myHipLogProb.pm_dec0_err).pdf(xs), color="k" ) - axes[2].set_xlabel("PM Dec [mas/yr]") + axes[2].set_xlabel("$\mu_{{\\delta_0}}$ [mas/yr]") # RA offset - axes[3].hist(sampler.flatchain[:, 3], bins=50, density=True, color="r") + axes[3].hist(sampler.flatchain[:, 3], bins=50, density=True, color="grey") xs = np.linspace( -3 * myHipLogProb.alpha0_err, 3 * myHipLogProb.alpha0_err, 1000 ) axes[3].plot(xs, norm(0, myHipLogProb.alpha0_err).pdf(xs), color="k") - axes[3].set_xlabel("RA Offset [mas]") + axes[3].set_xlabel("$\\alpha_0^*$ [mas]") # Dec offset xs = np.linspace( -3 * myHipLogProb.delta0_err, 3 * myHipLogProb.delta0_err, 1000 ) - axes[4].hist(sampler.flatchain[:, 4], bins=50, density=True, color="r") + axes[4].hist(sampler.flatchain[:, 4], bins=50, density=True, color="grey") axes[4].plot(xs, norm(0, myHipLogProb.delta0_err).pdf(xs), color="k") - axes[4].set_xlabel("Dec Offset [mas]") + axes[4].set_xlabel("$\delta_0$ [mas]") + + for a in axes: + a.set_ylabel('relative prob.') plt.tight_layout() plt.savefig(saveplot, dpi=250) diff --git a/orbitize/plot.py b/orbitize/plot.py index 0b5a1ef3..9eb85a3b 100644 --- a/orbitize/plot.py +++ b/orbitize/plot.py @@ -74,24 +74,24 @@ def plot_corner(results, param_list=None, **corner_kwargs): default_labels = { "sma": "$a_{{B}}$ [au]", - "ecc": "$ecc_{{B}}$", - "inc": "$inc_{{B}}$ [$^\\circ$]", + "ecc": "$ecc$", + "inc": "$i$ [$^\\circ$]", "aop": "$\\omega_{{B}}$ [$^\\circ$]", - "pan": "$\\Omega_{{B}}$ [$^\\circ$]", - "tau": "$\\tau_{{B}}$", + "pan": "$\\Omega$ [$^\\circ$]", + "tau": "$\\tau$", "tp": "$T_{{\\mathrm{{P}}}}$", "plx": "$\\pi$ [mas]", "gam": "$\\gamma$ [km/s]", "sig": "$\\sigma$ [km/s]", "mtot": "$M_T$ [M$_{{\\odot}}$]", - "m0": "$M_{{A}}$ [M$_{{\\odot}}$]", - "m": "$M_{{B}}$ [M$_{{\\odot}}$]", + "m0": "$M_{{a}}$ [M$_{{\\odot}}$]", + "m": "$M_{{b}}$ [M$_{{\\odot}}$]", "pm_ra": "$\\mu_{{\\alpha^{{*}}}}$ [mas/yr]", "pm_dec": "$\\mu_{{\\delta}}$ [mas/yr]", "alpha0": "$\\alpha^{{*}}_{{0}}$ [mas]", "delta0": "$\\delta_0$ [mas]", # "m": "$M_{0}$ [M$_{{\\rm Jup}}$]", - "per": "$P_{{B}}$ [yr]", + "per": "$P$ [yr]", "K": "$K_{0}$ [km/s]", "x": "$X_{0}$ [AU]", "y": "$Y_{0}$ [AU]", From 4432d9188e28b7adea0d7b8abb3f70f941f768b2 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 3 Sep 2024 14:35:37 -0700 Subject: [PATCH 09/12] fix error message that was printing incorrectly but not affecting code --- orbitize/hipparcos.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index fda8fde5..a3cab7d9 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -286,14 +286,14 @@ def __init__( else: iad = np.transpose(np.loadtxt(path_to_iad_file)) - n_lines = len(iad) - times = iad[1] + 1991.25 self.cos_phi = iad[3] # scan direction self.sin_phi = iad[4] self.R = iad[5] # abscissa residual [mas] self.eps = iad[6] # error on abscissa residual [mas] + n_lines = len(self.eps) + # reject negative errors (scans that were rejected by Hipparcos team) good_scans = np.where(self.eps > 0)[0] From e38fed075f359d4eb1359f22fa6a3c65d3f1b062 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Fri, 6 Sep 2024 16:06:49 -0700 Subject: [PATCH 10/12] add in rob's empirical var correction for type 1 sol --- orbitize/hipparcos.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 6c6a4a15..79828c19 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -262,8 +262,8 @@ def __init__( self.solution_type = solution_details["isol_n"].values[0] - if self.solution_type == 1: - self.var = astrometric_solution["var"].values[0] + if self.solution_type == 1: + self.var = (10 * astrometric_solution["var"].values[0])**2 # N.B. input is different units than var from Vizier catalog!! else: self.var = 0 From aa133fc765a243b5ced3e71826f951780466d795 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 10 Sep 2024 15:20:20 -0700 Subject: [PATCH 11/12] incorporate option for ad-hoc brandt+ 23 correction to iad --- orbitize/hipparcos.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 79828c19..2469cfd4 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -166,6 +166,8 @@ class HipparcosLogProb(object): using this renormalization. include_hip_iad_in_likelihood (bool): if False, then don't add the Hipparcos log(likelihood) to the overall log(likelihood computed in sampler.py) + brandt_correction (bool): if True, add delta_r = 0.140 mas and sigma_jit = 2.25 mas + to the residuals, following Brandt+ 2023. Written: Sarah Blunt & Rob de Rosa, 2021 """ @@ -177,7 +179,8 @@ def __init__( num_secondary_bodies, alphadec0_epoch=1991.25, renormalize_errors=False, - include_hip_iad_in_likelihood=True + include_hip_iad_in_likelihood=True, + brandt_correction=False, ): self.path_to_iad_file = path_to_iad_file self.renormalize_errors = renormalize_errors @@ -262,8 +265,10 @@ def __init__( self.solution_type = solution_details["isol_n"].values[0] - if self.solution_type == 1: - self.var = (10 * astrometric_solution["var"].values[0])**2 # N.B. input is different units than var from Vizier catalog!! + if self.solution_type == 1: + self.var = ( + 10 * astrometric_solution["var"].values[0] + ) ** 2 # N.B. input is different units than var from Vizier catalog!! else: self.var = 0 @@ -292,6 +297,8 @@ def __init__( self.cos_phi = iad[3] # scan direction self.sin_phi = iad[4] self.R = iad[5] # abscissa residual [mas] + if brandt_correction: + self.R += 0.140 # [mas] self.eps = iad[6] # error on abscissa residual [mas] # reject negative errors (scans that were rejected by Hipparcos team) @@ -308,6 +315,12 @@ def __init__( # if the star has a type 1 (stochastic) solution, we need to undo the addition of a jitter term in quadrature self.eps = np.sqrt(self.eps**2 - self.var**2) + if brandt_correction: + self.eps = np.sqrt(self.eps**2 + self.var**2) + import pdb + + pdb.set_trace() + epochs = Time(times, format="decimalyear") self.epochs = epochs.decimalyear self.epochs_mjd = epochs.mjd @@ -429,8 +442,8 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx): ) # compute chi2 (Nielsen+ 2020 Eq 7) - if 'sigma_ast' in param_idx: - eps = np.sqrt(self.eps**2 + samples[param_idx["sigma_ast"]]**2) + if "sigma_ast" in param_idx: + eps = np.sqrt(self.eps**2 + samples[param_idx["sigma_ast"]] ** 2) else: eps = self.eps chi2 = np.sum( @@ -568,7 +581,7 @@ def log_prob(model_pars): axes[4].set_xlabel("$\delta_0$ [mas]") for a in axes: - a.set_ylabel('relative prob.') + a.set_ylabel("relative prob.") plt.tight_layout() plt.savefig(saveplot, dpi=250) From 6e63d8e150d6b4c8cd75fd238a7ca4c2c1d1e202 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 10 Sep 2024 15:27:07 -0700 Subject: [PATCH 12/12] clean up --- orbitize/hipparcos.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py index 2469cfd4..64f61412 100644 --- a/orbitize/hipparcos.py +++ b/orbitize/hipparcos.py @@ -316,10 +316,7 @@ def __init__( self.eps = np.sqrt(self.eps**2 - self.var**2) if brandt_correction: - self.eps = np.sqrt(self.eps**2 + self.var**2) - import pdb - - pdb.set_trace() + self.eps = np.sqrt(self.eps**2 + 2.25**2) epochs = Time(times, format="decimalyear") self.epochs = epochs.decimalyear