diff --git a/bin/compare_centroids.py b/bin/compare_centroids.py
new file mode 100644
index 0000000..ae3a62f
--- /dev/null
+++ b/bin/compare_centroids.py
@@ -0,0 +1,29 @@
+import argparse
+from desc.twinkles import CentroidValidator
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser("Read in an InstanceCatalog and a centroid file. "
+                                     "Use CatSim to calculate the pixel positions of "
+                                     "objects in the InstanceCatalog.  Compare this "
+                                     "to the pixel positions as reported in the "
+                                     "centroid file.\n\nexample:\n\n"
+                                     "python compare_centroids.py --cat myInstanceCatalog.txt "
+                                     "--cent myCentroidFile.txt --out_dir my/output/director/")
+
+    parser.add_argument("--cat", type=str, help="path to the InstanceCatalog", default=None)
+    parser.add_argument("--cent", type=str, help="path to the centroid file", default=None)
+    parser.add_argument("--clean", type=bool, help="delete old files, if names conflict",
+                        default=False)
+    parser.add_argument("--out_dir", type=str, help="directory where we will put output files",
+                        default=".")
+
+    args = parser.parse_args()
+
+    if args.cat is None or args.cent is None:
+        raise RuntimeError("Must specify an InstanceCatalog and a centroid file.\n"
+                           "You specified %s and %s" % (args.cat, args.cent))
+
+
+    validator = CentroidValidator(args.cat, args.cent)
+    validator.create_tex_file(args.out_dir, args.clean)
+    print validator.get_scalars()
diff --git a/python/desc/twinkles/__init__.py b/python/desc/twinkles/__init__.py
index 3d82353..7189740 100644
--- a/python/desc/twinkles/__init__.py
+++ b/python/desc/twinkles/__init__.py
@@ -3,6 +3,7 @@
 with warnings.catch_warnings():
     warnings.filterwarnings('ignore', 'Duplicate object type id', UserWarning)
     warnings.filterwarnings('ignore', 'duplicate object identifie', UserWarning)
+    from .astrometryUtils import *
     from .analyseICat import *
     from .calc_snr import *
     from .cleanupspectra import *
@@ -17,6 +18,7 @@
     from .twinkles_io import *
     from .twinkles_sky import *
     from .obsHistIDOrdering import *
+    from .validateCentroidFile import *
     try:
         from .version import *
     except ImportError:
diff --git a/python/desc/twinkles/astrometryUtils.py b/python/desc/twinkles/astrometryUtils.py
new file mode 100644
index 0000000..d688648
--- /dev/null
+++ b/python/desc/twinkles/astrometryUtils.py
@@ -0,0 +1,76 @@
+import numpy as np
+import numbers
+from lsst.sims.utils import cartesianFromSpherical, sphericalFromCartesian
+from lsst.sims.utils import rotationMatrixFromVectors
+from lsst.sims.utils import _observedFromICRS, icrsFromObserved
+
+__all__ = ["_rePrecess", "icrsFromPhoSim"]
+
+def _rePrecess(ra_list, dec_list, obs):
+    """
+    Undo the 'dePrecession' done to mangle observed RA, Dec into the
+    RA, Dec coordinates expected by PhoSim
+
+    Parameters
+    ----------
+    ra_list a list of dePrecessed RA (in 'observed' frame)  in radians
+
+    dec_list is a list of dePrecessed Dec (in 'observed' frame) in radians.
+
+    *see lsst.sims.utils.AstrometryUtils.observedFromICRS for definition
+    of 'observed' RA, Dec
+
+    obs is an ObservationMetaData characterizing the telescope pointing
+
+    Returns
+    -------
+    RA and Dec in 'observed' frame in radians
+    """
+
+    xyz_bore = cartesianFromSpherical(np.array([obs._pointingRA]),
+                                      np.array([obs._pointingDec]))
+
+    precessedRA, precessedDec = _observedFromICRS(np.array([obs._pointingRA]),
+                                                  np.array([obs._pointingDec]),
+                                                  obs_metadata=obs,
+                                                  epoch=2000.0,
+                                                  includeRefraction=False)
+
+    xyz_precessed = cartesianFromSpherical(precessedRA, precessedDec)
+
+    rotMat = rotationMatrixFromVectors(xyz_bore[0], xyz_precessed[0])
+
+    xyz_list = cartesianFromSpherical(ra_list, dec_list)
+
+    if isinstance(ra_list, numbers.Number):
+        xyz_re_precessed = np.dot(rotMat, xyz_list)
+    else:
+        xyz_re_precessed = np.array([np.dot(rotMat, xx) for xx in xyz_list])
+
+    ra_re_precessed, dec_re_precessed = sphericalFromCartesian(xyz_re_precessed)
+
+    return ra_re_precessed, dec_re_precessed
+
+
+def icrsFromPhoSim(ra_list, dec_list, obs):
+    """
+    Parameters
+    ----------
+    ra_list is the PhoSim ra value in degrees
+
+    dec_list is the PhoSim dec value in degrees
+
+    obs is the ObservationMetaData characterizing the pointing
+
+    Returns
+    -------
+    RA, Dec in ICRS in degrees
+    """
+
+    ra_precessed, dec_precessed = _rePrecess(np.radians(ra_list),
+                                             np.radians(dec_list),
+                                             obs)
+
+    return icrsFromObserved(np.degrees(ra_precessed), np.degrees(dec_precessed),
+                            obs_metadata=obs, epoch=2000.0,
+                            includeRefraction=False)
diff --git a/python/desc/twinkles/validateCentroidFile.py b/python/desc/twinkles/validateCentroidFile.py
new file mode 100644
index 0000000..143d2ab
--- /dev/null
+++ b/python/desc/twinkles/validateCentroidFile.py
@@ -0,0 +1,438 @@
+from __future__ import with_statement
+
+import matplotlib
+matplotlib.use('Agg')
+
+import matplotlib.pyplot as plt
+from matplotlib.colors import LogNorm
+
+import numpy as np
+import pandas
+import os
+import warnings
+
+from lsst.sims.utils import ObservationMetaData
+from lsst.sims.coordUtils import pixelCoordsFromRaDec
+from lsst.obs.lsstSim import LsstSimMapper
+from desc.twinkles import icrsFromPhoSim
+
+__all__ = ["getPredictedCentroids", "CentroidValidator"]
+
+def getPredictedCentroids(cat_name):
+    """
+    Read in an InstanceCatalog.  Use CatSim to calculate the expected
+    pixel positions of all of the sources in that InstanceCatalog.
+    Return a pandas dataframe containing each source's id, xpix, and ypix
+    coordinates.
+
+    Parameters
+    ----------
+    cat_name is the path to the InstanceCatalog
+
+    Returns
+    -------
+    a pandas dataframe with columns 'id', 'x', and 'y'
+    """
+
+    target_chip = 'R:2,2 S:1,1'
+
+    camera = LsstSimMapper().camera
+
+    with open(cat_name, 'r') as input_catalog:
+        input_lines = input_catalog.readlines()
+
+    ra = None
+    dec = None
+    mjd = None
+    rotSkyPos = None
+    while ra is None or dec is None or mjd is None or rotSkyPos is None:
+        for line in input_lines:
+            line = line.split()
+            if line[0] == 'rightascension':
+                ra = float(line[1])
+            if line[0] == 'declination':
+                dec = float(line[1])
+            if line[0] == 'mjd':
+                mjd = float(line[1])
+            if line[0] == 'rotskypos':
+                rotSkyPos = float(line[1])
+
+    try:
+        assert ra is not None
+        assert dec is not None
+        assert mjd is not None
+        assert rotSkyPos is not None
+    except:
+        print 'ra: ',ra
+        print 'dec: ',dec
+        print 'mjd: ',mjd
+        print 'rotSkyPos: ',rotSkyPos
+
+    obs = ObservationMetaData(pointingRA=ra,
+                              pointingDec=dec,
+                              mjd=mjd,
+                              rotSkyPos=rotSkyPos)
+
+    id_list = []
+    ra_list = []
+    dec_list = []
+    for line in input_lines:
+        line = line.split()
+        if len(line) > 2:
+            id_list.append(int(line[1]))
+            ra_list.append(float(line[2]))
+            dec_list.append(float(line[3]))
+
+    id_list = np.array(id_list)
+    ra_list = np.array(ra_list)
+    dec_list = np.array(dec_list)
+
+    ra_icrs, dec_icrs = icrsFromPhoSim(ra_list, dec_list, obs)
+
+    x_pix, y_pix = pixelCoordsFromRaDec(ra_icrs, dec_icrs,
+                                        chipName=[target_chip]*len(id_list),
+                                        #chipName=chip_name_list,
+                                        obs_metadata=obs,
+                                        camera=camera)
+
+    return pandas.DataFrame({'id': id_list,
+                             'x': y_pix,
+                             'y': x_pix})  # need to reverse pixel order because
+                                           # DM and PhoSim have different
+                                           # conventions
+
+
+class CentroidValidator(object):
+    """
+    Read in an InstanceCatalog (specifiedy by cat_name)and a centroid file
+    (specified by centroid_name).  Calculate the predicted pixel positions
+    of the sources in the InstanceCatalog using Catsim.  Compare those to
+    the pixel positions reported in the centroid file and compute some summary
+    statistics.  Includes an optional method to create a .tex file
+    displaying the distribution of disagreements between PhoSim and CatSim
+    pixel positions.
+    """
+
+    def __init__(self, cat_name, centroid_name):
+        self.cat_name = cat_name
+        self.centroid_name = centroid_name
+        self.catsim_data = getPredictedCentroids(self.cat_name)
+
+        self._x_center = 2000.0
+        self._y_center = 2036.0  # central pixel on chip
+
+        phosim_dtype = np.dtype([('id', long), ('nphot', int),
+                                 ('x', float), ('y', float)])
+
+        _phosim_data = np.genfromtxt(self.centroid_name, dtype=phosim_dtype,
+                                     skip_header=1)
+
+        self.phosim_data = pandas.DataFrame({'id': _phosim_data['id'],
+                                             'nphot': _phosim_data['nphot'],
+                                             'x': _phosim_data['x'],
+                                             'y':_phosim_data['y']})
+
+        # make sure that any sources whicha appear in the PhoSim centroid file, but
+        # not in the CatSim-predicted centroid file have id==0
+        self.just_phosim = self.phosim_data[np.logical_not(self.phosim_data.id.isin(self.catsim_data.id.values).values)]
+
+        try:
+            assert self.just_phosim.id.max() == 0
+        except:
+            print 'a source with non-zero ID appears in PhoSim centroid file, but not CatSim centroid file'
+            raise
+
+        # find all of the CatSim sources that appeared in PhoSim
+        self.catsim_phosim = self.catsim_data.merge(self.phosim_data, left_on='id', right_on='id',
+                                                    how='left', suffixes=('_catsim', '_phosim'))
+
+        self.catsim_phosim['dx'] = pandas.Series(self.catsim_phosim['x_catsim']-self.catsim_phosim['x_phosim'],
+                                                 index=self.catsim_phosim.index)
+        self.catsim_phosim['dy'] = pandas.Series(self.catsim_phosim['y_catsim']-self.catsim_phosim['y_phosim'],
+                                                 index=self.catsim_phosim.index)
+
+        # select points that actually showed up on the PhoSim image
+        overlap_dexes = np.logical_not(np.logical_or(self.catsim_phosim['x_phosim'].isnull(),
+                                                     self.catsim_phosim['y_phosim'].isnull()))
+
+        overlap = self.catsim_phosim[overlap_dexes]
+
+        self.bright_sources = overlap.query('nphot>0')
+        self.bright_sources = self.bright_sources.sort_values(by='nphot')
+
+        self.just_catsim = self.catsim_data[np.logical_not(self.catsim_data.id.isin(self.phosim_data.id.values).values)]
+        self.min_d_just_catsim = np.sqrt(np.power(self.just_catsim.x-self._x_center,2) +
+                                         np.power(self.just_catsim.y-self._y_center,2)).min()
+        print 'minimum distance of just_catsim: ',self.min_d_just_catsim
+
+        nphot_sum = self.bright_sources.nphot.sum()
+        self.weighted_dx = (self.bright_sources.dx*self.bright_sources.nphot).sum()/nphot_sum
+        self.weighted_dy = (self.bright_sources.dy*self.bright_sources.nphot).sum()/nphot_sum
+        self.mean_dx = self.bright_sources.dx.mean()
+        self.mean_dy = self.bright_sources.dy.mean()
+        self.median_dx = self.bright_sources.dx.median()
+        self.median_dy = self.bright_sources.dy.median()
+
+        print 'weighted dx/dy: ',self.weighted_dx, self.weighted_dy
+        print 'mean dx/dy: ',self.mean_dx, self.mean_dy
+        print 'median dx/dy: ',self.median_dx, self.median_dy
+
+    def get_scalars(self):
+        """
+        Return a dict of the scalar values associated with validation of this InstanceCatalog
+        """
+        return {'mean_dx': (self.mean_dx, 'The mean value of x(CatSim)-x(PhoSim) in pixels'),
+                'mean_dy': (self.mean_dy, 'The mean value of y(CatSim)-y(PhoSim) in pixels'),
+                'weighted_dx': (self.weighted_dx, 'The mean value of x(CatSim)-x(PhoSim) in pixels weighted by number of photons'),
+                'weighted_dy': (self.weighted_dy, 'The mean value of y(CatSim)-y(PhoSim) in pixels weighted by number of photons'),
+                'median_dx': (self.median_dx, 'The median value of x(CatSim)-x(PhoSim) in pixels'),
+                'median_dy': (self.median_dy, 'The median value of y(CatSim)-y(PhoSim) in pixels'),
+                'd_just_catsim': (self.min_d_just_catsim,
+                                  'The minimum distance in pixels from the center of the chip to the sources in the InstanceCatalog '
+                                  'that do not appear in the centroid file')}
+
+    def create_tex_file(self, out_dir, clean=False):
+        """
+        Create a .tex file that can be compiled with pdflatex to display the distribution
+        of disagreements between PhoSim and CatSim pixel positions.
+
+        out_dir is the directory where we should write the figures
+
+        clean is a boolean which controls whether or not to overwrite existing files
+        (True means overwrite).  If clean is False and the files you want to create
+        already exist, this method will append an integer to the file names so that
+        they are unique.  A message will be printed to stdout, informing you of this
+        name mangling.
+        """
+
+        if not os.path.exists(out_dir):
+            os.mkdir(out_dir)
+
+        self.scatter_fig_name = os.path.join(out_dir, 'dx_dy_scatter.png')
+        self.radial_fig_name = os.path.join(out_dir, 'dx_dy_radius.png')
+        self.displacement_fig_name = os.path.join(out_dir, 'max_displacement.png')
+        self.tex_name = os.path.join(out_dir, self.centroid_name.replace('.txt','_validation.tex'))
+
+        if os.path.exists(self.tex_name) and clean:
+            os.unlink(self.tex_name)
+
+        if os.path.exists(self.scatter_fig_name) and clean:
+            os.unlink(self.scatter_fig_name)
+
+        if os.path.exists(self.displacement_fig_name) and clean:
+            os.unlink(self.displacement_fig_name)
+
+        if os.path.exists(self.radial_fig_name) and clean:
+            os.unlink(self.radial_fig_name)
+
+        if os.path.exists(self.scatter_fig_name) or os.path.exists(self.displacement_fig_name) \
+        or os.path.exists(self.tex_name) or os.path.exists(self.radial_fig_name):
+
+            scatter_root = self.scatter_fig_name.replace('.png', '')
+            displacement_root = self.displacement_fig_name.replace('.png', '')
+            tex_root = self.tex_name.replace('.tex', '')
+            radial_root = self.radial_fig_name.replace('.png', '')
+
+            ix = 1
+            while os.path.exists(self.scatter_fig_name) or os.path.exists(self.displacement_fig_name) \
+            or os.path.exists(self.tex_name):
+
+                self.scatter_fig_name = scatter_root+'_%d.png' % ix
+                self.displacement_fig_name = displacement_root+'_%d.png' % ix
+                self.tex_name = tex_root+'_%d.tex' % ix
+                self.radial_fig_name = radial_root+'_%d.png' % ix
+                ix += 1
+
+            warnings.warn('Needed to rename figures to %s, %s, %s, %s to avoid overwriting older files' %
+                          (self.scatter_fig_name, self.displacement_fig_name,
+                           self.radial_fig_name, self.tex_name))
+
+        self._create_figures(out_dir, clean)
+
+        if os.path.exists(self.tex_name):
+            raise RuntimeError('%s already exists; not going to overwrite it\n' % tex_name)
+
+        with open(self.tex_name, 'w') as tex_file:
+            self.tex_opening_boilerplate(tex_file)
+            self.tex_figure(tex_file, self.scatter_fig_name.split('/')[-1],
+                            'The displacement between CatSim and PhoSim in pixel coordinates for each source. '
+                            'Color bar indicates number of photons in the source.')
+            self.tex_figure(tex_file, self.displacement_fig_name.split('/')[-1],
+                            'Maximum displacement in pixel coordinates as a function of number of '
+                            'photons in the source')
+            self.tex_figure(tex_file, self.radial_fig_name.split('/')[-1],
+                            'Displacement between PhoSim and CatSim as a function of distance from '
+                            'chip center (in PhoSim)')
+            self.tex_scalar(tex_file, self.weighted_dx, 'Weighted (by nphoton) mean displacement in x')
+            self.tex_scalar(tex_file, self.weighted_dy, 'Weighted (by nphoton) mean displacement in y')
+            self.tex_scalar(tex_file, self.mean_dx, 'Mean displacement in x')
+            self.tex_scalar(tex_file, self.mean_dy, 'Mean displacement in y')
+            self.tex_scalar(tex_file, self.median_dx, 'Median displacement in x')
+            self.tex_scalar(tex_file, self.median_dy, 'Median displacement in y')
+            self.tex_scalar(tex_file, self.min_d_just_catsim,
+                            'Minimum distance (in pixels) from the center of the chip '
+                            'of objects that appear in CatSim but not PhoSim')
+            self.tex_closing_boilerplate(tex_file)
+
+        print 'created file: %s\ncompile it with pdflatex' % self.tex_name
+
+    def _create_figures(self, out_dir, clean=False):
+        """
+        Create the figures expected by create_tex_file()
+
+        out_dir is the directory where we should write the figures
+
+        clean is a boolean which controls whether or not to overwrite existing files
+        (True means overwrite)
+        """
+        if os.path.exists(self.scatter_fig_name):
+            raise RuntimeError('%s already exists; not going to overwrite it' % scatter_fig_name)
+
+        _dx = self.bright_sources.dx.max()-self.bright_sources.dx.min()
+        _dy = self.bright_sources.dy.max()-self.bright_sources.dy.min()
+
+        plt.figsize=(30,30)
+        for i_fig, limit in enumerate([([-20+self.weighted_dx, 20+self.weighted_dx],[-20+self.weighted_dy, 20+self.weighted_dy]),
+                                       ([-50+self.weighted_dx, 50+self.weighted_dx],[-50+self.weighted_dy, 50+self.weighted_dy]),
+                                       ([-200+self.weighted_dx,200+self.weighted_dx],[-200+self.weighted_dy,200+self.weighted_dy]),
+                                       ([self.bright_sources.dx.min()-0.01*_dx, self.bright_sources.dx.max()+0.01*_dx],
+                                        [self.bright_sources.dy.min()-0.01*_dy, self.bright_sources.dy.max()+0.01*_dy])]):
+            plt.subplot(2,2,i_fig+1)
+            plt.scatter(self.bright_sources.dx, self.bright_sources.dy, c=self.bright_sources.nphot,
+                        s=10,edgecolor='',cmap=plt.cm.gist_ncar,norm=LogNorm())
+
+
+            plt.plot([limit[0][0], limit[0][1]], [0.0, 0.0], color='k', linestyle='--')
+            plt.plot([0.0, 0.0], [limit[1][0], limit[1][1]], color='k', linestyle='--')
+
+            plt.xlim(limit[0])
+            plt.ylim(limit[1])
+            if limit[0][1]<0.0 or limit[0][0]>0.0:
+                dx_tick = (limit[0][1]-limit[0][0])/5.0
+            else:
+                if 0.0 - limit[0][0] > 0.5*limit[0][1]:
+                    dx_tick = (0.0-limit[0][0])/2.0
+                else:
+                    dx_tick = limit[0][1]/2.0
+                    limit[0][0] = -dx_tick
+
+            xticks = np.arange(limit[0][0],limit[0][1],dx_tick)
+            xtick_labels = ['%.2e' % vv if ix%2==0 else '' for ix, vv in enumerate(xticks)]
+
+            if limit[1][1]<0.0 or limit[1][0]>0.0:
+                dy_tick = (limit[1][1]-limit[1][0])/5.0
+            else:
+                if 0.0 - limit[1][0] > 0.5*limit[1][1]:
+                    dy_tick = (0.0-limit[1][0])/2.0
+                else:
+                    dy_tick = limit[1][1]/2.0
+                    limit[1][0] = -dy_tick
+
+            yticks = np.arange(limit[1][0],limit[1][1],dy_tick)
+            ytick_labels = ['%.2e' % vv if ix%2==0 else '' for ix, vv in enumerate(yticks)]
+
+            plt.xticks(xticks, xtick_labels, fontsize=10)
+            plt.yticks(yticks, ytick_labels, fontsize=10)
+            plt.minorticks_on()
+            plt.tick_params(axis='both', length=10)
+
+            if i_fig==0:
+                cb = plt.colorbar()
+                cb.set_label('photons in source')
+                plt.xlabel('dx (pixels)')
+                plt.ylabel('dy (pixels)')
+
+        plt.tight_layout()
+        plt.savefig(self.scatter_fig_name)
+        plt.close()
+
+        if os.path.exists(self.radial_fig_name):
+            raise RuntimeError("%s exists; not willing to overwrite it" % self.radial_fig_name)
+
+        d_rr = np.sqrt(np.power(self.bright_sources.dx, 2) + np.power(self.bright_sources.dy, 2))
+        center_dist = np.sqrt(np.power(self.bright_sources.x_phosim-self._x_center,2) +
+                              np.power(self.bright_sources.y_phosim-self._y_center,2))
+
+        plt.figsize = (30, 30)
+        fig = plt.figure()
+        ax = plt.gca()
+        ax.scatter(center_dist, d_rr)
+        plt.xlabel('distance from center in PhoSim (pixels)')
+        plt.ylabel('displacement between PhoSim and CatSim (pixels)')
+        ax.set_yscale('log')
+        ax.set_xscale('log')
+        plt.savefig(self.radial_fig_name)
+
+        nphot_unique = np.unique(self.bright_sources.nphot)
+        sorted_dex = np.argsort(-1.0*nphot_unique)
+        nphot_unique = nphot_unique[sorted_dex]
+
+        dx_of_nphot = np.array([np.abs(self.bright_sources.query('nphot>=%e' % nn).dx).max() for nn in nphot_unique[1:]])
+        dy_of_nphot = np.array([np.abs(self.bright_sources.query('nphot>=%e' % nn).dy).max() for nn in nphot_unique[1:]])
+
+        if os.path.exists(self.displacement_fig_name):
+            raise RuntimeError('%s already exists; will not overwrite it' % self.displacement_fig_name)
+
+        plt.figsize = (30, 30)
+        plt.subplot(2,2,1)
+        plt.semilogx(nphot_unique[1:], dx_of_nphot, label='dx')
+        plt.semilogx(nphot_unique[1:], dy_of_nphot, label='dy')
+        plt.xlabel('minimum number of photons')
+        plt.ylabel('max pixel displacement')
+        plt.subplot(2,2,2)
+        plt.semilogx(nphot_unique[1:], dx_of_nphot, label='dx')
+        plt.semilogx(nphot_unique[1:], dy_of_nphot, label='dy')
+        plt.ylim((0,200))
+        plt.subplot(2,2,3)
+        plt.semilogx(nphot_unique[1:], dx_of_nphot, label='dx')
+        plt.semilogx(nphot_unique[1:], dy_of_nphot, label='dy')
+        plt.ylim((0,100))
+        plt.subplot(2,2,4)
+        plt.semilogx(nphot_unique[1:], dx_of_nphot, label='dx')
+        plt.semilogx(nphot_unique[1:], dy_of_nphot, label='dy')
+        plt.ylim((0,50))
+        plt.legend(loc=0)
+        plt.tight_layout()
+        plt.savefig(self.displacement_fig_name)
+
+
+    def tex_opening_boilerplate(self, file_handle):
+        """
+        Write the boilerplate needed at the opening of the tex file.
+        file_handle is a file object pointing to the tex file.
+        """
+        file_handle.write('\documentclass[preprint]{aastex}\n')
+        file_handle.write('\\topmargin 0.0cm\n\\textheight 8.5in\n')
+        file_handle.write('\\begin{document}\n\n')
+
+
+    def tex_closing_boilerplate(self, file_handle):
+        """
+        Write the boilerplate needed at the end of the tex file.
+        file_handle is a file object pointing to the tex file.
+        """
+        file_handle.write('\n\end{document}\n')
+
+
+    def tex_figure(self, file_handle, fig_name, caption):
+        """
+        Add a figure to the tex file.
+        file_handle is a file object pointing to the tex file.
+        fig_name is the name of the file containing the figure.
+        caption is the caption that should appear with the figure (a str).
+        """
+        file_handle.write('\n\\begin{figure}\n')
+        file_handle.write('\includegraphics[scale=0.9]{%s}\n' % fig_name)
+        file_handle.write('\caption{\n')
+        file_handle.write('%s\n' % caption)
+        file_handle.write('}\n')
+        file_handle.write('\end{figure}\n\n')
+
+    def tex_scalar(self, file_handle, value, caption):
+        """
+        Add a scalar to the tex file.
+        file_handle is a file object pointing to the texfile.
+        value is the numerical value of the scalar.
+        captin is a string that will appear before the scalar value in the tex file.
+        """
+        file_handle.write('\n\n%s: $%.12e$\n\n' % (caption, value))
diff --git a/tests/scratchSpace/.gitignore b/tests/scratchSpace/.gitignore
new file mode 100644
index 0000000..ccf1cce
--- /dev/null
+++ b/tests/scratchSpace/.gitignore
@@ -0,0 +1,2 @@
+*.db
+*.txt
diff --git a/tests/test_astrometry_utils.py b/tests/test_astrometry_utils.py
new file mode 100644
index 0000000..d0281ee
--- /dev/null
+++ b/tests/test_astrometry_utils.py
@@ -0,0 +1,161 @@
+from __future__ import with_statement
+import unittest
+import numpy as np
+import os
+
+from lsst.sims.utils import ObservationMetaData, raDecFromAltAz
+from lsst.sims.utils import cartesianFromSpherical
+from lsst.sims.catUtils.mixins import PhoSimAstrometryBase
+from lsst.utils import getPackageDir
+from lsst.sims.catalogs.db import fileDBObject
+from lsst.sims.catalogs.definitions import InstanceCatalog
+from lsst.sims.catalogs.decorators import compound
+from lsst.sims.catUtils.mixins import PhoSimAstrometryStars
+from lsst.sims.utils import _icrsFromObserved
+from desc.twinkles import _rePrecess, icrsFromPhoSim
+from desc.twinkles import StarCacheDBObj
+
+class PrecessionTestCase(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        """
+        This method creates a dummy database that looks like our
+        cache of stars from fatboy, but only contains astrometrically
+        relevant quantities, ignoring photometrically relevant quantities.
+        """
+        cls.scratch_dir = os.path.join(getPackageDir('twinkles'),
+                                       'tests', 'scratchSpace')
+
+        rng = np.random.RandomState(81822)
+        mjd = 59781.2
+        ra, dec = raDecFromAltAz(60.0, 112.0, ObservationMetaData(mjd=mjd))
+
+        rot = 12.0
+        cls.cat_obs = ObservationMetaData(pointingRA=ra,
+                                          pointingDec=dec,
+                                          rotSkyPos=rot,
+                                          mjd=mjd,
+                                          boundType='circle',
+                                          boundLength=1.75)
+
+        n_obj = 2000
+        rr = rng.random_sample(n_obj)*1.75
+        theta = rng.random_sample(n_obj)*2.0*np.pi
+        ra_list = ra + rr*np.cos(theta)
+        dec_list = dec + rr*np.sin(theta)
+        mura_list = (rng.random_sample(n_obj)-0.5)*100.0
+        mudec_list = (rng.random_sample(n_obj)-0.5)*100.0
+        px_list = rng.random_sample(n_obj)*50.0
+        vrad_list = (rng.random_sample(n_obj)-0.5)*600.0
+
+        cat_txt_name = os.path.join(cls.scratch_dir,
+                                    'test_stars_cat_source.txt')
+        if os.path.exists(cat_txt_name):
+            os.unlink(cat_txt_name)
+
+        with open(cat_txt_name, 'w') as output_file:
+            for ix, (ra, dec, mura, mudec, px, vrad) in \
+            enumerate(zip(ra_list, dec_list, mura_list,
+                          mudec_list, px_list, vrad_list)):
+
+                output_file.write('%d;%.16g;%.16g;%.16g;%.16g;%.16g;%.16g\n'
+                                  % (ix, ra, dec, mura, mudec, px, vrad))
+
+        dtype = np.dtype([('simobjid', int), ('ra', float), ('decl', float),
+                          ('mura', float), ('mudecl', float),
+                          ('parallax', float), ('vrad', float)])
+
+        cls.db_name = os.path.join(cls.scratch_dir, 'test_stars_cata_sqlite.db')
+
+        if os.path.exists(cls.db_name):
+            os.unlink(cls.db_name)
+
+        fileDBObject(cat_txt_name, database=cls.db_name, dtype=dtype,
+                     delimiter=';', driver='sqlite', idColKey='simobjid',
+                     runtable='star_cache_table')
+
+        if os.path.exists(cat_txt_name):
+            os.unlink(cat_txt_name)
+
+    @classmethod
+    def tearDownClass(cls):
+        if os.path.exists(cls.db_name):
+            os.unlink(cls.db_name)
+
+    def test_re_precess(self):
+        """
+        Test that _rePrecess undoes the _dePrecess method from
+        PhoSimAstrometryBase
+        """
+
+        ra = 25.0
+        dec = -113.0
+        rotSkyPos = 32.1
+        mjd = 59571.2
+        obs = ObservationMetaData(pointingRA=ra, pointingDec=dec,
+                                  rotSkyPos=rotSkyPos, mjd=mjd)
+
+        rng = np.random.RandomState(33)
+        n_obj = 1000
+        rr = rng.random_sample(n_obj)*3.0
+        theta = rng.random_sample(n_obj)*2.0*np.pi
+        ra_list = np.radians(ra + rr*np.cos(theta))
+        dec_list = np.radians(dec + rr*np.sin(theta))
+        ast = PhoSimAstrometryBase()
+        ra_precessed, dec_precessed = ast._dePrecess(ra_list, dec_list, obs)
+        ra_test, dec_test = _rePrecess(ra_precessed, dec_precessed, obs)
+
+        xyz0 = cartesianFromSpherical(ra_list, dec_list)
+        xyz1 = cartesianFromSpherical(ra_test, dec_test)
+        for control, test in zip(xyz0, xyz1):
+            dd = np.sqrt(np.power(xyz0-xyz1,2).sum())
+            self.assertLess(dd, 1.0e-10)
+
+        # try it one at a time
+        for ix, (ra_in, dec_in) in enumerate(zip(ra_precessed, dec_precessed)):
+            ra_test, dec_test = _rePrecess(ra_in, dec_in, obs)
+            xyz0 = cartesianFromSpherical(ra_list[ix], dec_list[ix])
+            xyz1 = cartesianFromSpherical(ra_test, dec_test)
+            xyz2 = cartesianFromSpherical(ra_precessed[ix], dec_precessed[ix])
+            self.assertEqual(xyz0.shape, (3,))
+            self.assertEqual(xyz1.shape, (3,))
+            dd = np.sqrt(np.power(xyz0-xyz1, 2).sum())
+            self.assertLess(dd, 1.0e-10)
+            dd2 = np.sqrt(np.power(xyz0-xyz2, 2).sum())
+            self.assertGreater(dd2,1.0e-6)
+
+    def test_end_to_end(self):
+        """
+        Test that we can transform from PhoSim RA, Dec back to
+        ICRS RA, Dec
+        """
+
+        class testRaDecCat(PhoSimAstrometryStars, InstanceCatalog):
+            column_outputs = ['raPhoSim', 'decPhoSim', 'raControl', 'decControl']
+
+            transformations = {'raPhoSim': np.degrees, 'decPhoSim': np.degrees,
+                               'raControl': np.degrees, 'decControl': np.degrees}
+
+            @compound('raControl', 'decControl')
+            def get_controlCoords(self):
+                return _icrsFromObserved(self.column_by_name('raObserved'),
+                                         self.column_by_name('decObserved'),
+                                         obs_metadata=self.obs_metadata,
+                                         epoch=2000.0)
+
+        StarCacheDBObj.database = self.db_name
+        db = StarCacheDBObj()
+        cat = testRaDecCat(db, obs_metadata=self.cat_obs)
+        ct = 0
+        for star in cat.iter_catalog():
+            ra_test, dec_test = icrsFromPhoSim(star[0], star[1], self.cat_obs)
+            xyz0 = cartesianFromSpherical(np.radians(star[2]), np.radians(star[3]))
+            xyz1 = cartesianFromSpherical(np.radians(ra_test), np.radians(dec_test))
+            dd = np.sqrt(np.power(xyz0-xyz1,2).sum())
+            self.assertLess(dd, 1.0e-10)
+            ct += 1
+        self.assertGreater(ct, 0)
+
+if __name__ == "__main__":
+    unittest.main()