diff --git a/pint/models/astrometry.py b/pint/models/astrometry.py
index c013f5991..2c9e8619b 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.
 
@@ -190,6 +189,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 +358,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..50086dbca 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,10 @@ 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]
+            DMEPOCH = tbl[self.UNITS.value.lower()+'ld'][0]
         else:
             DMEPOCH = self.DMEPOCH.value
-        dt = (tbl['tdbld'] - DMEPOCH) * u.day
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - 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 +158,10 @@ 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]
+            DMEPOCH = tbl[self.UNITS.value.lower()+'ld'][0]
         else:
             DMEPOCH = self.DMEPOCH.value
-        dt = (tbl['tdbld'] - DMEPOCH) * u.day
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - 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/glitch.py b/pint/models/glitch.py
index 75ce37a8e..d7ac10840 100644
--- a/pint/models/glitch.py
+++ b/pint/models/glitch.py
@@ -76,6 +76,7 @@ def setup(self):
             if not hasattr(self, 'GLEP_%d' % idx):
                 msg = 'Glicth Epoch is needed for Glicth %d.' % idx
                 raise MissingParameter("Glitch", 'GLEP_%d' % idx, msg)
+            getattr(self, 'GLEP_%d' % idx).scale = self.UNITS.value.lower()
             for param in self.glitch_prop:
                 if not hasattr(self, param + '%d' % idx):
                     param0 = getattr(self, param + '1')
@@ -121,7 +122,7 @@ def glitch_phase(self, toas, delay):
                 dF0 = getattr(self, "GLF0_%d" % idx).quantity
                 dF1 = getattr(self, "GLF1_%d" % idx).quantity
                 dF2 = getattr(self, "GLF2_%d" % idx).quantity
-                dt = (tbl['tdbld'] - eph) * u.day - delay
+                dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
                 dt = dt.to(u.second)
                 affected = dt > 0.0  # TOAs affected by glitch
                 # decay term
@@ -146,7 +147,7 @@ def d_phase_d_GLPH(self, toas, param, delay):
             raise ValueError("Can not calculate d_phase_d_GLPH with respect to %s." % param)
         eph = time_to_longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLPH = getattr(self, param)
-        dt = (tbl['tdbld'] - eph) * u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = numpy.where(dt > 0.0)[0]
         dpdGLPH = numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLPH.units
@@ -163,7 +164,7 @@ def d_phase_d_GLF0(self, toas, param, delay):
             raise ValueError("Can not calculate d_phase_d_GLF0 with respect to %s." % param)
         eph = time_to_longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLF0 = getattr(self, param)
-        dt = (tbl['tdbld'] - eph) * u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = numpy.where(dt > 0.0)[0]
         dpdGLF0 = numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLF0.units
@@ -179,7 +180,7 @@ def d_phase_d_GLF1(self, toas,  param, delay):
             raise ValueError("Can not calculate d_phase_d_GLF1 with respect to %s." % param)
         eph = time_to_longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLF1 = getattr(self, param)
-        dt = (tbl['tdbld'] - eph) * u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = numpy.where(dt > 0.0)[0]
         dpdGLF1 = numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLF1.units
@@ -195,7 +196,7 @@ def d_phase_d_GLF2(self, toas,  param, delay):
             raise ValueError("Can not calculate d_phase_d_GLF2 with respect to %s." % param)
         eph = time_to_longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLF2 = getattr(self, param)
-        dt = (tbl['tdbld'] - eph) * u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = numpy.where(dt > 0.0)[0]
         dpdGLF2 = numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLF2.units
@@ -213,7 +214,7 @@ def d_phase_d_GLF0D(self, toas, param, delay):
         eph = time_to_longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLF0D = getattr(self, param)
         tau = getattr(self, "GLTD_%d" % idv).quantity
-        dt = (tbl['tdbld'] - eph) * u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = numpy.where(dt > 0.0)[0]
         dpdGLF0D = numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLF0D.units
@@ -235,7 +236,7 @@ def d_phase_d_GLTD(self, toas, param, delay):
             return numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLTD.units
         glf0d = getattr(self, 'GLF0D_'+ids).quantity
         tau = par_GLTD.quantity
-        dt = (tbl['tdbld'] - eph) * u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = numpy.where(dt > 0.0)[0]
         dpdGLTD = numpy.zeros(len(tbl), dtype=numpy.longdouble) * u.cycle/par_GLTD.units
diff --git a/pint/models/noise_model.py b/pint/models/noise_model.py
index fd991233a..ba454adda 100644
--- a/pint/models/noise_model.py
+++ b/pint/models/noise_model.py
@@ -181,7 +181,7 @@ def ecorr_basis_weight_pair(self, toas):
 
         """
         tbl = toas.table
-        t = (tbl['tdbld'].quantity * u.day).to(u.s).value
+        t = (tbl[self.UNITS.value.lower()+'ld'].quantity * u.day).to(u.s).value
         ecorrs = self.get_ecorrs()
         Umats = []
         for ec in ecorrs:
@@ -263,7 +263,7 @@ def pl_rn_basis_weight_pair(self, toas):
 
         """
         tbl = toas.table
-        t = (tbl['tdbld'].quantity * u.day).to(u.s).value
+        t = (tbl[self.UNITS.value.lower()+'ld'].quantity * u.day).to(u.s).value
         amp, gam, nf = self.get_pl_vals()
         Fmat, f = create_fourier_design_matrix(t, nf)
         weight = powerlaw(f, amp, gam) * f[0]
diff --git a/pint/models/pulsar_binary.py b/pint/models/pulsar_binary.py
index e606f1c27..3c482714e 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,17 @@ 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
+        self.barycentric_time = tbl[self.UNITS.value.lower()+'ld'] * u.day - acc_delay
+        updates['psr_pos'] = self.ssb_to_psb_xyz_ICRS(epoch=tbl[self.UNITS.value.lower()+'ld'].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 +187,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/solar_system_shapiro.py b/pint/models/solar_system_shapiro.py
index 256df2df4..7fde460ab 100644
--- a/pint/models/solar_system_shapiro.py
+++ b/pint/models/solar_system_shapiro.py
@@ -78,7 +78,7 @@ def solar_system_shapiro_delay(self, toas, acc_delay=None):
             if key['obs'].lower() == 'barycenter':
                 log.debug("Skipping Shapiro delay for Barycentric TOAs")
                 continue
-            psr_dir = self.ssb_to_psb_xyz_ICRS(epoch=grp['tdbld'].astype(numpy.float64))
+            psr_dir = self.ssb_to_psb_xyz_ICRS(epoch=grp[self.UNITS.value.lower()+'ld'].astype(numpy.float64))
             delay[loind:hiind] += self.ss_obj_shapiro_delay(grp['obs_sun_pos'],
                                     psr_dir, self._ss_mass_sec['sun'])
             if self.PLANET_SHAPIRO.value:
diff --git a/pint/models/solar_wind_dispersion.py b/pint/models/solar_wind_dispersion.py
index 3f47d5467..3545be36b 100644
--- a/pint/models/solar_wind_dispersion.py
+++ b/pint/models/solar_wind_dispersion.py
@@ -42,7 +42,7 @@ def solar_wind_delay(self, toas, acc_delay=None):
                 bfreq = tbl['freq']
 
             rsa = tbl['obs_sun_pos'].quantity
-            pos = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(np.float64))
+            pos = self.ssb_to_psb_xyz_ICRS(epoch=tbl[self.UNITS.value.lower()+'ld'].astype(np.float64))
             r = np.sqrt(np.sum(rsa*rsa, axis=1))
             cos_theta = np.sum(rsa*pos, axis=1) / r
             ret =  const.au**2.0 * np.arccos(cos_theta) * DMconst * self.NE_SW.quantity / \
@@ -63,7 +63,7 @@ def d_delay_d_ne_sw(self, toas, param_name, acc_delay=None):
                 bfreq = tbl['freq']
 
             rsa = tbl['obs_sun_pos'].quantity
-            pos = self.ssb_to_psb_xyz_ICRS(epoch=tbl['tdbld'].astype(np.float64))
+            pos = self.ssb_to_psb_xyz_ICRS(epoch=tbl[self.UNITS.value.lower()+'ld'].astype(np.float64))
             r = np.sqrt(np.sum(rsa*rsa, axis=1))
             cos_theta = np.sum(rsa*pos, axis=1) / r
             
diff --git a/pint/models/spindown.py b/pint/models/spindown.py
index c1d4c0e70..4d64bcdef 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()
@@ -88,11 +92,11 @@ def get_dt(self, toas, delay):
         """
         tbl = toas.table
         if self.PEPOCH.value is None:
-            phsepoch_ld = time_to_longdouble(tbl['tdb'][0] - delay[0])
+            phsepoch_ld = time_to_longdouble(tbl[self.UNITS.value.lower()][0] - delay[0])
         else:
             phsepoch_ld = time_to_longdouble(self.PEPOCH.quantity)
 
-        dt = (tbl['tdbld'] - phsepoch_ld)*u.day - delay
+        dt = (tbl[self.UNITS.value.lower()+'ld'] - phsepoch_ld)*u.day - delay
 
         return dt
 
diff --git a/pint/models/timing_model.py b/pint/models/timing_model.py
index ffda084af..5eb831283 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,12 @@ 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 in ['TDB', 'TCB']:
+            return tbl[timescale.lower()+'ld'] * 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 +795,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 +1052,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 +1099,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..604c9022d 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,8 @@ 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
+        time = self.barycentric_time = toas.table[self.UNITS.value.lower()+'ld'] * 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..ff5e7af93 100644
--- a/pint/observatory/observatory.py
+++ b/pint/observatory/observatory.py
@@ -44,10 +44,46 @@ 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=None,
+                 tt2tcb_mode=None):
         if aliases is not None:
             Observatory._add_aliases(self,aliases)
         self.tt2tdb_mode = tt2tdb_mode
+        self.tt2tcb_mode = tt2tcb_mode
+
+    @property
+    def tt2tdb_mode(self):
+        return self._tt2tdb_mode
+
+    @tt2tdb_mode.setter
+    def tt2tdb_mode(self, val):
+        if val is None:
+            # set default to 'pint'
+            self._tt2tdb_mode = 'pint'
+        elif isinstance(val, str):
+            if val.lower() not in ['pint', 'astropy']:
+                raise ValueError("tt2tdb mode must be 'astropy' or 'pint'")
+            else:
+                self._tt2tdb_mode = val.lower()
+        else:
+            raise TypeError("tt2tdb mode is the wrong type")
+
+    @property
+    def tt2tcb_mode(self):
+        return self._tt2tcb_mode
+
+    @tt2tcb_mode.setter
+    def tt2tcb_mode(self, val):
+        if val is None:
+            # set default to 'pint'
+            self._tt2tcb_mode = 'pint'
+        elif isinstance(val, str):
+            if val.lower() not in ['pint', 'astropy']:
+                raise ValueError("tt2tcb mode must be 'astropy' or 'pint'")
+            else:
+                self._tt2tcb_mode = val.lower()
+        else:
+            raise TypeError("tt2tcb mode is the wrong type")
 
     @classmethod
     def _register(cls,obs,name):
@@ -156,7 +192,7 @@ def clock_corrections(self,t):
         # TOA metadata which may be necessary in some cases.
         raise NotImplementedError
 
-    def get_TDBs(self, t,  method='default', ephem=None, options=None, grp=None):
+    def get_TDBs(self, t,  method='default', tt2tdbmode=None, ephem=None, options=None):
         """This is a high level function for converting TOAs to TDB time scale.
             Different method can be applied to obtain the result. Current supported
             methods are ['astropy', 'ephemeris']
@@ -166,12 +202,15 @@ def get_TDBs(self, t,  method='default', ephem=None, options=None, grp=None):
                 The time need for converting toas
             method: str or callable, optional
                 Method of computing TDB
-
                 - 'default': Astropy time.Time object built-in converter, use FB90.
                     Also uses topocentric correction term if self.tt2tdbmethod is
                     pint.
                 - 'ephemeris': JPL ephemeris included TDB-TT correction.
-            ephme: str, optional
+            tt2tdbmode: str
+                If set this overrides the default tt2tdb_mode set on
+                initialization of the observatory. This can be 'astropy' or
+                'pint', which 'pint' computes the topocentric correction term.
+            ephem: str, optional
                 The ephemeris to get he TDB-TT correction. Required for the
                 'ephemeris' method.
         """
@@ -190,15 +229,19 @@ def get_TDBs(self, t,  method='default', ephem=None, options=None, grp=None):
             options = dict(options)
             return method(t, **options)
         if meth == 'default':
-            if self.tt2tdb_mode.lower().startswith('astropy'):
+            tt2tdb_mode = self.tt2tdb_mode.lower() if tt2tdbmode is None else tt2tdbmode.lower()
+
+            if tt2tdb_mode.startswith('astropy'):
                 log.info('Doing astropy mode TDB conversion')
                 return self._get_TDB_astropy(t)
-            elif self.tt2tdb_mode.lower().startswith('pint'):
+            elif tt2tdb_mode.startswith('pint'):
                 log.info('Doing PINT mode TDB conversion')
                 if ephem is None:
                     raise ValueError("A ephemeris file should be provided to get"
                                         " the TDB-TT corrections, or use tt2tdb_mode=astropy")
-                return self._get_TDB_PINT(t, ephem, grp)
+                return self._get_TDB_PINT(t, ephem)
+            else:
+                raise ValueError("The TDB-TT mode should be 'astropy' or 'pint'")
         elif meth == "ephemeris":
             if ephem is None:
                 raise ValueError("A ephemeris file should be provided to get"
@@ -210,7 +253,7 @@ def get_TDBs(self, t,  method='default', ephem=None, options=None, grp=None):
     def _get_TDB_astropy(self, t):
         return t.tdb
 
-    def _get_TDB_PINT(self, t, ephem, grp=None):
+    def _get_TDB_PINT(self, t, ephem):
         """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
@@ -223,7 +266,7 @@ def _get_TDB_PINT(self, t, ephem, grp=None):
 
         #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)
+        pos = self.get_gcrs(t, ephem=ephem)
         dnom = const.c * const.c 
 
         corr = ((pos[0] * vel[0] + pos[1] * vel[1] + pos[2] * vel[2])/dnom).to(u.s)
@@ -237,6 +280,94 @@ def _get_TDB_ephem(self, t, ephem):
         """
         raise NotImplementedError
 
+    def get_TCBs(self, t,  method='default', tt2tcbmode=None, 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.
+            tt2tcbmode: str
+                If set this overrides the default tt2tcb_mode set on
+                initialization of the observatory. This can be 'astropy' or
+                'pint', which 'pint' computes the topocentric correction term.
+            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':
+            tt2tcb_mode = self.tt2tcb_mode.lower() if tt2tcbmode is None else tt2tcbmode.lower()
+
+            if tt2tcb_mode.startswith('astropy'):
+                log.info('Doing astropy mode TCB conversion')
+                return self._get_TCB_astropy(t)
+            elif tt2tcb_mode.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)
+            else:
+                raise ValueError("The TCB-TT mode should be 'astropy' or 'pint'")
+        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)
+        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/scripts/pintempo.py b/pint/scripts/pintempo.py
index acf8b081a..9344ad04d 100755
--- a/pint/scripts/pintempo.py
+++ b/pint/scripts/pintempo.py
@@ -35,7 +35,7 @@ def main(argv=None):
     log.warning(m.params)
 
     log.info("Reading TOAs")
-    t = pint.toa.get_TOAs(args.timfile)
+    t = pint.toa.get_TOAs(args.timfile, ephem=m.EPHEM.value)
     prefit_resids = pint.residuals.resids(t, m).time_resids
 
     log.info("Fitting...")
diff --git a/pint/toa.py b/pint/toa.py
index 88d67238d..5db6143e6 100644
--- a/pint/toa.py
+++ b/pint/toa.py
@@ -32,7 +32,8 @@
 
 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", tt2tdb_mode=None,
+             tt2tcb_mode=None):
     """Convenience function to load and prepare TOAs for PINT use.
 
     Loads TOAs from a '.tim' file, applies clock corrections, computes
@@ -58,7 +59,11 @@ def get_TOAs(timfile, ephem=None, include_bipm=True, bipm_version='BIPM2015',
                                   include_bipm=include_bipm,
                                   bipm_version=bipm_version)
     if 'tdb' not in t.table.colnames:
-        t.compute_TDBs(method=tdb_method, ephem=ephem)
+        t.compute_TDBs(method=tdb_method, tt2tdb_mode=tt2tdb_mode,
+                       ephem=ephem)
+    if 'tcb' not in t.table.colnames:
+        t.compute_TCBs(method=tcb_method, tt2tcb_mode=tt2tcb_mode,
+                       ephem=ephem)
     if 'ssb_obs_pos' not in t.table.colnames:
         t.compute_posvels(ephem, planets)
     # Update pickle if needed:
@@ -101,7 +106,8 @@ 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", tt2tdb_mode=None,
+                  tt2tcb_mode=None):
     """Load TOAs from a list of TOA objects.
 
     Compute the TDB time and observatory positions and velocity
@@ -118,7 +124,9 @@ def get_TOAs_list(toa_list,ephem=None, include_bipm=True,
                                   include_bipm=include_bipm,
                                   bipm_version=bipm_version)
     if 'tdb' not in t.table.colnames:
-        t.compute_TDBs(method=tdb_method, ephem=ephem)
+        t.compute_TDBs(method=tdb_method, ephem=ephem, tt2tdb_mode=tt2tdb_mode)
+    if 'tcb' not in t.table.colnames:
+        t.compute_TCBs(method=tcb_method, ephem=ephem, tt2tcb_mode=tt2tcb_mode)
     if 'ssb_obs_pos' not in t.table.colnames:
         t.compute_posvels(ephem, planets)
     return t
@@ -688,6 +696,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'):
@@ -805,7 +814,7 @@ def apply_clock_corrections(self, include_bipm=True,
                                      'bipm_version':bipm_version,
                                      'include_gps':include_gps})
 
-    def compute_TDBs(self, method="default", ephem=None):
+    def compute_TDBs(self, method="default", ephem=None, tt2tdb_mode=None):
         """Compute and add TDB and TDB long double columns to the TOA table.
         This routine creates new columns 'tdb' and 'tdbld' in a TOA table
         for TDB times, using the Observatory locations and IERS A Earth
@@ -838,6 +847,7 @@ def compute_TDBs(self, method="default", ephem=None):
             obs = self.table.groups.keys[ii]['obs']
             loind, hiind = self.table.groups.indices[ii:ii+2]
             site = get_observatory(obs)
+            site.tt2tdb_mode = tt2tdb_mode
             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
@@ -871,6 +881,73 @@ 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, tt2tcb_mode=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)
+            site.tt2tcb_mode = tt2tcb_mode
+            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.
 
diff --git a/tests/datafile/test_TCB.par b/tests/datafile/test_TCB.par
new file mode 100644
index 000000000..bd4ac3c5f
--- /dev/null
+++ b/tests/datafile/test_TCB.par
@@ -0,0 +1,26 @@
+PSRJ     J1010+2020
+RAJ      10:10:54.21
+DECJ     20:20:01.05
+#PMRA     12.3
+#PMDEC    -9.8
+#POSEPOCH 56789
+#DM       101.1
+#DMEPOCH  56789
+F0       19.67297587421
+F1       -5.48585e-12
+PEPOCH   56789
+#BINARY   BT
+#ECC      0.0001
+#T0       57123
+#A1       23.4
+#PB       10.0
+#OM       98.7
+#GLEP_1   57000
+#GLPH_1   0
+#GLF0_1   8.7e-07
+#GLF1_1   -6.6e-15
+TZRMJD   56990
+TZRFRQ   1440
+TZRSITE  pks
+EPHEM    DE405
+UNITS    TDB
diff --git a/tests/datafile/test_TCB.simulate b/tests/datafile/test_TCB.simulate
new file mode 100644
index 000000000..e51c90556
--- /dev/null
+++ b/tests/datafile/test_TCB.simulate
@@ -0,0 +1,368 @@
+FORMAT 1
+ fake.rf 1440.00000000 56990.82474824677843728 0.00001 7 
+ fake.rf 1440.00000000 56991.82201806510150277 0.00001 7 
+ fake.rf 1440.00000000 56992.81928765426196648 0.00001 7 
+ fake.rf 1440.00000000 56993.81655704281042674 0.00001 7 
+ fake.rf 1440.00000000 56994.81382625881313686 0.00001 7 
+ fake.rf 1440.00000000 56995.81109591825106264 0.00001 7 
+ fake.rf 1440.00000000 56996.80836546035216728 0.00001 7 
+ fake.rf 1440.00000000 56997.80563550080942292 0.00001 7 
+ fake.rf 1440.00000000 56998.80290489054827319 0.00001 7 
+ fake.rf 1440.00000000 56999.80017424543495252 0.00001 7 
+ fake.rf 1440.00000000 57000.79744418151563679 0.00001 7 
+ fake.rf 1440.00000000 57001.79471355021483348 0.00001 7 
+ fake.rf 1440.00000000 57002.79198296795236089 0.00001 7 
+ fake.rf 1440.00000000 57003.78925246303314722 0.00001 7 
+ fake.rf 1440.00000000 57004.78652206389039492 0.00001 7 
+ fake.rf 1440.00000000 57005.78379179906246321 0.00001 7 
+ fake.rf 1440.00000000 57006.78106110890541203 0.00001 7 
+ fake.rf 1440.00000000 57007.77833061041031115 0.00001 7 
+ fake.rf 1440.00000000 57008.77560033238588488 0.00001 7 
+ fake.rf 1440.00000000 57009.77287030373106091 0.00001 7 
+ fake.rf 1440.00000000 57010.77013937686564304 0.00001 7 
+ fake.rf 1440.00000000 57011.76740934560760010 0.00001 7 
+ fake.rf 1440.00000000 57012.76467847403396405 0.00001 7 
+ fake.rf 1440.00000000 57013.76194855563046815 0.00001 7 
+ fake.rf 1440.00000000 57014.75921785375821571 0.00001 7 
+ fake.rf 1440.00000000 57015.75648757247618548 0.00001 7 
+ fake.rf 1440.00000000 57016.75375715012821942 0.00001 7 
+ fake.rf 1440.00000000 57017.75102661233427170 0.00001 7 
+ fake.rf 1440.00000000 57018.74829598368354411 0.00001 7 
+ fake.rf 1440.00000000 57019.74556587606709002 0.00001 7 
+ fake.rf 1440.00000000 57020.74283513565772452 0.00001 7 
+ fake.rf 1440.00000000 57021.74010496100964573 0.00001 7 
+ fake.rf 1440.00000000 57022.73737419693858186 0.00001 7 
+ fake.rf 1440.00000000 57023.73464404090834279 0.00001 7 
+ fake.rf 1440.00000000 57024.73191333688029658 0.00001 7 
+ fake.rf 1440.00000000 57025.72918328168281477 0.00001 7 
+ fake.rf 1440.00000000 57026.72645271881614448 0.00001 7 
+ fake.rf 1440.00000000 57027.72372225650813959 0.00001 7 
+ fake.rf 1440.00000000 57028.72099191461168388 0.00001 7 
+ fake.rf 1440.00000000 57029.71826112461666369 0.00001 7 
+ fake.rf 1440.00000000 57030.71553108280073729 0.00001 7 
+ fake.rf 1440.00000000 57031.71280063219183987 0.00001 7 
+ fake.rf 1440.00000000 57032.71006979227609435 0.00001 7 
+ fake.rf 1440.00000000 57033.70733975898214041 0.00001 7 
+ fake.rf 1440.00000000 57034.70460937491583309 0.00001 7 
+ fake.rf 1440.00000000 57035.70187865909457869 0.00001 7 
+ fake.rf 1440.00000000 57036.69914821865699395 0.00001 7 
+ fake.rf 1440.00000000 57037.69641807227216290 0.00001 7 
+ fake.rf 1440.00000000 57038.69368765011937938 0.00001 7 
+ fake.rf 1440.00000000 57039.69095697043252358 0.00001 7 
+ fake.rf 1440.00000000 57040.68822663940961704 0.00001 7 
+ fake.rf 1440.00000000 57041.68549608616237734 0.00001 7 
+ fake.rf 1440.00000000 57042.68276591568706024 0.00001 7 
+ fake.rf 1440.00000000 57043.68003555542286520 0.00001 7 
+ fake.rf 1440.00000000 57044.67730501993095160 0.00001 7 
+ fake.rf 1440.00000000 57045.67457432247167048 0.00001 7 
+ fake.rf 1440.00000000 57046.67184406330909852 0.00001 7 
+ fake.rf 1440.00000000 57047.66911366491628144 0.00001 7 
+ fake.rf 1440.00000000 57048.66638313706305041 0.00001 7 
+ fake.rf 1440.00000000 57049.66365248870872051 0.00001 7 
+ fake.rf 1440.00000000 57050.66092231648193334 0.00001 7 
+ fake.rf 1440.00000000 57051.65819203989390473 0.00001 7 
+ fake.rf 1440.00000000 57052.65546166635476055 0.00001 7 
+ fake.rf 1440.00000000 57053.65273120293800702 0.00001 7 
+ fake.rf 1440.00000000 57054.65000065642915317 0.00001 7 
+ fake.rf 1440.00000000 57055.64727003336199118 0.00001 7 
+ fake.rf 1440.00000000 57056.64453992835570517 0.00001 7 
+ fake.rf 1440.00000000 57057.64180917087610112 0.00001 7 
+ fake.rf 1440.00000000 57058.63907894343209648 0.00001 7 
+ fake.rf 1440.00000000 57059.63634807505798818 0.00001 7 
+ fake.rf 1440.00000000 57060.63361774782078228 0.00001 7 
+ fake.rf 1440.00000000 57061.63088737859302313 0.00001 7 
+ fake.rf 1440.00000000 57062.62815697229916978 0.00001 7 
+ fake.rf 1440.00000000 57063.62542653358397260 0.00001 7 
+ fake.rf 1440.00000000 57064.62269606679950584 0.00001 7 
+ fake.rf 1440.00000000 57065.61996557599091418 0.00001 7 
+ fake.rf 1440.00000000 57066.61723506487265567 0.00001 7 
+ fake.rf 1440.00000000 57067.61450512510856825 0.00001 7 
+ fake.rf 1440.00000000 57068.61177458292954512 0.00001 7 
+ fake.rf 1440.00000000 57069.60904402893933707 0.00001 7 
+ fake.rf 1440.00000000 57070.60631346462255209 0.00001 7 
+ fake.rf 1440.00000000 57071.60358289043410096 0.00001 7 
+ fake.rf 1440.00000000 57072.60085289391258101 0.00001 7 
+ fake.rf 1440.00000000 57073.59812229623100421 0.00001 7 
+ fake.rf 1440.00000000 57074.59539168218080007 0.00001 7 
+ fake.rf 1440.00000000 57075.59266163533512994 0.00001 7 
+ fake.rf 1440.00000000 57076.58993097328163913 0.00001 7 
+ fake.rf 1440.00000000 57077.58720027786698381 0.00001 7 
+ fake.rf 1440.00000000 57078.58447013041608997 0.00001 7 
+ fake.rf 1440.00000000 57079.58173993520653866 0.00001 7 
+ fake.rf 1440.00000000 57080.57900909622576435 0.00001 7 
+ fake.rf 1440.00000000 57081.57627878223720685 0.00001 7 
+ fake.rf 1440.00000000 57082.57354839677957514 0.00001 7 
+ fake.rf 1440.00000000 57083.57081793152965332 0.00001 7 
+ fake.rf 1440.00000000 57084.56808737797342701 0.00001 7 
+ fake.rf 1440.00000000 57085.56535731575814197 0.00001 7 
+ fake.rf 1440.00000000 57086.56262655932592764 0.00001 7 
+ fake.rf 1440.00000000 57087.55989627631699790 0.00001 7 
+ fake.rf 1440.00000000 57088.55716586916413746 0.00001 7 
+ fake.rf 1440.00000000 57089.55443532845691479 0.00001 7 
+ fake.rf 1440.00000000 57090.55170464458324275 0.00001 7 
+ fake.rf 1440.00000000 57091.54897439605806042 0.00001 7 
+ fake.rf 1440.00000000 57092.54624398445398370 0.00001 7 
+ fake.rf 1440.00000000 57093.54351339941338850 0.00001 7 
+ fake.rf 1440.00000000 57094.54078321861362610 0.00001 7 
+ fake.rf 1440.00000000 57095.53805284265682474 0.00001 7 
+ fake.rf 1440.00000000 57096.53532226003633454 0.00001 7 
+ fake.rf 1440.00000000 57097.53259204701807761 0.00001 7 
+ fake.rf 1440.00000000 57098.52986160241573899 0.00001 7 
+ fake.rf 1440.00000000 57099.52713091245715660 0.00001 7 
+ fake.rf 1440.00000000 57100.52440055060904101 0.00001 7 
+ fake.rf 1440.00000000 57101.52166991239362659 0.00001 7 
+ fake.rf 1440.00000000 57102.51893956882526382 0.00001 7 
+ fake.rf 1440.00000000 57103.51620950148709355 0.00001 7 
+ fake.rf 1440.00000000 57104.51347910277480580 0.00001 7 
+ fake.rf 1440.00000000 57105.51074835289218967 0.00001 7 
+ fake.rf 1440.00000000 57106.50801782009445517 0.00001 7 
+ fake.rf 1440.00000000 57107.50528748413348268 0.00001 7 
+ fake.rf 1440.00000000 57108.50255732470939662 0.00001 7 
+ fake.rf 1440.00000000 57109.49982673314420367 0.00001 7 
+ fake.rf 1440.00000000 57110.49709627747996521 0.00001 7 
+ fake.rf 1440.00000000 57111.49436593738313661 0.00001 7 
+ fake.rf 1440.00000000 57112.49163510413378830 0.00001 7 
+ fake.rf 1440.00000000 57113.48890493409028224 0.00001 7 
+ fake.rf 1440.00000000 57114.48617423012598238 0.00001 7 
+ fake.rf 1440.00000000 57115.48344414857082896 0.00001 7 
+ fake.rf 1440.00000000 57116.48071349227289772 0.00001 7 
+ fake.rf 1440.00000000 57117.47798341754286966 0.00001 7 
+ fake.rf 1440.00000000 57118.47525272718723954 0.00001 7 
+ fake.rf 1440.00000000 57119.47252257744231585 0.00001 7 
+ fake.rf 1440.00000000 57120.46979177097041358 0.00001 7 
+ fake.rf 1440.00000000 57121.46706146377923830 0.00001 7 
+ fake.rf 1440.00000000 57122.46433104655453761 0.00001 7 
+ fake.rf 1440.00000000 57123.46160049806836412 0.00001 7 
+ fake.rf 1440.00000000 57124.45887038508634959 0.00001 7 
+ fake.rf 1440.00000000 57125.45613950875284104 0.00001 7 
+ fake.rf 1440.00000000 57126.45340902307911435 0.00001 7 
+ fake.rf 1440.00000000 57127.45067890455610993 0.00001 7 
+ fake.rf 1440.00000000 57128.44794854040256382 0.00001 7 
+ fake.rf 1440.00000000 57129.44521790525256577 0.00001 7 
+ fake.rf 1440.00000000 57130.44248756118716770 0.00001 7 
+ fake.rf 1440.00000000 57131.43975689274062901 0.00001 7 
+ fake.rf 1440.00000000 57132.43702646064068773 0.00001 7 
+ fake.rf 1440.00000000 57133.43429623693887009 0.00001 7 
+ fake.rf 1440.00000000 57134.43156560524821330 0.00001 7 
+ fake.rf 1440.00000000 57135.42883512605483531 0.00001 7 
+ fake.rf 1440.00000000 57136.42610477167707472 0.00001 7 
+ fake.rf 1440.00000000 57137.42337451467974674 0.00001 7 
+ fake.rf 1440.00000000 57138.42064432787096706 0.00001 7 
+ fake.rf 1440.00000000 57139.41791359590109067 0.00001 7 
+ fake.rf 1440.00000000 57140.41518346876954126 0.00001 7 
+ fake.rf 1440.00000000 57141.41245274313939717 0.00001 7 
+ fake.rf 1440.00000000 57142.40972256938984231 0.00001 7 
+ fake.rf 1440.00000000 57143.40699174456038634 0.00001 7 
+ fake.rf 1440.00000000 57144.40426141941808424 0.00001 7 
+ fake.rf 1440.00000000 57145.40153097977729502 0.00001 7 
+ fake.rf 1440.00000000 57146.39880040002007178 0.00001 7 
+ fake.rf 1440.00000000 57147.39607024306716099 0.00001 7 
+ fake.rf 1440.00000000 57148.39333989517342260 0.00001 7 
+ fake.rf 1440.00000000 57149.39060933100360273 0.00001 7 
+ fake.rf 1440.00000000 57150.38787911353901805 0.00001 7 
+ fake.rf 1440.00000000 57151.38514862881413592 0.00001 7 
+ fake.rf 1440.00000000 57152.38241785095470604 0.00001 7 
+ fake.rf 1440.00000000 57153.37968734206395993 0.00001 7 
+ fake.rf 1440.00000000 57154.37695707534511058 0.00001 7 
+ fake.rf 1440.00000000 57155.37422702339852165 0.00001 7 
+ fake.rf 1440.00000000 57156.37149656978159129 0.00001 7 
+ fake.rf 1440.00000000 57157.36876568577820734 0.00001 7 
+ fake.rf 1440.00000000 57158.36603551885218977 0.00001 7 
+ fake.rf 1440.00000000 57159.36330486247270244 0.00001 7 
+ fake.rf 1440.00000000 57160.36057486340225253 0.00001 7 
+ fake.rf 1440.00000000 57161.35784431483504520 0.00001 7 
+ fake.rf 1440.00000000 57162.35511377530802690 0.00001 7 
+ fake.rf 1440.00000000 57163.35238321531575068 0.00001 7 
+ fake.rf 1440.00000000 57164.34965260579076229 0.00001 7 
+ fake.rf 1440.00000000 57165.34692250652519618 0.00001 7 
+ fake.rf 1440.00000000 57166.34419171262252135 0.00001 7 
+ fake.rf 1440.00000000 57167.34146137317316061 0.00001 7 
+ fake.rf 1440.00000000 57168.33873087252517919 0.00001 7 
+ fake.rf 1440.00000000 57169.33600077219014679 0.00001 7 
+ fake.rf 1440.00000000 57170.33327045728925953 0.00001 7 
+ fake.rf 1440.00000000 57171.33053990170782654 0.00001 7 
+ fake.rf 1440.00000000 57172.32780966809857404 0.00001 7 
+ fake.rf 1440.00000000 57173.32507914272756011 0.00001 7 
+ fake.rf 1440.00000000 57174.32234830061952735 0.00001 7 
+ fake.rf 1440.00000000 57175.31961829390883167 0.00001 7 
+ fake.rf 1440.00000000 57176.31688792146841749 0.00001 7 
+ fake.rf 1440.00000000 57177.31415715912092779 0.00001 7 
+ fake.rf 1440.00000000 57178.31142657112490824 0.00001 7 
+ fake.rf 1440.00000000 57179.30869613326469292 0.00001 7 
+ fake.rf 1440.00000000 57180.30596582109366466 0.00001 7 
+ fake.rf 1440.00000000 57181.30323560981777220 0.00001 7 
+ fake.rf 1440.00000000 57182.30050488583250967 0.00001 7 
+ fake.rf 1440.00000000 57183.29777480022211833 0.00001 7 
+ fake.rf 1440.00000000 57184.29504415009731488 0.00001 7 
+ fake.rf 1440.00000000 57185.29231349733948520 0.00001 7 
+ fake.rf 1440.00000000 57186.28958340356430412 0.00001 7 
+ fake.rf 1440.00000000 57187.28685266512862739 0.00001 7 
+ fake.rf 1440.00000000 57188.28412243197474751 0.00001 7 
+ fake.rf 1440.00000000 57189.28139208914294045 0.00001 7 
+ fake.rf 1440.00000000 57190.27866161045896121 0.00001 7 
+ fake.rf 1440.00000000 57191.27593097027883573 0.00001 7 
+ fake.rf 1440.00000000 57192.27320073194664118 0.00001 7 
+ fake.rf 1440.00000000 57193.27047028269020856 0.00001 7 
+ fake.rf 1440.00000000 57194.26773959875565012 0.00001 7 
+ fake.rf 1440.00000000 57195.26500924538061099 0.00001 7 
+ fake.rf 1440.00000000 57196.26227920000711791 0.00001 7 
+ fake.rf 1440.00000000 57197.25954826387802044 0.00001 7 
+ fake.rf 1440.00000000 57198.25681818065502071 0.00001 7 
+ fake.rf 1440.00000000 57199.25408775264335404 0.00001 7 
+ fake.rf 1440.00000000 57200.25135695942677572 0.00001 7 
+ fake.rf 1440.00000000 57201.24862695786059419 0.00001 7 
+ fake.rf 1440.00000000 57202.24589655181621950 0.00001 7 
+ fake.rf 1440.00000000 57203.24316572239561296 0.00001 7 
+ fake.rf 1440.00000000 57204.24043523285127222 0.00001 7 
+ fake.rf 1440.00000000 57205.23770526662140767 0.00001 7 
+ fake.rf 1440.00000000 57206.23497482229969791 0.00001 7 
+ fake.rf 1440.00000000 57207.23224447018639438 0.00001 7 
+ fake.rf 1440.00000000 57208.22951360364263351 0.00001 7 
+ fake.rf 1440.00000000 57209.22678338081665927 0.00001 7 
+ fake.rf 1440.00000000 57210.22405319437213578 0.00001 7 
+ fake.rf 1440.00000000 57211.22132243662387907 0.00001 7 
+ fake.rf 1440.00000000 57212.21859226472768967 0.00001 7 
+ fake.rf 1440.00000000 57213.21586148224117707 0.00001 7 
+ fake.rf 1440.00000000 57214.21313124621691770 0.00001 7 
+ fake.rf 1440.00000000 57215.21040094879483817 0.00001 7 
+ fake.rf 1440.00000000 57216.20767057082290208 0.00001 7 
+ fake.rf 1440.00000000 57217.20494009362570864 0.00001 7 
+ fake.rf 1440.00000000 57218.20220949911987418 0.00001 7 
+ fake.rf 1440.00000000 57219.19947876989915869 0.00001 7 
+ fake.rf 1440.00000000 57220.19674847764921211 0.00001 7 
+ fake.rf 1440.00000000 57221.19401801808714580 0.00001 7 
+ fake.rf 1440.00000000 57222.19128796440257290 0.00001 7 
+ fake.rf 1440.00000000 57223.18855712543703262 0.00001 7 
+ fake.rf 1440.00000000 57224.18582666416810767 0.00001 7 
+ fake.rf 1440.00000000 57225.18309656752269987 0.00001 7 
+ fake.rf 1440.00000000 57226.18036623471621027 0.00001 7 
+ fake.rf 1440.00000000 57227.17763565394961489 0.00001 7 
+ fake.rf 1440.00000000 57228.17490540239573349 0.00001 7 
+ fake.rf 1440.00000000 57229.17217488112762425 0.00001 7 
+ fake.rf 1440.00000000 57230.16944408017004520 0.00001 7 
+ fake.rf 1440.00000000 57231.16671357846125900 0.00001 7 
+ fake.rf 1440.00000000 57232.16398336707794314 0.00001 7 
+ fake.rf 1440.00000000 57233.16125284912355298 0.00001 7 
+ fake.rf 1440.00000000 57234.15852260461240419 0.00001 7 
+ fake.rf 1440.00000000 57235.15579203686940346 0.00001 7 
+ fake.rf 1440.00000000 57236.15306172573456323 0.00001 7 
+ fake.rf 1440.00000000 57237.15033107403378665 0.00001 7 
+ fake.rf 1440.00000000 57238.14760066093270652 0.00001 7 
+ fake.rf 1440.00000000 57239.14487047696728439 0.00001 7 
+ fake.rf 1440.00000000 57240.14213992420591026 0.00001 7 
+ fake.rf 1440.00000000 57241.13940958147928484 0.00001 7 
+ fake.rf 1440.00000000 57242.13667885120131729 0.00001 7 
+ fake.rf 1440.00000000 57243.13394890125039893 0.00001 7 
+ fake.rf 1440.00000000 57244.13121795836891081 0.00001 7 
+ fake.rf 1440.00000000 57245.12848778001116301 0.00001 7 
+ fake.rf 1440.00000000 57246.12575718266680980 0.00001 7 
+ fake.rf 1440.00000000 57247.12302674862520036 0.00001 7 
+ fake.rf 1440.00000000 57248.12029647264227350 0.00001 7 
+ fake.rf 1440.00000000 57249.11756576195469037 0.00001 7 
+ fake.rf 1440.00000000 57250.11483578961291840 0.00001 7 
+ fake.rf 1440.00000000 57251.11210537611639992 0.00001 7 
+ fake.rf 1440.00000000 57252.10937451939823006 0.00001 7 
+ fake.rf 1440.00000000 57253.10664439478508569 0.00001 7 
+ fake.rf 1440.00000000 57254.10391382497633117 0.00001 7 
+ fake.rf 1440.00000000 57255.10118339834330570 0.00001 7 
+ fake.rf 1440.00000000 57256.09845311558003544 0.00001 7 
+ fake.rf 1440.00000000 57257.09572238968852886 0.00001 7 
+ fake.rf 1440.00000000 57258.09299239926487957 0.00001 7 
+ fake.rf 1440.00000000 57259.09026197016670778 0.00001 7 
+ fake.rf 1440.00000000 57260.08753110540921938 0.00001 7 
+ fake.rf 1440.00000000 57261.08480098504233169 0.00001 7 
+ fake.rf 1440.00000000 57262.08207043604262765 0.00001 7 
+ fake.rf 1440.00000000 57263.07934005038140413 0.00001 7 
+ fake.rf 1440.00000000 57264.07660924320334672 0.00001 7 
+ fake.rf 1440.00000000 57265.07387919428828127 0.00001 7 
+ fake.rf 1440.00000000 57266.07114872977169284 0.00001 7 
+ fake.rf 1440.00000000 57267.06841844049096224 0.00001 7 
+ fake.rf 1440.00000000 57268.06568774059424243 0.00001 7 
+ fake.rf 1440.00000000 57269.06295722106677815 0.00001 7 
+ fake.rf 1440.00000000 57270.06022688500947737 0.00001 7 
+ fake.rf 1440.00000000 57271.05749673613633988 0.00001 7 
+ fake.rf 1440.00000000 57272.05476619057472476 0.00001 7 
+ fake.rf 1440.00000000 57273.05203584187026422 0.00001 7 
+ fake.rf 1440.00000000 57274.04930510775883334 0.00001 7 
+ fake.rf 1440.00000000 57275.04657458342537524 0.00001 7 
+ fake.rf 1440.00000000 57276.04384427656298229 0.00001 7 
+ fake.rf 1440.00000000 57277.04111360736748892 0.00001 7 
+ fake.rf 1440.00000000 57278.03838317345126185 0.00001 7 
+ fake.rf 1440.00000000 57279.03565298489357716 0.00001 7 
+ fake.rf 1440.00000000 57280.03292246421952072 0.00001 7 
+ fake.rf 1440.00000000 57281.03019221128980831 0.00001 7 
+ fake.rf 1440.00000000 57282.02746165004108292 0.00001 7 
+ fake.rf 1440.00000000 57283.02473138167335165 0.00001 7 
+ fake.rf 1440.00000000 57284.02200083139569387 0.00001 7 
+ fake.rf 1440.00000000 57285.01927060160479854 0.00001 7 
+ fake.rf 1440.00000000 57286.01654011863154281 0.00001 7 
+ fake.rf 1440.00000000 57287.01380939759170730 0.00001 7 
+ fake.rf 1440.00000000 57288.01107904232643975 0.00001 7 
+ fake.rf 1440.00000000 57289.00834848035344393 0.00001 7 
+ fake.rf 1440.00000000 57290.00561831595653572 0.00001 7 
+ fake.rf 1440.00000000 57291.00288797678801700 0.00001 7 
+ fake.rf 1440.00000000 57292.00015747857959880 0.00001 7 
+ fake.rf 1440.00000000 57292.99742683670141830 0.00001 7 
+ fake.rf 1440.00000000 57293.99469665442014943 0.00001 7 
+ fake.rf 1440.00000000 57294.99196635811814460 0.00001 7 
+ fake.rf 1440.00000000 57295.98923537411346629 0.00001 7 
+ fake.rf 1440.00000000 57296.98650548211231381 0.00001 7 
+ fake.rf 1440.00000000 57297.98377493254149684 0.00001 7 
+ fake.rf 1440.00000000 57298.98104432966871613 0.00001 7 
+ fake.rf 1440.00000000 57299.97831369024260439 0.00001 7 
+ fake.rf 1440.00000000 57300.97558362011129063 0.00001 7 
+ fake.rf 1440.00000000 57301.97285296105766506 0.00001 7 
+ fake.rf 1440.00000000 57302.97012290879930063 0.00001 7 
+ fake.rf 1440.00000000 57303.96739230664508824 0.00001 7 
+ fake.rf 1440.00000000 57304.96466176348654287 0.00001 7 
+ fake.rf 1440.00000000 57305.96193130061832832 0.00001 7 
+ fake.rf 1440.00000000 57306.95920094000327794 0.00001 7 
+ fake.rf 1440.00000000 57307.95647070424731240 0.00001 7 
+ fake.rf 1440.00000000 57308.95374002828503279 0.00001 7 
+ fake.rf 1440.00000000 57309.95100952419700135 0.00001 7 
+ fake.rf 1440.00000000 57310.94827921631057777 0.00001 7 
+ fake.rf 1440.00000000 57311.94554854115097697 0.00001 7 
+ fake.rf 1440.00000000 57312.94281811223243395 0.00001 7 
+ fake.rf 1440.00000000 57313.94008795514355725 0.00001 7 
+ fake.rf 1440.00000000 57314.93735750747520541 0.00001 7 
+ fake.rf 1440.00000000 57315.93462679529532977 0.00001 7 
+ fake.rf 1440.00000000 57316.93189643303754721 0.00001 7 
+ fake.rf 1440.00000000 57317.92916585852914224 0.00001 7 
+ fake.rf 1440.00000000 57318.92643568596198733 0.00001 7 
+ fake.rf 1440.00000000 57319.92370535262528719 0.00001 7 
+ fake.rf 1440.00000000 57320.92097488363430458 0.00001 7 
+ fake.rf 1440.00000000 57321.91824430362532183 0.00001 7 
+ fake.rf 1440.00000000 57322.91551422512785408 0.00001 7 
+ fake.rf 1440.00000000 57323.91278349566498918 0.00001 7 
+ fake.rf 1440.00000000 57324.91005331597205696 0.00001 7 
+ fake.rf 1440.00000000 57325.90732253404712182 0.00001 7 
+ fake.rf 1440.00000000 57326.90459235156415474 0.00001 7 
+ fake.rf 1440.00000000 57327.90186161776680507 0.00001 7 
+ fake.rf 1440.00000000 57328.89913153570139670 0.00001 7 
+ fake.rf 1440.00000000 57329.89640095598140235 0.00001 7 
+ fake.rf 1440.00000000 57330.89367049466433812 0.00001 7 
+ fake.rf 1440.00000000 57331.89094018009364007 0.00001 7 
+ fake.rf 1440.00000000 57332.88820945285193886 0.00001 7 
+ fake.rf 1440.00000000 57333.88547893054593985 0.00001 7 
+ fake.rf 1440.00000000 57334.88274864295358668 0.00001 7 
+ fake.rf 1440.00000000 57335.88001803199797735 0.00001 7 
+ fake.rf 1440.00000000 57336.87728771653977589 0.00001 7 
+ fake.rf 1440.00000000 57337.87455772751810557 0.00001 7 
+ fake.rf 1440.00000000 57338.87182691962787118 0.00001 7 
+ fake.rf 1440.00000000 57339.86909650089430102 0.00001 7 
+ fake.rf 1440.00000000 57340.86636591467999935 0.00001 7 
+ fake.rf 1440.00000000 57341.86363578095915372 0.00001 7 
+ fake.rf 1440.00000000 57342.86090554312374579 0.00001 7 
+ fake.rf 1440.00000000 57343.85817464442280667 0.00001 7 
+ fake.rf 1440.00000000 57344.85544429266866473 0.00001 7 
+ fake.rf 1440.00000000 57345.85271393048631339 0.00001 7 
+ fake.rf 1440.00000000 57346.84998358829571075 0.00001 7 
+ fake.rf 1440.00000000 57347.84725329594654397 0.00001 7 
+ fake.rf 1440.00000000 57348.84452249439187810 0.00001 7 
+ fake.rf 1440.00000000 57349.84179238879033846 0.00001 7 
+ fake.rf 1440.00000000 57350.83906183069199614 0.00001 7 
+ fake.rf 1440.00000000 57351.83633143613122840 0.00001 7 
+ fake.rf 1440.00000000 57352.83360123277189047 0.00001 7 
+ fake.rf 1440.00000000 57353.83087066014566346 0.00001 7 
+ fake.rf 1440.00000000 57354.82814033466054937 0.00001 7 
+ fake.rf 1440.00000000 57355.82540969661322450 0.00001 7 
diff --git a/tests/test_pintempo2tcb.py b/tests/test_pintempo2tcb.py
new file mode 100644
index 000000000..882cdcd4d
--- /dev/null
+++ b/tests/test_pintempo2tcb.py
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+from __future__ import division, print_function
+import sys, os
+try:
+    from StringIO import StringIO
+except ImportError:
+    from io import StringIO
+import unittest
+import numpy as np
+import pint.scripts.pintempo as pintempo
+from pinttestdata import testdir, datadir
+
+parfile = os.path.join(datadir, 'test_TCB.par')
+timfile = os.path.join(datadir, 'test_TCB.simulate')
+
+"""
+Note: test_TCB.simulate has been produced using TEMPO2, with the fake plugin,
+using:
+
+tempo2 -gr fake -start 56990 -end 57355 -f test_TCB.par -ndobs 1 -nobsd 1 -ha 6 -randha n -rms 0.00000001
+"""
+
+class TestPintempo2tcb(unittest.TestCase):
+
+    def test_result(self):
+        saved_stdout, pintempo.sys.stdout = pintempo.sys.stdout, StringIO('_')
+        cmd = '{0} {1}'.format(parfile,timfile)
+        pintempo.main(cmd.split())
+        lines = pintempo.sys.stdout.getvalue()
+        v = 999.0
+        for l in lines.split('\n'):
+            if l.startswith('RMS in time is'):
+                v = float(l.split()[4])
+        # Check that RMS is less than 10 microseconds
+        from astropy import log
+        log.warning('%f' % v)
+        print(v)
+        self.assertTrue(v<10.0)
+        pintempo.sys.stdout = saved_stdout
+
+if __name__ == '__main__':
+    unittest.main()