From 99e7a688736919ee642ba3f89b691804c485ad7f Mon Sep 17 00:00:00 2001
From: Matthew Pitkin <matthew.pitkin@ligo.org>
Date: Mon, 4 Feb 2019 23:33:53 +0000
Subject: [PATCH] An attempt at allowing the 'TCB' units in PINT:  - refs
 nanograv/PINT#37 and nanograv/PINT#151

---
 pint/models/astrometry.py       | 18 +++++--
 pint/models/binary_ell1.py      |  3 ++
 pint/models/dispersion_model.py | 22 +++++++--
 pint/models/pulsar_binary.py    | 18 +++++--
 pint/models/spindown.py         | 20 ++++++--
 pint/models/timing_model.py     | 41 +++++++++++-----
 pint/models/wave.py             |  6 ++-
 pint/observatory/observatory.py | 87 ++++++++++++++++++++++++++++++++-
 pint/observatory/topo_obs.py    |  3 +-
 pint/toa.py                     | 75 +++++++++++++++++++++++++++-
 10 files changed, 258 insertions(+), 35 deletions(-)

diff --git a/pint/models/astrometry.py b/pint/models/astrometry.py
index c013f5991..08ecef031 100644
--- a/pint/models/astrometry.py
+++ b/pint/models/astrometry.py
@@ -42,7 +42,6 @@ def __init__(self):
     def setup(self):
         super(Astrometry, self).setup()
 
-
     def ssb_to_psb_xyz_ICRS(self, epoch=None):
         """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS.
 
@@ -54,7 +53,10 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None):
     def barycentric_radio_freq(self, toas):
         """Return radio frequencies (MHz) of the toas corrected for Earth motion"""
         tbl = toas.table
-        L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(numpy.float64))
+        if self.UNITS.value == 'TCB':
+            L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tcbld'].astype(numpy.float64))
+        else:
+            L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(numpy.float64))
         v_dot_L_array = numpy.sum(tbl['ssb_obs_vel']*L_hat, axis=1)
         return tbl['freq'] * (1.0 - v_dot_L_array / const.c)
 
@@ -66,7 +68,10 @@ def solar_system_geometric_delay(self, toas, acc_delay=None):
         available as 3-vector toa.xyz, in units of light-seconds.
         """
         tbl = toas.table
-        L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(numpy.float64))
+        if self.UNITS.value == 'TCB':
+            L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tcbld'].astype(numpy.float64))
+        else:
+            L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(numpy.float64))
         re_dot_L = numpy.sum(tbl['ssb_obs_pos']*L_hat, axis=1)
         delay = -re_dot_L.to(ls).value
         if self.PX.value != 0.0 \
@@ -88,7 +93,10 @@ def get_d_delay_quantities(self, toas):
         # TODO: tbl['tdbld'].quantity should have units of u.day
         # NOTE: Do we need to include the delay here?
         tbl = toas.table
-        rd['epoch'] = tbl['tdbld'].quantity * u.day #- delay * u.second
+        if self.UNITS.value == 'TCB':
+            rd['epoch'] = tbl['tcbld'].quantity * u.day #- delay * u.second
+        else:
+            rd['epoch'] = tbl['tdbld'].quantity * u.day #- delay * u.second
 
         # Distance from SSB to observatory, and from SSB to psr
         ssb_obs = tbl['ssb_obs_pos'].quantity
@@ -190,6 +198,7 @@ def setup(self):
                             "POSEPOCH or PEPOCH are required if PM is set.")
                 else:
                     self.POSEPOCH.quantity = self.PEPOCH.quantity
+            self.POSEPOCH.scale = self.UNITS.value.lower()
 
     def print_par(self):
         result = ''
@@ -358,6 +367,7 @@ def setup(self):
                             "POSEPOCH or PEPOCH are required if PM is set.")
                 else:
                     self.POSEPOCH.quantity = self.PEPOCH.quantity
+                    self.POSEPOCH.scale = self.UNITS.value.lower()
 
     def get_psr_coords(self, epoch=None):
         """Returns pulsar sky coordinates as an astropy ecliptic oordinates
diff --git a/pint/models/binary_ell1.py b/pint/models/binary_ell1.py
index ca876230f..24397a434 100644
--- a/pint/models/binary_ell1.py
+++ b/pint/models/binary_ell1.py
@@ -61,10 +61,13 @@ def setup(self):
                 warn("Since ECC is 0.0, using T0 as TASC.")
                 if self.T0.value is not None:
                     self.TASC.value = self.T0.value
+                    self.TASC.scale = self.UNITS.value.lower()
                 else:
                     raise MissingParameter("ELL1", 'T0', "T0 or TASC is required for ELL1 model.")
             else:
                 raise MissingParameter("ELL1", 'TASC', "TASC is required for ELL1 model.")
+        else:
+            self.TASC.scale = self.UNITS.value.lower()
 
 
 
diff --git a/pint/models/dispersion_model.py b/pint/models/dispersion_model.py
index 0205cd89e..b6afba0c7 100644
--- a/pint/models/dispersion_model.py
+++ b/pint/models/dispersion_model.py
@@ -80,6 +80,8 @@ def setup(self):
             if self.DMEPOCH.value is None:
                 raise MissingParameter("Dispersion", "DMEPOCH",
                         "DMEPOCH is required if DM1 or higher are set")
+            else:
+                self.DMEPOCH.scale = self.UNITS.value.lower()
         base_dms = list(self.get_prefix_mapping_component('DM').values())
         base_dms += ['DM',]
 
@@ -105,10 +107,16 @@ def base_dm(self, toas):
         dm = np.zeros(len(tbl))
         dm_terms = self.get_DM_terms()
         if self.DMEPOCH.value is None:
-            DMEPOCH = tbl['tdbld'][0]
+            if self.UNITS.value == 'TCB':
+                DMEPOCH = tbl['tcbld'][0]
+            else:
+                DMEPOCH = tbl['tdbld'][0]
         else:
             DMEPOCH = self.DMEPOCH.value
-        dt = (tbl['tdbld'] - DMEPOCH) * u.day
+        if self.UNITS.value == 'TCB':
+            dt = (tbl['tcbld'] - DMEPOCH) * u.day
+        else:
+            dt = (tbl['tdbld'] - DMEPOCH) * u.day
         dt_value = (dt.to(u.yr)).value
         dm_terms_value = [d.value for d in dm_terms]
         dm = taylor_horner(dt_value, dm_terms_value)
@@ -156,10 +164,16 @@ def d_delay_d_DMs(self, toas, param_name, acc_delay=None): # NOTE we should have
         dm_terms = np.longdouble(np.zeros(len(dms)))
         dm_terms[order] = np.longdouble(1.0)
         if self.DMEPOCH.value is None:
-            DMEPOCH = tbl['tdbld'][0]
+            if self.UNITS.value == 'TCB':
+                DMEPOCH = tbl['tcbld'][0]
+            else:
+                DMEPOCH = tbl['tdbld'][0]
         else:
             DMEPOCH = self.DMEPOCH.value
-        dt = (tbl['tdbld'] - DMEPOCH) * u.day
+        if self.UNITS.value == 'TCB':
+            dt = (tbl['tcbld'] - DMEPOCH) * u.day
+        else:
+            dt = (tbl['tdbld'] - DMEPOCH) * u.day
         dt_value = (dt.to(u.yr)).value
         d_dm_d_dm_param = taylor_horner(dt_value, dm_terms)* (self.DM.units/par.units)
         return DMconst * d_dm_d_dm_param/ bfreq**2.0
diff --git a/pint/models/pulsar_binary.py b/pint/models/pulsar_binary.py
index e606f1c27..b9bbb65ff 100644
--- a/pint/models/pulsar_binary.py
+++ b/pint/models/pulsar_binary.py
@@ -32,7 +32,6 @@ def __init__(self,):
             units=u.day,
             description="Orbital period", long_double=True))
 
-
         self.add_param(p.floatParameter(name="PBDOT",
             units = u.day/u.day,
             description="Orbital period derivitve respect to time", \
@@ -89,6 +88,11 @@ def __init__(self,):
 
     def setup(self):
         super(PulsarBinary, self).setup()
+
+        if self.T0.value is not None:
+            # set the units
+            self.T0.scale = self.UNITS.value.lower()
+
         for bpar in self.params:
             self.register_deriv_funcs(self.d_binary_delay_d_xxxx, bpar)
         # Setup the model isinstance
@@ -137,17 +141,21 @@ def update_binary_object(self, toas, acc_delay=None):
         """
         # Don't need to fill P0 and P1. Translate all the others to the format
         # that is used in bmodel.py
-        # Get barycnetric toa first
+        # Get barycentric toa first
         updates = {}
         tbl = toas.table
         if acc_delay is None:
             # If the accumulate delay is not provided, it will try to get
             # the barycentric correction.
             acc_delay = self.delay(toas, self.__class__.__name__, False)
-        self.barycentric_time = tbl['tdbld'] * u.day - acc_delay
+        if self.UNITS.value == 'TCB':
+            self.barycentric_time = tbl['tcbld'] * u.day - acc_delay
+            updates['psr_pos'] = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tcbld'].astype(np.float64))
+        else:
+            self.barycentric_time = tbl['tdbld'] * u.day - acc_delay
+            updates['psr_pos'] = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(np.float64))
         updates['barycentric_toa'] = self.barycentric_time
         updates['obs_pos'] = tbl['ssb_obs_pos'].quantity
-        updates['psr_pos'] = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(np.float64))
         for par in self.binary_instance.binary_params:
             binary_par_names = [par,]
             if par in self.binary_instance.param_aliases.keys():
@@ -183,7 +191,7 @@ def binarymodel_delay(self, toas, acc_delay=None):
         return self.binary_instance.binary_delay()
 
     def d_binary_delay_d_xxxx(self, toas, param, acc_delay):
-        """Return the bianry model delay derivtives"""
+        """Return the binary model delay derivtives"""
         self.update_binary_object(toas, acc_delay)
         return self.binary_instance.d_binarydelay_d_par(param)
 
diff --git a/pint/models/spindown.py b/pint/models/spindown.py
index c1d4c0e70..6deb36b8f 100644
--- a/pint/models/spindown.py
+++ b/pint/models/spindown.py
@@ -44,6 +44,10 @@ def setup(self):
             if getattr(self, p).value is None:
                 raise MissingParameter("Spindown", p)
 
+        if self.PEPOCH.value is not None:
+            # set the units
+            self.PEPOCH.scale = self.UNITS.value.lower()
+
         # Check continuity
         F_terms = list(self.get_prefix_mapping_component('F').keys())
         F_terms.sort()
@@ -87,12 +91,20 @@ def get_dt(self, toas, delay):
         pulse phase will be handled at a higher level in the code.
         """
         tbl = toas.table
-        if self.PEPOCH.value is None:
-            phsepoch_ld = time_to_longdouble(tbl['tdb'][0] - delay[0])
+        if self.PEPOCH.scale == 'tcb':
+            if self.PEPOCH.value is None:
+                phsepoch_ld = time_to_longdouble(tbl['tcb'][0] - delay[0])
+            else:
+                phsepoch_ld = time_to_longdouble(self.PEPOCH.quantity)
+
+            dt = (tbl['tcbld'] - phsepoch_ld)*u.day - delay
         else:
-            phsepoch_ld = time_to_longdouble(self.PEPOCH.quantity)
+            if self.PEPOCH.value is None:
+                phsepoch_ld = time_to_longdouble(tbl['tdb'][0] - delay[0])
+            else:
+                phsepoch_ld = time_to_longdouble(self.PEPOCH.quantity)
 
-        dt = (tbl['tdbld'] - phsepoch_ld)*u.day - delay
+            dt = (tbl['tdbld'] - phsepoch_ld)*u.day - delay
 
         return dt
 
diff --git a/pint/models/timing_model.py b/pint/models/timing_model.py
index ffda084af..e6e43cd9c 100644
--- a/pint/models/timing_model.py
+++ b/pint/models/timing_model.py
@@ -63,8 +63,8 @@ class TimingModel(object):
         A list for register the timing model component type name. For example,
         delay components will be register as 'DelayComponent'.
     top_level_params : list
-        A parameter name list for thoes parameters belong to the top timing
-        model class rather than a specific component.
+        A parameter name list for those parameters that belong to the top
+        timing model class rather than a specific component.
 
     Properties
     ----------
@@ -96,7 +96,7 @@ def __init__(self, name='', components=[]):
         self.add_param_from_top(strParameter(name="EPHEM",
             description="Ephemeris to use"), '')
         self.add_param_from_top(strParameter(name="UNITS",
-            description="Units (TDB assumed)"), '')
+            description="Units", value="TDB"), '')  # defaults to TDB if not given
 
         self.setup_components(components)
 
@@ -378,7 +378,7 @@ def add_param_from_top(self, param, target_component):
         else:
             if target_component not in list(self.components.keys()):
                 raise AttributeError("Can not find component '%s' in "
-                                     "timging model." % target_component)
+                                     "timing model." % target_component)
             self.components[target_component].add_param(param)
 
     def remove_param(self, param):
@@ -424,6 +424,7 @@ def get_params_of_type(self, param_type):
                 param_type.upper() == par_prefix.upper():
                 result.append(par.name)
         return result
+
     def get_prefix_mapping(self, prefix):
         """Get the index mapping for the prefix parameters.
            Parameter
@@ -564,8 +565,6 @@ def noise_model_basis_weight(self, toas):
             result.append(nf(toas)[1])
         return np.hstack((r for r in result))
 
-
-
     def get_barycentric_toas(self, toas, cutoff_component=''):
         """This is a convenient function for calculate the barycentric TOAs.
            Parameter
@@ -587,7 +586,14 @@ def get_barycentric_toas(self, toas, cutoff_component=''):
                 if cp.category == 'pulsar_system':
                     cutoff_component = cp.__class__.__name__
         corr = self.delay(toas, cutoff_component, False)
-        return tbl['tdbld'] * u.day - corr
+
+        timescale = 'TDB' if not hasattr(self, 'UNITS') else getattr(self, 'UNITS').value
+        if timescale == 'TDB':
+            return tbl['tdbld'] * u.day - corr
+        elif timescale == 'TCB':
+            return tbl['tcbld'] * u.day - corr
+        else:
+            raise ValueError("Unknown 'UNITS' value")
 
     def d_phase_d_toa(self, toas, sample_step=None):
         """Return the derivative of phase wrt TOA
@@ -791,12 +797,18 @@ def read_parfile(self, filename):
 
             k = l.split()
             name = k[0].upper()
-            
-            if name == 'UNITS' and len(k) > 1 and k[1] != 'TDB':
-                log.error("UNITS %s not yet supported by PINT" % k[1])
-                raise Exception("UNITS %s not yet supported by PINT" % k[1])
 
-            
+            if name == 'UNITS' and len(k) > 1:
+                if k[1] != 'TDB' and k[1] != 'TCB' and k[1] != 'SI':
+                    log.error("UNITS {} not yet supported by "
+                              "PINT".format(k[1]))
+                    raise Exception("UNITS {} not yet supported by "
+                                    "PINT".format(k[1]))
+
+                if k[1] == 'SI':
+                    # SI is just an alias for TCB so set to 'TCB'
+                    k[1] = 'TCB'
+
             if name in checked_param:
                 if name in repeat_param.keys():
                     repeat_param[name] += 1
@@ -1042,6 +1054,10 @@ def is_in_parfile(self,para_dict):
         """
         pNames_inpar = list(para_dict.keys())
         pNames_inModel = self.params
+        
+        # remove UNITS as many components contain it
+        if 'UNITS' in pNames_inModel:
+            del pNames_inModel[pNames_inModel.index('UNITS')]
 
         # For solar system Shapiro delay component
         if hasattr(self,'PLANET_SHAPIRO'):
@@ -1085,6 +1101,7 @@ def print_par(self,):
             result += getattr(self, p).as_parfile_line()
         return result
 
+
 class DelayComponent(Component):
     def __init__(self,):
         super(DelayComponent, self).__init__()
diff --git a/pint/models/wave.py b/pint/models/wave.py
index 5556a9616..205c9e384 100644
--- a/pint/models/wave.py
+++ b/pint/models/wave.py
@@ -33,6 +33,7 @@ def setup(self):
                                        "WAVE_OM is set.")
             else:
                 self.WAVEEPOCH = self.PEPOCH
+            self.WAVEEPOCH.scale = self.UNITS.value.lower()
 
         wave_terms = list(self.get_prefix_mapping_component('WAVE').keys())
         wave_terms.sort()
@@ -62,7 +63,10 @@ def wave_delay(self, toas, acc_delay=None):
                       range(1, self.num_wave_terms + 1)]
         wave_terms = [getattr(self, name) for name in wave_names]
         wave_om = self.WAVE_OM.quantity
-        time = self.barycentric_time = toas.table['tdbld'] * u.day
+        if self.UNITS.value == 'TCB':
+            time = self.barycentric_time = toas.table['tcbld'] * u.day
+        else:
+            time = self.barycentric_time = toas.table['tdbld'] * u.day
         for k, wave_term in enumerate(wave_terms):
             wave_a, wave_b = wave_term.quantity
             k = k + 1
diff --git a/pint/observatory/observatory.py b/pint/observatory/observatory.py
index c182908b3..50122b1f9 100644
--- a/pint/observatory/observatory.py
+++ b/pint/observatory/observatory.py
@@ -44,10 +44,12 @@ def __new__(cls, name, *args, **kwargs):
         cls._register(obs, name)
         return obs
 
-    def __init__(self,name,aliases=None,tt2tdb_mode='pint'):
+    def __init__(self, name, aliases=None, tt2tdb_mode='pint',
+                 tt2tcb_mode='pint'):
         if aliases is not None:
             Observatory._add_aliases(self,aliases)
         self.tt2tdb_mode = tt2tdb_mode
+        self.tt2tcb_mode = tt2tcb_mode
 
     @classmethod
     def _register(cls,obs,name):
@@ -171,7 +173,7 @@ def get_TDBs(self, t,  method='default', ephem=None, options=None, grp=None):
                     Also uses topocentric correction term if self.tt2tdbmethod is
                     pint.
                 - 'ephemeris': JPL ephemeris included TDB-TT correction.
-            ephme: str, optional
+            ephem: str, optional
                 The ephemeris to get he TDB-TT correction. Required for the
                 'ephemeris' method.
         """
@@ -237,6 +239,87 @@ def _get_TDB_ephem(self, t, ephem):
         """
         raise NotImplementedError
 
+    def get_TCBs(self, t,  method='default', ephem=None, options=None, grp=None):
+        """This is a high level function for converting TOAs to TCB time scale.
+            Different method can be applied to obtain the result. Current supported
+            methods are ['astropy', 'ephemeris']
+            Parameters
+            ----------
+            t: astropy.time.Time object
+                The time need for converting toas
+            method: str or callable, optional
+                Method of computing TCB
+
+                - 'default': Astropy time.Time object built-in converter, use FB90.
+                    Also uses topocentric correction term if self.tt2tcbmethod is
+                    pint.
+                - 'ephemeris': JPL ephemeris included TCB-TT correction.
+            ephem: str, optional
+                The ephemeris to get the TCB-TT correction. Required for the
+                'ephemeris' method.
+        """
+
+        if t.isscalar: t = Time([t])
+        if t.scale == 'tcb':
+            return t
+        # Check the method. This pattern is from numpy minize
+        if callable(method):
+            meth = "_custom"
+        else:
+            meth = method.lower()
+        if options is None:
+            options = {}
+        if meth == "_custom":
+            options = dict(options)
+            return method(t, **options)
+        if meth == 'default':
+            if self.tt2tcb_mode.lower().startswith('astropy'):
+                log.info('Doing astropy mode TCB conversion')
+                return self._get_TCB_astropy(t)
+            elif self.tt2tcb_mode.lower().startswith('pint'):
+                log.info('Doing PINT mode TCB conversion')
+                if ephem is None:
+                    raise ValueError("A ephemeris file should be provided to get"
+                                     " the TCB-TT corrections, or use tt2tcb_mode=astropy")
+                return self._get_TCB_PINT(t, ephem, grp)
+        elif meth == "ephemeris":
+            if ephem is None:
+                raise ValueError("A ephemeris file should be provided to get"
+                                 " the TCB-TT corrections.")
+            return self._get_TCB_ephem(t, ephem)
+        else:
+            raise ValueError("Unknown method '%s'." % method)
+
+    def _get_TCB_astropy(self, t):
+        return t.tcb
+
+    def _get_TCB_PINT(self, t, ephem, grp=None):
+        """Uses astropy.Time location to add the topocentric correction term to
+            the Time object. The topocentric correction is given as (r/c).(v/c),
+            with r equal to the geocentric position of the observer, v being the
+            barycentric velocity of the earth, and c being the speed of light.
+
+            The geocentric observer position can be obtained from Time object.
+            The barycentric velocity can be obtained using solar_system_ephemerides
+            objPosVel_wrt_SSB
+        """
+
+        #Add in correction term to t.tdb equal to r.v / c^2
+        vel = sse.objPosVel_wrt_SSB('earth', t, ephem).vel
+        pos = self.get_gcrs(t, ephem=ephem, grp=grp)
+        dnom = const.c * const.c 
+
+        corr = ((pos[0] * vel[0] + pos[1] * vel[1] + pos[2] * vel[2])/dnom).to(u.s)
+        log.debug('\tTopocentric Correction:\t%s' % corr)
+
+        return t.tcb + corr
+
+    def _get_TCB_ephem(self, t, ephem):
+        """This is a function that reads the ephem TCB-TT column. This column is
+            provided by DE4XXt version of ephemeris.
+        """
+        raise NotImplementedError
+
     def posvel(self,t,ephem):
         """Returns observatory position and velocity relative to solar system
         barycenter for the given times (astropy array-valued Time objects)."""
diff --git a/pint/observatory/topo_obs.py b/pint/observatory/topo_obs.py
index 415dde850..2dbbfb832 100644
--- a/pint/observatory/topo_obs.py
+++ b/pint/observatory/topo_obs.py
@@ -112,7 +112,8 @@ def __init__(self, name, tempo_code=None, itoa_code=None, aliases=None,
         for code in (tempo_code, itoa_code):
             if code is not None: aliases.append(code)
 
-        super(TopoObs,self).__init__(name,aliases=aliases, tt2tdb_mode='astropy')
+        super(TopoObs,self).__init__(name, aliases=aliases, tt2tdb_mode='astropy',
+                                     tt2tcb_mode='astropy')
 
     @property
     def clock_fullpath(self):
diff --git a/pint/toa.py b/pint/toa.py
index 88d67238d..4fa21f658 100644
--- a/pint/toa.py
+++ b/pint/toa.py
@@ -32,7 +32,7 @@
 
 def get_TOAs(timfile, ephem=None, include_bipm=True, bipm_version='BIPM2015',
              include_gps=True, planets=False, usepickle=False,
-             tdb_method="default"):
+             tdb_method="default", tcb_method="default"):
     """Convenience function to load and prepare TOAs for PINT use.
 
     Loads TOAs from a '.tim' file, applies clock corrections, computes
@@ -59,6 +59,8 @@ def get_TOAs(timfile, ephem=None, include_bipm=True, bipm_version='BIPM2015',
                                   bipm_version=bipm_version)
     if 'tdb' not in t.table.colnames:
         t.compute_TDBs(method=tdb_method, ephem=ephem)
+    if 'tcb' not in t.table.colnames:
+        t.compute_TCBs(method=tcb_method, ephem=ephem)
     if 'ssb_obs_pos' not in t.table.colnames:
         t.compute_posvels(ephem, planets)
     # Update pickle if needed:
@@ -101,7 +103,7 @@ def _check_pickle(toafilename, picklefilename=None):
 
 def get_TOAs_list(toa_list,ephem=None, include_bipm=True,
                   bipm_version='BIPM2015', include_gps=True, planets=False,
-                  tdb_method="default"):
+                  tdb_method="default", tcb_method="default"):
     """Load TOAs from a list of TOA objects.
 
     Compute the TDB time and observatory positions and velocity
@@ -119,6 +121,8 @@ def get_TOAs_list(toa_list,ephem=None, include_bipm=True,
                                   bipm_version=bipm_version)
     if 'tdb' not in t.table.colnames:
         t.compute_TDBs(method=tdb_method, ephem=ephem)
+    if 'tcb' not in t.table.colnames:
+        t.compute_TCBs(method=tcb_method, ephem=ephem)
     if 'ssb_obs_pos' not in t.table.colnames:
         t.compute_posvels(ephem, planets)
     return t
@@ -688,6 +692,7 @@ def adjust_TOAs(self, delta):
         # Changed high_precision from False to True to avoid self referential get_mjds()
         self.table['mjd_float'] = [t.mjd for t in self.get_mjds(high_precision=True)] * u.day
         self.compute_TDBs()
+        self.compute_TCBs()
         self.compute_posvels(self.ephem, self.planets)
 
     def write_TOA_file(self,filename,name='pint', format='Princeton'):
@@ -871,6 +876,72 @@ def compute_TDBs(self, method="default", ephem=None):
                 data=[utils.time_to_longdouble(t) for t in tdbs])
         self.table.add_columns([col_tdb, col_tdbld])
 
+    def compute_TCBs(self, method="default", ephem=None):
+        """Compute and add TCB and TCB long double columns to the TOA table.
+        This routine creates new columns 'tcb' and 'tcbld' in a TOA table
+        for TCB times, using the Observatory locations and IERS A Earth
+        rotation corrections for UT1.
+        """
+        log.info('Computing TCB columns.')
+        if 'tcb' in self.table.colnames:
+            log.info('tcb column already exists. Deleting...')
+            self.table.remove_column('tcb')
+        if 'tcbld' in self.table.colnames:
+            log.info('tcbld column already exists. Deleting...')
+            self.table.remove_column('tcbld')
+
+        if ephem is None:
+            if self.ephem is not None:
+                ephem = self.ephem
+            else:
+                log.warning('No ephemeris provided to TOAs object or compute_TCBs. Using DE421')
+                ephem = 'DE421'
+        else:
+            # If user specifies an ephemeris, make sure it is the same as the one already in the TOA object, to prevent mixing.
+            if (self.ephem is not None) and (ephem != self.ephem):
+                log.error('Ephemeris provided to compute_TCBs {0} is different than TOAs object ephemeris {1}!'.format(ephem,self.ephem))
+        self.ephem = ephem
+
+        # Compute in observatory groups
+        tcbs = numpy.zeros_like(self.table['mjd'])
+        for ii, key in enumerate(self.table.groups.keys):
+            grp = self.table.groups[ii]
+            obs = self.table.groups.keys[ii]['obs']
+            loind, hiind = self.table.groups.indices[ii:ii+2]
+            site = get_observatory(obs)
+            if isinstance(site, TopoObs):
+                # For TopoObs, it is safe to assume that all TOAs have same location
+                # I think we should report to astropy that initializing
+                # a Time from a list (or Column) of Times throws away the location information
+                grpmjds = time.Time(grp['mjd'], location=grp['mjd'][0].location)
+            else:
+                # Grab locations for each TOA
+                # It is crazy that I have to deconstruct the locations like
+                # this to build a single EarthLocation object with an array
+                # of locations contained in it.
+                # Is there a more efficient way to convert a list of EarthLocations
+                # into a single EarthLocation object with an array of values internally?
+                loclist = [t.location for t in grp['mjd']]
+                if loclist[0] is None:
+                    grpmjds = time.Time(grp['mjd'],location=None)
+                else:
+                    locs = EarthLocation(numpy.array([l.x.value for l in loclist])*u.m,
+                                         numpy.array([l.y.value for l in loclist])*u.m,
+                                         numpy.array([l.z.value for l in loclist])*u.m)
+                    grpmjds = time.Time(grp['mjd'],location=locs)
+
+            if isinstance(site, SpacecraftObs): #spacecraft-topocentric toas
+                grptcbs = site.get_TCBs(grpmjds, method=method, ephem=ephem, grp=grp)
+            else:
+                grptcbs = site.get_TCBs(grpmjds, method=method, ephem=ephem)
+            tcbs[loind:hiind] = numpy.asarray([t for t in grptcbs])
+
+        # Now add the new columns to the table
+        col_tcb = table.Column(name='tcb', data=tcbs)
+        col_tcbld = table.Column(name='tcbld',
+                data=[utils.time_to_longdouble(t) for t in tcbs])
+        self.table.add_columns([col_tcb, col_tcbld])
+
     def compute_posvels(self, ephem=None, planets=False):
         """Compute positions and velocities of the observatories and Earth.