Updated astrometry routines w/ output dictionaries.

Author: iancze

psoap/orbit.py

@@ -7,6 +7,10 @@
       v_{r,2} = K_2 (\cos (\omega_2 + f) + e \cos \omega_2)
 
 where a positive :math:`v_r` indicates a redshifted source and :math:`\omega_2 = \omega_1 + \pi`. In an RV-only paper, typically authors report :math:`\omega = \omega_1`.
+
+In general, quantities like ``vA``, ``vB``, and ``vC`` refer to radial velocities relative to the heliocentric frame of the earth, i.e., quantities to directly compare to actual radial velocity measurements.
+
+(Private) quantities like ``v1``, ``v2``, ``v1_in``, and ``v1_out`` have more subtle meanings---these are generally partial velocities used in the process of constructing ``vA``, etc...
 '''
 
 import numpy as np
@@ -40,9 +44,15 @@ def __init__(self, K, e, omega, P, T0, gamma, obs_dates=None, **kwargs):
         self.obs_dates = obs_dates
 
 
-    def _theta(self, t):
-        '''Calculate the true anomoly for the A-B orbit.
-        Input is in days.'''
+    def _f(self, t):
+        '''Calculate the true anomoly :math:`f` for the A-B orbit.
+
+        Args:
+            t (float): time to evaluate true anomaly at [days].
+
+        Returns:
+            float: true anomaly
+            '''
 
         # t is input in seconds
 
@@ -63,7 +73,10 @@ def _theta(self, t):
 
     def _v1_f(self, f):
         '''Calculate the component of A's velocity based on only the inner orbit.
-        f is the true anomoly of this inner orbit.'''
+
+        Args: f (float): true anomaly of inner orbit.
+
+        Returns: velocity of A (no systemic velocity included).'''
 
         return self.K * (np.cos(self.omega * np.pi/180 + f) + self.e * np.cos(self.omega * np.pi/180))
 
@@ -72,7 +85,7 @@ def _vA_t(self, t):
 
         '''
         # Get the true anomoly "f" from time
-        f = self._theta(t)
+        f = self._f(t)
 
         # Feed this into the orbit equation and add the systemic velocity
         return self._v1_f(f) + self.gamma
@@ -127,7 +140,7 @@ def _v2_f(self, f):
         return self.K/self.q * (np.cos(self.omega_2 * np.pi/180 + f) + self.e * np.cos(self.omega_2 * np.pi/180))
 
     def _vB_t(self, t):
-        f = self._theta(t)
+        f = self._f(t)
 
         # Feed this into the orbit equation and add the systemic velocity
         return self._v2_f(f) + self.gamma
@@ -193,7 +206,7 @@ def __init__(self, K_in, e_in, omega_in, P_in, T0_in, K_out, e_out, omega_out, P
         # just store them to the object.
         self.obs_dates = obs_dates
 
-    def _theta_in(self, t):
+    def _f_in(self, t):
         '''Calculate the true anomoly for the A-B orbit.'''
 
         # t is input in seconds
@@ -213,7 +226,7 @@ def _theta_in(self, t):
         else:
             return th + 2 * np.pi
 
-    def _theta_out(self, t):
+    def _f_out(self, t):
         '''Calculate the true anomoly for the (A-B) - C orbit.'''
 
         # t is input in seconds
@@ -239,24 +252,50 @@ def _v1_f(self, f):
 
         return self.K_in * (np.cos(self.omega_in * np.pi/180 + f) + self.e_in * np.cos(self.omega_in * np.pi/180))
 
-
     def _v3_f(self, f):
         '''Calculate the velocity of (A-B) based only on the outer orbit.
         f is the true anomoly of the outer orbit'''
         return  self.K_out * (np.cos(self.omega_out * np.pi/180 + f) + self.e_out * np.cos(self.omega_out * np.pi/180))
 
-
     def _vA_t(self, t):
 
         # Get the true anomoly "f" from time
-        f_in = self._theta_in(t)
-        f_out = self._theta_out(t)
+        f_in = self._f_in(t)
+        f_out = self._f_out(t)
 
         v1 = self._v1_f(f_in)
         v3 = self._v3_f(f_out)
 
         return v1 + v3 + self.gamma
 
+    def get_component_velocities(self, dates=None):
+        '''
+        Routine to grab v1_inner and v1_outer.
+
+        Args:
+            dates (optional): if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized.
+
+        Returns:
+            np.array: A ``(1, nvels)`` shape array of the primary velocities for the outer motion
+        '''
+
+        if dates is None and self.obs_dates is None:
+            raise RuntimeError("Must provide input dates or specify observation dates upon creation of orbit object.")
+
+        if dates is None and self.obs_dates is not None:
+            dates = self.obs_dates
+
+        dates = np.atleast_1d(dates)
+
+        # calculate f_in and f_out for all dates
+        f_ins = np.array([self._f_in(date) for date in dates])
+        f_outs = np.array([self._f_out(date) for date in dates])
+
+        v1s = np.array([self._v1_f(f) for f in f_ins])
+        v3s = np.array([self._v3_f(f) for f in f_outs])
+
+        return np.vstack((v1s, v3s))
+
     def get_velocities(self, dates=None):
         '''
         Calculate the velocities of the primary (``vA``) for all dates provided.
@@ -280,6 +319,7 @@ def get_velocities(self, dates=None):
 
         return np.atleast_2d(vAs)
 
+
 class ST2(ST1):
     '''
     Orbital model for a double-line (inner orbit) hierarchical spectroscopic triple.
@@ -314,14 +354,43 @@ def _v2_f(self, f):
     def _vB_t(self, t):
 
         # Get the true anolomy "f" from time
-        f_in = self._theta_in(t)
-        f_out = self._theta_out(t)
+        f_in = self._f_in(t)
+        f_out = self._f_out(t)
 
         v2 = self._v2_f(f_in)
         v3 = self._v3_f(f_out)
 
         return v2 + v3 + self.gamma
 
+    def get_component_velocities(self, dates=None):
+        '''
+        Routine to grab v1_inner and v1_outer.
+
+        Args:
+            dates (optional): if provided, calculate velocities at this new vector of dates, rather than the one provided when the object was initialized.
+
+        Returns:
+            np.array: A ``(1, nvels)`` shape array of the primary velocities for the outer motion
+        '''
+
+        if dates is None and self.obs_dates is None:
+            raise RuntimeError("Must provide input dates or specify observation dates upon creation of orbit object.")
+
+        if dates is None and self.obs_dates is not None:
+            dates = self.obs_dates
+
+        dates = np.atleast_1d(dates)
+
+        # calculate f_in and f_out for all dates
+        f_ins = np.array([self._f_in(date) for date in dates])
+        f_outs = np.array([self._f_out(date) for date in dates])
+
+        v1s = np.array([self._v1_f(f) for f in f_ins])
+        v2s = np.array([self._v2_f(f) for f in f_ins])
+        v3s = np.array([self._v3_f(f) for f in f_outs])
+
+        return np.vstack((v1s, v2s, v3s))
+
     def get_velocities(self, dates=None):
         '''
         Calculate the velocities of the primary (``vA``) and secondary (``vB``) for all dates provided.
@@ -384,7 +453,7 @@ def _v3_f_C(self, f):
     def _vC_t(self, t):
 
         # Get the true anolomy "f" from time
-        f_out = self._theta_out(t)
+        f_out = self._f_out(t)
 
         v3 = self._v3_f_C(f_out)
 

psoap/orbit_astrometry.py

@@ -53,14 +53,15 @@ class Binary:
         omega (float): argument of periastron of the *primary*, i.e. :math:`\omega_1` [degrees]
         Omega (float): position angle of the ascending node (going into sky) [deg] east of north
         T0 (float): epoch of periastron passage [JD]
-        M_tot (float): sum of the masses [:math:`M_\odot`]
+        M_tot (float): sum of the masses :math:`M_\mathrm{A} + M_\mathrm{B}` [:math:`M_\odot`]
         M_2 (float): mass of B [:math:`M_\odot`]
         gamma (float): systemic velocity (km/s)
         obs_dates (1D np.array): dates of observation (JD)
     '''
 
     def __init__(self, a, e, i, omega, Omega, T0, M_tot, M_2, gamma, obs_dates=None, **kwargs):
         assert (e >= 0.0) and (e < 1.0), "Eccentricity must be between [0, 1)"
+        assert (i >= 0.0) and (i <= 180.), "Inclination must be between [0, 180]"
         self.a = a # [AU] semi-major axis
         self.e = e # eccentricity
         self.i = i # [deg] inclination
@@ -349,6 +350,8 @@ class Triple:
     def __init__(self, a_in, e_in, i_in, omega_in, Omega_in, T0_in, a_out, e_out, i_out, omega_out, Omega_out, T0_out, M_1, M_2, M_3, gamma, obs_dates=None, **kwargs):
         assert (e_in >= 0.0) and (e_in < 1.0), "Inner eccentricity must be between [0, 1)"
         assert (e_out >= 0.0) and (e_out < 1.0), "Outer eccentricity must be between [0, 1)"
+        assert (i_in >= 0.0) and (i_in <= 180.), "Inner inclination must be between [0, 180]"
+        assert (i_out >= 0.0) and (i_out <= 180.), "Outur inclination must be between [0, 180]"
         self.a_in = a_in # [AU]
         self.e_in = e_in #
         self.i_in = i_in # [deg]