diff --git a/flystar/align.py b/flystar/align.py
index db37954..eabe887 100755
--- a/flystar/align.py
+++ b/flystar/align.py
@@ -1,10 +1,7 @@
 import numpy as np
-from flystar import match
-from flystar import transforms
-from flystar import plots
+from flystar import match, transforms, plots, motion_model
 from flystar.starlists import StarList
 from flystar.startables import StarTable
-from flystar import motion_model
 from astropy.table import Table, Column, vstack
 import datetime
 import copy
@@ -28,7 +25,7 @@ def __init__(self, list_of_starlists, ref_index=0, iters=2,
                  default_motion_model='Fixed',
                  motion_model_dict = {},
                  use_scipy=True, 
-                 absolute_sigma=False, 
+                 absolute_sigma=True, 
                  save_path=None, 
                  verbose=True):
         """
@@ -207,10 +204,10 @@ def = None. If not None, then this should contain an array or list of transform
         self.verbose = verbose
 
         # For backwards compatibility.
-        if self.verbose is True:
-            self.verbose = 9
-        if self.verbose is False:
-            self.verbose = 0
+        # if self.verbose is True:
+        #     self.verbose = 9
+        # if self.verbose is False:
+        #     self.verbose = 0
             
         self.N_lists = len(self.star_lists)
 
@@ -255,7 +252,7 @@ def fix_iterable_conditions(self):
 
         if self.mag_lim is None:
             self.mag_lim = np.repeat([[None, None]], len(self.star_lists), axis=0)
-        elif (len(self.mag_lim) == 2):
+        elif (len(self.mag_lim) == 2) and (np.ndim(self.mag_lim) == 1):
             self.mag_lim = np.repeat([self.mag_lim], len(self.star_lists), axis=0)
         assert len(self.mag_lim) == len(self.star_lists)
 
@@ -420,10 +417,24 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar
             else:
                 star_list_T.transform_xy(trans) 
 
-            # Match stars between the transformed, trimmed lists.
-            idx1, idx2, dr, dm = match.match(star_list_T['x'], star_list_T['y'], star_list_T['m'],
-                                             ref_list['x'], ref_list['y'], ref_list['m'],
-                                             dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose)
+            if 'use_in_trans' in ref_list.colnames:
+                # Only use stars specified by "use_in_trans" column.
+                use_in_trans = ref_list['use_in_trans']
+                # Match stars between the transformed, trimmed lists.
+                idx1, idx2, dr, dm = match.match(
+                    star_list_T['x'], star_list_T['y'], star_list_T['m'],
+                    ref_list['x'][use_in_trans], ref_list['y'][use_in_trans], ref_list['m'][use_in_trans],
+                    dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose
+                )
+                # Restore idx2 to the full reference list indices
+                idx2 = np.where(use_in_trans)[0][idx2]
+            else:
+                idx1, idx2, dr, dm = match.match(
+                    star_list_T['x'], star_list_T['y'], star_list_T['m'],
+                    ref_list['x'], ref_list['y'], ref_list['m'],
+                    dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose
+                )
+
             if self.verbose > 1:
                 print( '  Match 1: Found ', len(idx1), ' matches out of ', len(star_list_T),
                        '. If match count is low, check dr_tol, dm_tol.' )
@@ -438,16 +449,6 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar
                 idx1 = idx1[keepers]
                 idx2 = idx2[keepers]
 
-            # Only use stars specified by "use_in_trans" column.
-            if 'use_in_trans' in ref_list.colnames:
-                keepers = np.where(ref_list[idx2]['use_in_trans'] == True)[0]
-                
-                if self.verbose > 1:
-                    print( '  Rejected ', len(idx1) - len(keepers), ' with use_in_trans=False.' )
-                    
-                idx1 = idx1[keepers]
-                idx2 = idx2[keepers]
-
             # Determine weights in the fit.
             weight = self.get_weights_for_lists(ref_list[idx2], star_list_T[idx1])
 
@@ -506,11 +507,11 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar
                                      dy=(star_t['y'] - star_r['y']) * 1e3,
                                      dm=(star_t['m'] - star_r['m']),
                                      xo=star_s['x'], yo=star_s['y'], mo=star_s['m']))
-                    
+
             idx_lis, idx_ref, dr, dm = match.match(star_list_T['x'], star_list_T['y'], star_list_T['m'],
                                                    ref_list['x'], ref_list['y'], ref_list['m'],
                                                    dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose)
-            
+
             if self.verbose > 1:
                 print( '  Match 2: After trans, found ', len(idx_lis), ' matches out of ', len(star_list_T),
                        '. If match count is low, check dr_tol, dm_tol.' )
@@ -858,7 +859,7 @@ def update_ref_table_aggregates(self, keep_orig=None, n_boot=0):
             fit_star_idxs = [idx for idx in range(len(self.ref_table)) if idx not in keep_orig]
         else:
             fit_star_idxs = None
-        #pdb.set_trace()
+
         # Figure out whether motion fits are necessary
         all_fixed = np.all(self.ref_table['motion_model_input']=='Fixed')
         if all_fixed:
@@ -1240,7 +1241,6 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
                                                                    m=starlist_boot['m'], mref=ref_boot['m'],
                                                                    weights=weight, mag_trans=self.mag_trans)
                 #print(jj)
-                #pdb.set_trace()
 
                 # Apply transformation to *all* orig positions in this epoch. Need to make a new
                 # FLYSTAR starlist object with the original positions for this. We don't
@@ -1378,7 +1378,6 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
                 
                 col[idx_good] = data_dict[ff]
                 self.ref_table.add_column(col)
-        #pdb.set_trace()
         
         print('===============================')
         print('Done with bootstrap')
@@ -1416,7 +1415,7 @@ def __init__(self, ref_list, list_of_starlists, iters=2,
                  default_motion_model='Fixed',
                  motion_model_dict={},
                  use_scipy=True,
-                 absolute_sigma=False,
+                 absolute_sigma=True,
                  save_path=None,
                  verbose=True):
 
@@ -2403,8 +2402,7 @@ def write_transform(transform, starlist, reference, N_trans, deltaMag=0, restric
         Xcoeff = transform.px.parameters
         Ycoeff = transform.py.parameters
     else:
-        print(( '{0} not yet supported!'.format(transType)))
-        return
+        raise TypeError(( '{0} not yet supported!'.format(trans_name)))
         
     # Write output
     _out = open(outFile, 'w')
@@ -2632,7 +2630,7 @@ def position_transform_from_object(x, y, xe, ye, transform):
         order = transform.order
     else:
         txt = 'Transform not yet supported by position_transform_from_object'
-        raise StandardError(txt)
+        raise TypeError(txt)
         
     # How the transformation is applied depends on the type of transform.
     # This can be determined by the length of Xcoeff, Ycoeff
@@ -2731,7 +2729,7 @@ def velocity_transform_from_object(x0, y0, x0e, y0e, vx, vy, vxe, vye, transform
         order = transform.order
     else:
         txt = 'Transform not yet supported by velocity_transform_from_object'
-        raise StandardError(txt)
+        raise TypeError(txt)
         
     # How the transformation is applied depends on the type of transform.
     # This can be determined by the length of Xcoeff, Ycoeff
@@ -2894,8 +2892,8 @@ def trans_initial_guess(ref_list, star_list, trans_args, motion_model_dict, mode
     if mode == 'name':
         # First trim the two lists down to only those that don't contain
         # the "ignore_contains" string.
-        idx_r = np.flatnonzero(np.char.find(ref_list['name'], ignore_contains) == -1)
-        idx_s = np.flatnonzero(np.char.find(star_list['name'], ignore_contains) == -1)
+        idx_r = np.flatnonzero(np.char.find(ref_list['name'].astype(str), ignore_contains) == -1)
+        idx_s = np.flatnonzero(np.char.find(star_list['name'].astype(str), ignore_contains) == -1)
 
         # Match the star names
         name_matches, ndx_r, ndx_s = np.intersect1d(ref_list['name'][idx_r],
@@ -2954,7 +2952,11 @@ def trans_initial_guess(ref_list, star_list, trans_args, motion_model_dict, mode
         trans.mag_offset = 0
         
     if verbose > 1:
-        print('init guess: ', trans.px.parameters, trans.py.parameters)
+        # print('init guess: ', trans.px.parameters, trans.py.parameters)
+        print('Initial guess:')
+        print(f'{trans.px.parameters=}')
+        print(f'{trans.py.parameters=}')
+        print(f'{trans.mag_offset=}')
 
     warnings.filterwarnings('default', category=AstropyUserWarning)
         
@@ -3018,53 +3020,6 @@ def copy_and_rename_for_ref(star_list):
 
     return ref_list
 
-def outlier_rejection_indices(star_list, ref_list, outlier_tol, verbose=True):
-    """
-    Determine the outliers based on the residual positions between two different
-    starlists and some threshold (in sigma). Return the indices of the stars 
-    to keep (that shouldn't be rejected as outliers). 
-
-    Note that we assume that the star_list and ref_list are already transformed and
-    matched. 
-
-    Parameters
-    ----------
-    star_list : StarList
-        starlist with 'x', 'y'
-
-    ref_list : StarList
-        starlist with 'x0', 'y0'
-
-    outlier_tol : float
-        Number of sigma inside which we keep stars and outside of which we 
-        reject stars as outliers. 
-
-    Optional Parameters
-    --------------------
-    verbose : boolean
-
-    Returns
-    ----------
-    keepers : nd.array
-        The indicies of the stars to keep. 
-    """
-    # Optionally propogate the reference positions forward in time.
-    xref, yref = get_pos_in_time(star_list['t'][0], ref_list)
-    
-    # Residuals
-    x_resid_on_old_trans = star_list['x'] - xref
-    y_resid_on_old_trans = star_list['y'] - yref
-    resid_on_old_trans = np.hypot(x_resid_on_old_trans, y_resid_on_old_trans)
-
-    threshold = outlier_tol * resid_on_old_trans.std()
-    keepers = np.where(resid_on_old_trans < threshold)[0]
-
-    if verbose > 0:
-        msg = '  Outlier Rejection: Keeping {0:d} of {1:d}'
-        print(msg.format(len(keepers), len(resid_on_old_trans)))
-        
-    return keepers
-
 def setup_trans_info(trans_input, trans_args, N_lists, iters):
     """ Setup transformation info into a usable format.
 
@@ -3193,3 +3148,210 @@ def logger(logfile, message, verbose = 9):
         print(message)
     logfile.write(message + '\n')
     return
+
+
+def generic_match(sl1, sl2, init_mode='triangle',
+                  model=transforms.PolyTransform, order_dr=(1, 1.0),
+                  dr_final=1.0,
+                  xy_match=(None, None, None, None, None, None, None, None),
+                  m_match=(None, None, None, None), sigma_match=None,
+                  n_bright=100, verbose=True, **kwargs):
+    """
+    Finds the transformation between two starlists using the first one
+    as reference frame. Different matching methods can be used. If no
+    transformation is found, it returns an error message.
+
+
+    Parameters
+    sl1 : StarList
+        starlist used for reference frame
+    sl2 : StarList
+        starlist transformed
+    init_mode : str
+        Initial matching method.
+        If 'triangle', uses the blind triangle method.
+        If 'match_name', uses match by name
+        If 'load', uses the transformation from a loaded file
+    model : str
+        Transformation model to be used with the 'triangle' initial mode
+    poly_order : int
+        Order of the transformation model
+    order_dr : int, float [n, 2]
+        Combinations of polinomial order (first column) and search radius
+        (second column) to refine the transformation. Rows are executed in
+        orders
+    dr_final: float
+        Search radius used for the final matching
+    n_bright : int
+        Number of bright stars used in the initial blind triangles matching
+    xy_match : array
+        Area of the images to remove in the matching [reference catalog min x,
+        reference catalog max x, reference catalog min y, reference catalog max y,
+        transformed catalog min x, transformed catalog max x,
+        transformed catalog min y, transformed catalog max y]. Use None for values not used.
+    m_match : array
+        Magnitude limits of matching stars used to find transformations
+        [reference catalog min mag, reference catalog max mag, transformed
+        catalog min mag, transformed catalog max mag]. Use None for values not
+        used
+    sigma_match : array
+        Number of Deltap movement sigmas [0] used for sigma-cutting matched
+        stars for a number of times [1]. Use None for no sigma-cut. The last
+        polynomial order and search radius in 'order_dr' are used
+    transf_file : str
+        File name and path of the transformation file used with the 'load'
+        init_mode
+    verbose : bool, optional
+        Prints on screen information on the matching
+    
+    Returns
+    -------
+    transf : Transform2D
+        Transformation of the second starlist respect to the first
+    st : StarTable
+        Startable of the two matched catalogs
+
+    """
+    from flystar import starlists, startables
+    #  Check the input StarLists and transform them into astropy Tables
+    if not isinstance(sl1, starlists.StarList):
+        raise TypeError("The first catalog has to be a StarList")
+    if not isinstance(sl2, starlists.StarList):
+        raise TypeError("The second catalog has to be a StarList")
+    
+    #  Find the initial transformation
+    if init_mode == 'triangle': #  Blind triangles method
+        
+        #  Prepare the reduced starlists for matching
+        sl1_cut = copy.deepcopy(sl1)
+        sl2_cut = copy.deepcopy(sl2)
+        sl1_cut.restrict_by_value(x_min=xy_match[0], x_max=xy_match[1],
+                                  y_min=xy_match[2], y_max=xy_match[3])
+        sl2_cut.restrict_by_value(x_min=xy_match[4], x_max=xy_match[5],
+                                  y_min=xy_match[6], y_max=xy_match[7])
+        sl1_cut.restrict_by_value(m_min=m_match[0], m_max=m_match[1])
+        sl2_cut.restrict_by_value(m_min=m_match[2], m_max=m_match[3])
+        
+        # Find the transformation
+        # TODO: test 'initial_align' with StarList input
+        transf = initial_align(sl1_cut, sl2_cut, briteN=n_bright,
+                                     transformModel=model, order=order_dr[0]) #order_dr[i_loop][0] ?
+        
+    elif init_mode == 'match_name': #  Name match
+        sl1_idx_init, sl2_idx_init, _ = starlists.restrict_by_name(sl1, sl2)
+        transf = model(sl2['x'][sl2_idx_init], sl2['y'][sl2_idx_init],
+                       sl1['x'][sl1_idx_init], sl1['y'][sl1_idx_init],
+                       order=int(order_dr[0][0]))
+        
+    elif init_mode == 'load': #  Load a transformation file
+        transf = transforms.Transform2D.from_file(kwargs['transf_file'])
+        
+    else: #  None of the above
+        raise TypeError("Unrecognized initial matching method")
+
+    # Restrict the matching catalogs
+    sl1_match = copy.deepcopy(sl1)
+    sl2_match = copy.deepcopy(sl2)
+    sl1_match.restrict_by_value(m_min=m_match[0], m_max=m_match[1])
+    sl2_match.restrict_by_value(m_min=m_match[2], m_max=m_match[3])
+    
+    #  Refine the transformation
+    if sigma_match:
+        order_dr_len = len(order_dr)
+        
+        for i_loop in range(sigma_match[1]):
+            order_dr = np.vstack((np.array(order_dr), np.array(order_dr[-1])))
+    
+    for i_loop in range(len(order_dr)):
+        
+        #  Transform and match the catalog to the reference frame
+#        sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf,
+#                                                     dr_tol=order_dr[i_loop][1],
+#                                                     verbose=verbose)
+
+        sl2_idx, sl1_idx = transform_and_match(sl2_match, sl1_match, transf,
+                                                     dr_tol=order_dr[1],
+                                                     verbose=verbose)
+
+        #  Transform the catalog to the reference frame
+        sl2_transf_match = transform_from_object(sl2_match, transf)
+        
+        # Sigma-rejection
+        if sigma_match and (i_loop >= order_dr_len):
+            resid = np.sqrt((sl1_match['x'][sl1_idx] -
+                            sl2_transf_match['x'][sl2_idx])**2 +
+                            (sl1_match['y'][sl1_idx] -
+                            sl2_transf_match['y'][sl2_idx])**2)
+            sl1_idx = sl1_idx[resid <= (sigma_match[0] * np.std(resid))]
+            sl2_idx = sl2_idx[resid <= (sigma_match[0] * np.std(resid))]
+        
+        # Test section to observe the matching catalogs before refining the transformation
+        """
+        from matplotlib import pyplot
+        
+        _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10,10))
+        axarr.scatter(sl1_match['x'][sl1_idx], sl1_match['y'][sl1_idx])
+        xlim = axarr.get_xlim()
+        ylim = axarr.get_ylim()
+        
+        _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10, 10))
+        axarr.scatter(sl2_transf_match['x'][sl2_idx], sl2_transf_match['y'][sl2_idx])
+        axarr.set_xlim(xlim)
+        axarr.set_ylim(ylim)
+        """
+        
+        #  Find a better transformation
+        transf, _ = find_transform(sl2_match[sl2_idx],
+                                         sl2_transf_match[sl2_idx],
+                                         sl1_match[sl1_idx], transModel=model,
+                                         order=order_dr[0], verbose=verbose)
+#                                         order=int(order_dr[i_loop][0]), verbose=verbose)
+        
+        # This section was used for testing transformations with normalized
+        # coordinates. Only several catalogs had reduced residuals when using
+        # high order polynomials (>3), some of them became unstable
+        """sl1_match_norm = sl1_match[sl1_idx]
+        sl2_match_norm = sl2_match[sl2_idx]
+        sl2_transf_match_norm = sl2_transf_match[sl2_idx]
+        mm = max(max(sl1_match_norm['x']), max(sl1_match_norm['y']),
+                 max(sl2_transf_match_norm['x']), max(sl2_transf_match_norm['y']))
+        sl1_match_norm['x'] = sl1_match_norm['x'] / mm
+        sl1_match_norm['y'] = sl1_match_norm['y'] / mm
+        sl2_match_norm['x'] = sl2_match_norm['x'] / mm
+        sl2_match_norm['y'] = sl2_match_norm['y'] / mm
+        sl2_transf_match_norm['x'] = sl2_transf_match_norm['x'] / mm
+        sl2_transf_match_norm['y'] = sl2_transf_match_norm['y'] / mm
+        transf, _ = align.find_transform(sl2_match_norm, sl2_transf_match_norm,
+                                         sl1_match_norm, transModel=model,
+                                         order=poly_order, verbose=verbose)
+        c_exp = np.zeros(len(transf.px._parameters))
+        
+        for i_c in range(len(transf.px._parameters)):
+            c_exp[i_c] = int(transf.px._param_names[i_c][1:].split('_')[0]) +\
+                         int(transf.px._param_names[i_c][1:].split('_')[1])
+        
+        c_corr = mm ** (1 - c_exp)
+        transf.px._parameters = transf.px._parameters * c_corr
+        transf.py._parameters = transf.py._parameters * c_corr"""
+    
+    # Do the final transformation and matching using
+    sl2_idx, sl1_idx = transform_and_match(sl2, sl1, transf, dr_tol=dr_final,
+                                                 verbose=verbose)
+    #  StarTable output
+    sl2_transf = transform_from_object(sl2, transf)
+    unames = np.array(range(len(sl1_idx)))
+    st = startables.StarTable(name=unames,
+         x=np.column_stack((np.array(sl1['x'][sl1_idx]), np.array(sl2_transf['x'][sl2_idx]))),
+         y=np.column_stack((np.array(sl1['y'][sl1_idx]), np.array(sl2_transf['y'][sl2_idx]))),
+         m=np.column_stack((np.array(sl1['m'][sl1_idx]), np.array(sl2_transf['m'][sl2_idx]))),
+         ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))))
+#         ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))),
+#         list_times=[sl1.meta['list_time'], sl2.meta['list_time']],
+#         list_names=[sl1.meta['list_name'], sl2.meta['list_name']])
+    
+    for col in sl1.colnames:
+        if col in sl2.colnames:
+            if col not in ['name', 'x', 'y', 'm']:
+                st.add_column(Column(np.column_stack((np.array(sl1[col][sl1_idx]),np.array(sl2_transf[col][sl2_idx]))), name=col))
+    
+    return transf, st
diff --git a/flystar/analysis.py b/flystar/analysis.py
index f300723..1a2ea82 100644
--- a/flystar/analysis.py
+++ b/flystar/analysis.py
@@ -1,19 +1,12 @@
 import numpy as np
 import pylab as plt
-from flystar import starlists
-from flystar import startables
-from flystar import align
-from flystar import match
-from flystar import transforms
+from flystar import starlists, match
 from astropy import table
 from astropy.table import Table, Column
 from astropy.coordinates import SkyCoord
 from astropy import units as u
-from astropy.wcs import WCS
 from astroquery.gaia import Gaia
-from astroquery.mast import Observations, Catalogs
 import pdb, copy
-import math
 from scipy.stats import f
 
 ##################################################
@@ -49,7 +42,7 @@ def query_gaia(ra, dec, search_radius=30.0, table_name='gaiadr3'):
     search_radius *= u.arcsec
 
     Gaia.ROW_LIMIT = 50000
-    gaia_job = Gaia.cone_search_async(target_coords, search_radius, table_name = table_name + '.gaia_source')
+    gaia_job = Gaia.cone_search_async(target_coords, radius=search_radius, table_name=table_name + '.gaia_source')
     gaia = gaia_job.get_results()
 
     #Change new 'SOURCE_ID' column header back to lowercase 'source_id' so all subsequent functions still work:
@@ -475,85 +468,6 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False):
 # Old codes.
 ##################################################
 
-def calc_chi2(ref_mat, starlist_mat, transform, errs='both'):
-    """
-    calculate the chi2 and reduced chi2 of the position 
-    between two matched starlists.
-    Input:
-    ref_mat: astropy table
-        Reference starlist only containing matched stars that were used in the
-        transformation. Standard column headers are assumed.
-        
-    starlist_mat: astropy table
-        Transformed starlist only containing the matched stars used in
-        the transformation. Standard column headers are assumed.
-
-    transform: transformation object
-        Transformation object of final transform. Used in chi-square
-        determination
-
-    errs: string; 'both', 'reference', or 'starlist'
-        If both, add starlist errors in quadrature with reference errors.
-
-        If reference, only consider reference errors. This should be used if the starlist
-        does not have valid errors
-
-        If starlist, only consider starlist errors. This should be used if the reference
-        does not have valid errors
-
-    Output:
-    chi_sq: float
-        chi2 = sum (diff_x**2 / xerr**2 + diff_y**2 /yerr**2)
-    chi_sq_red: float
-        reduced chi2 = chi2/ degree of freedom
-    deg_freedom: int
-        degree of freedom
-
-    """
-    diff_x = ref_mat['x'] - starlist_mat['x']
-    diff_y = ref_mat['y'] - starlist_mat['y']
-
-    # Set errors as per user input
-    if errs == 'both':
-        xerr = np.hypot(ref_mat['xe'], starlist_mat['xe'])
-        yerr = np.hypot(ref_mat['ye'], starlist_mat['ye'])
-    elif errs == 'reference':
-        xerr = ref_mat['xe']
-        yerr = ref_mat['ye']
-    elif errs == 'starlist':
-        xerr = starlist_mat['xe']
-        yerr = starlist_mat['ye']
-          
-
-    # For both X and Y, calculate chi-square. Combine arrays to get combined
-    # chi-square
-    chi_sq_x = diff_x**2. / xerr**2.
-    chi_sq_y = diff_y**2. / yerr**2.
-
-    chi_sq = np.append(chi_sq_x, chi_sq_y)
-    
-    # Calculate degrees of freedom in transformation
-    num_mod_params = calc_nparam(transform)
-    deg_freedom = len(chi_sq) - num_mod_params
-    
-    # Calculate reduced chi-square
-    chi_sq = np.sum(chi_sq)
-    chi_sq_red = chi_sq / deg_freedom
-
-    return chi_sq, chi_sq_red, deg_freedom
-
-
-def calc_nparam(transformation):
-    """
-    calculate the degree of freedom for a transformation
-    """
-    # Read transformation: Extract X, Y coefficients from transform
-    if transformation.__class__.__name__ == 'four_paramNW':
-        nparam = 4
-    elif transformation.__class__.__name__ == 'PolyTransform':
-        order = transformation.order
-        nparam = (order+1) * (order+2) 
-    return nparam
 
 def calc_F(red_chi2_1, red_chi2_2, v1, v2):
     """
diff --git a/flystar/examples.py b/flystar/examples.py
index 8059562..0e0a042 100644
--- a/flystar/examples.py
+++ b/flystar/examples.py
@@ -1,8 +1,4 @@
-from flystar import transforms
-from flystar import match
-from flystar import align
-from flystar import starlists
-from flystar import plots
+from . import transforms, match, align, starlists, plots
 import numpy as np
 import copy
 import pdb
diff --git a/flystar/match.py b/flystar/match.py
index bba108a..55b707e 100644
--- a/flystar/match.py
+++ b/flystar/match.py
@@ -1,5 +1,4 @@
 import numpy as np
-from flystar import starlists, transforms, startables, align
 from collections import Counter
 from scipy.spatial import cKDTree as KDT
 from astropy.table import Column, Table
@@ -462,210 +461,3 @@ def add_votes(votes, match1, match2):
     votes.flat[unique_idx] += deltas
     
     return
-
-
-def generic_match(sl1, sl2, init_mode='triangle',
-                  model=transforms.PolyTransform, order_dr=(1, 1.0),
-                  dr_final=1.0,
-                  xy_match=(None, None, None, None, None, None, None, None),
-                  m_match=(None, None, None, None), sigma_match=None,
-                  n_bright=100, verbose=True, **kwargs):
-    """
-    Finds the transformation between two starlists using the first one
-    as reference frame. Different matching methods can be used. If no
-    transformation is found, it returns an error message.
-
-
-    Parameters
-    sl1 : StarList
-        starlist used for reference frame
-    sl2 : StarList
-        starlist transformed
-    init_mode : str
-        Initial matching method.
-        If 'triangle', uses the blind triangle method.
-        If 'match_name', uses match by name
-        If 'load', uses the transformation from a loaded file
-    model : str
-        Transformation model to be used with the 'triangle' initial mode
-    poly_order : int
-        Order of the transformation model
-    order_dr : int, float [n, 2]
-        Combinations of polinomial order (first column) and search radius
-        (second column) to refine the transformation. Rows are executed in
-        orders
-    dr_final: float
-        Search radius used for the final matching
-    n_bright : int
-        Number of bright stars used in the initial blind triangles matching
-    xy_match : array
-        Area of the images to remove in the matching [reference catalog min x,
-        reference catalog max x, reference catalog min y, reference catalog max y,
-        transformed catalog min x, transformed catalog max x,
-        transformed catalog min y, transformed catalog max y]. Use None for values not used.
-    m_match : array
-        Magnitude limits of matching stars used to find transformations
-        [reference catalog min mag, reference catalog max mag, transformed
-        catalog min mag, transformed catalog max mag]. Use None for values not
-        used
-    sigma_match : array
-        Number of Deltap movement sigmas [0] used for sigma-cutting matched
-        stars for a number of times [1]. Use None for no sigma-cut. The last
-        polynomial order and search radius in 'order_dr' are used
-    transf_file : str
-        File name and path of the transformation file used with the 'load'
-        init_mode
-    verbose : bool, optional
-        Prints on screen information on the matching
-    
-    Returns
-    -------
-    transf : Transform2D
-        Transformation of the second starlist respect to the first
-    st : StarTable
-        Startable of the two matched catalogs
-
-    """
-    
-    #  Check the input StarLists and transform them into astropy Tables
-    if not isinstance(sl1, starlists.StarList):
-        raise TypeError("The first catalog has to be a StarList")
-    if not isinstance(sl2, starlists.StarList):
-        raise TypeError("The second catalog has to be a StarList")
-    
-    #  Find the initial transformation
-    if init_mode == 'triangle': #  Blind triangles method
-        
-        #  Prepare the reduced starlists for matching
-        sl1_cut = copy.deepcopy(sl1)
-        sl2_cut = copy.deepcopy(sl2)
-        sl1_cut.restrict_by_value(x_min=xy_match[0], x_max=xy_match[1],
-                                  y_min=xy_match[2], y_max=xy_match[3])
-        sl2_cut.restrict_by_value(x_min=xy_match[4], x_max=xy_match[5],
-                                  y_min=xy_match[6], y_max=xy_match[7])
-        sl1_cut.restrict_by_value(m_min=m_match[0], m_max=m_match[1])
-        sl2_cut.restrict_by_value(m_min=m_match[2], m_max=m_match[3])
-        
-        # Find the transformation
-        # TODO: test 'initial_align' with StarList input
-        transf = align.initial_align(sl1_cut, sl2_cut, briteN=n_bright,
-                                     transformModel=model, order=order_dr[0]) #order_dr[i_loop][0] ?
-        
-    elif init_mode == 'match_name': #  Name match
-        sl1_idx_init, sl2_idx_init, _ = starlists.restrict_by_name(sl1, sl2)
-        transf = model(sl2['x'][sl2_idx_init], sl2['y'][sl2_idx_init],
-                       sl1['x'][sl1_idx_init], sl1['y'][sl1_idx_init],
-                       order=int(order_dr[0][0]))
-        
-    elif init_mode == 'load': #  Load a transformation file
-        transf = transforms.Transform2D.from_file(kwargs['transf_file'])
-        
-    else: #  None of the above
-        raise TypeError("Unrecognized initial matching method")
-
-    # Restrict the matching catalogs
-    sl1_match = copy.deepcopy(sl1)
-    sl2_match = copy.deepcopy(sl2)
-    sl1_match.restrict_by_value(m_min=m_match[0], m_max=m_match[1])
-    sl2_match.restrict_by_value(m_min=m_match[2], m_max=m_match[3])
-    
-    #  Refine the transformation
-    if sigma_match:
-        order_dr_len = len(order_dr)
-        
-        for i_loop in range(sigma_match[1]):
-            order_dr = np.vstack((np.array(order_dr), np.array(order_dr[-1])))
-    
-    for i_loop in range(len(order_dr)):
-        
-        #  Transform and match the catalog to the reference frame
-#        sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf,
-#                                                     dr_tol=order_dr[i_loop][1],
-#                                                     verbose=verbose)
-
-        sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf,
-                                                     dr_tol=order_dr[1],
-                                                     verbose=verbose)
-
-        #  Transform the catalog to the reference frame
-        sl2_transf_match = align.transform_from_object(sl2_match, transf)
-        
-        # Sigma-rejection
-        if sigma_match and (i_loop >= order_dr_len):
-            resid = np.sqrt((sl1_match['x'][sl1_idx] -
-                            sl2_transf_match['x'][sl2_idx])**2 +
-                            (sl1_match['y'][sl1_idx] -
-                            sl2_transf_match['y'][sl2_idx])**2)
-            sl1_idx = sl1_idx[resid <= (sigma_match[0] * np.std(resid))]
-            sl2_idx = sl2_idx[resid <= (sigma_match[0] * np.std(resid))]
-        
-        # Test section to observe the matching catalogs before refining the transformation
-        """
-        from matplotlib import pyplot
-        
-        _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10,10))
-        axarr.scatter(sl1_match['x'][sl1_idx], sl1_match['y'][sl1_idx])
-        xlim = axarr.get_xlim()
-        ylim = axarr.get_ylim()
-        
-        _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10, 10))
-        axarr.scatter(sl2_transf_match['x'][sl2_idx], sl2_transf_match['y'][sl2_idx])
-        axarr.set_xlim(xlim)
-        axarr.set_ylim(ylim)
-        """
-        
-        #  Find a better transformation
-        transf, _ = align.find_transform(sl2_match[sl2_idx],
-                                         sl2_transf_match[sl2_idx],
-                                         sl1_match[sl1_idx], transModel=model,
-                                         order=order_dr[0], verbose=verbose)
-#                                         order=int(order_dr[i_loop][0]), verbose=verbose)
-        
-        # This section was used for testing transformations with normalized
-        # coordinates. Only several catalogs had reduced residuals when using
-        # high order polynomials (>3), some of them became unstable
-        """sl1_match_norm = sl1_match[sl1_idx]
-        sl2_match_norm = sl2_match[sl2_idx]
-        sl2_transf_match_norm = sl2_transf_match[sl2_idx]
-        mm = max(max(sl1_match_norm['x']), max(sl1_match_norm['y']),
-                 max(sl2_transf_match_norm['x']), max(sl2_transf_match_norm['y']))
-        sl1_match_norm['x'] = sl1_match_norm['x'] / mm
-        sl1_match_norm['y'] = sl1_match_norm['y'] / mm
-        sl2_match_norm['x'] = sl2_match_norm['x'] / mm
-        sl2_match_norm['y'] = sl2_match_norm['y'] / mm
-        sl2_transf_match_norm['x'] = sl2_transf_match_norm['x'] / mm
-        sl2_transf_match_norm['y'] = sl2_transf_match_norm['y'] / mm
-        transf, _ = align.find_transform(sl2_match_norm, sl2_transf_match_norm,
-                                         sl1_match_norm, transModel=model,
-                                         order=poly_order, verbose=verbose)
-        c_exp = np.zeros(len(transf.px._parameters))
-        
-        for i_c in range(len(transf.px._parameters)):
-            c_exp[i_c] = int(transf.px._param_names[i_c][1:].split('_')[0]) +\
-                         int(transf.px._param_names[i_c][1:].split('_')[1])
-        
-        c_corr = mm ** (1 - c_exp)
-        transf.px._parameters = transf.px._parameters * c_corr
-        transf.py._parameters = transf.py._parameters * c_corr"""
-    
-    # Do the final transformation and matching using
-    sl2_idx, sl1_idx = align.transform_and_match(sl2, sl1, transf, dr_tol=dr_final,
-                                                 verbose=verbose)
-    #  StarTable output
-    sl2_transf = align.transform_from_object(sl2, transf)
-    unames = np.array(range(len(sl1_idx)))
-    st = startables.StarTable(name=unames,
-         x=np.column_stack((np.array(sl1['x'][sl1_idx]), np.array(sl2_transf['x'][sl2_idx]))),
-         y=np.column_stack((np.array(sl1['y'][sl1_idx]), np.array(sl2_transf['y'][sl2_idx]))),
-         m=np.column_stack((np.array(sl1['m'][sl1_idx]), np.array(sl2_transf['m'][sl2_idx]))),
-         ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))))
-#         ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))),
-#         list_times=[sl1.meta['list_time'], sl2.meta['list_time']],
-#         list_names=[sl1.meta['list_name'], sl2.meta['list_name']])
-    
-    for col in sl1.colnames:
-        if col in sl2.colnames:
-            if col not in ['name', 'x', 'y', 'm']:
-                st.add_column(Column(np.column_stack((np.array(sl1[col][sl1_idx]),np.array(sl2_transf[col][sl2_idx]))), name=col))
-    
-    return transf, st
diff --git a/flystar/motion_model.py b/flystar/motion_model.py
index 0b86d07..c4964d1 100644
--- a/flystar/motion_model.py
+++ b/flystar/motion_model.py
@@ -61,7 +61,16 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var',
         """
         # Run a single fit (used both for overall fit + bootstrap iterations)
         pass
-        
+    
+    def calc_sigma(self, xe, ye, weighting='var'):
+        if weighting=='std':
+            return np.sqrt(np.abs(xe)), np.sqrt(np.abs(ye))
+        elif weighting=='var':
+            return np.abs(xe), np.abs(ye)
+        else:
+            warnings.warn("Invalid weighting, using default weighting scheme var.", UserWarning)
+            return np.abs(xe), np.abs(ye)
+            
     def get_weights(self, xe, ye, weighting='var'):
         """
         Get the weights for each data point for fitting. Options are 'var' (default)
@@ -74,18 +83,7 @@ def get_weights(self, xe, ye, weighting='var'):
         else:
             warnings.warn("Invalid weighting, using default weighting scheme var.", UserWarning)
             return 1./xe**2, 1./ye**2
-            
-    def scale_errors(self, errs, weighting='var'):
-        """
-        Rescale the fit result errors as needed, according to the weighting scheme used.
-        """
-        if weighting=='std':
-            return np.array(errs)**2
-        elif weighting=='var':
-            return errs
-        else:
-            warnings.warn("Invalid weighting, using default weighting scheme var.", UserWarning)
-            return errs
+
 
     def fit_motion_model(self, t, x, y, xe, ye, t0, bootstrap=0, weighting='var',
                             use_scipy=True, absolute_sigma=True):
@@ -174,11 +172,16 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None,
             x0,y0,x0e,y0e = x[0],y[0],xe[0],ye[0]
     
         else:
-            x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting)
+            sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting)
+            x_wt, y_wt = 1. / sigma_x**2, 1. / sigma_y**2
+            x_wt /= np.sum(x_wt)
+            y_wt /= np.sum(y_wt)
             x0 = np.average(x, weights=x_wt)
-            x0e = np.sqrt(np.average((x-x0)**2,weights=x_wt))
+            # x0e = np.sqrt(np.average((x-x0)**2,weights=x_wt))
+            x0e = np.sum(x_wt**2 * xe**2)**0.5 # Error propagation
             y0 = np.average(y, weights=y_wt)
-            y0e = np.sqrt(np.average((y-y0)**2,weights=y_wt))
+            # y0e = np.sqrt(np.average((y-y0)**2,weights=y_wt))
+            y0e = np.sum(y_wt**2 * ye**2)**0.5 # Error propagation
         
         params = [x0, y0]
         param_errors = [x0e, y0e]
@@ -226,67 +229,67 @@ def get_batch_pos_at_time(self, t, x0=[],vx=[], y0=[],vy=[], t0=[],
     def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None,
                             use_scipy=True, absolute_sigma=True):
         dt = t-t0
-        x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting)
+        sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting)
+        x_wt, y_wt = 1. / sigma_x**2, 1. / sigma_y**2
         if params_guess is None:
             params_guess = [x.mean(),0.0,y.mean(),0.0]
 
-        # Handle 2-data point case
-        if len(np.unique(dt))==2:
-            if len(x)>2: # Catch case where bootstrap sends only 2 unique epochs
-                _,idx=np.unique(dt, return_index=True)
-                dt = dt[idx]
-                x = x[idx]
-                y = y[idx]
-                xe = xe[idx]
-                ye = ye[idx]
-            dx = np.diff(x)[0]
-            dy = np.diff(y)[0]
-            dt_diff = np.diff(dt)[0]
-            vx = dx / dt_diff
-            vy = dy / dt_diff
-            # TODO: still not sure about the error handling here
-            x0 = x[0] - dt[0]*vx # np.average(x, weights=x_wt) #
-            y0 = y[0] - dt[0]*vy # np.average(y, weights=y_wt) #
-            x0e = np.abs(dx) / 2**0.5 # np.sqrt(np.sum(xe**2)/2) #
-            y0e = np.abs(dy) / 2**0.5 # np.sqrt(np.sum(ye**2)/2) #
-            vxe = 0.0 #np.abs(vx) * np.sqrt(np.sum(xe**2/x**2))
-            vye = 0.0 #np.abs(vy) * np.sqrt(np.sum(ye**2/y**2))
+        # # Handle 2-data point case
+        # if len(np.unique(dt))==2:
+        #     if len(x)>2: # Catch case where bootstrap sends only 2 unique epochs
+        #         _,idx=np.unique(dt, return_index=True)
+        #         dt = dt[idx]
+        #         x = x[idx]
+        #         y = y[idx]
+        #         xe = xe[idx]
+        #         ye = ye[idx]
+        #     dx = np.diff(x)[0]
+        #     dy = np.diff(y)[0]
+        #     dt_diff = np.diff(dt)[0]
+        #     vx = dx / dt_diff
+        #     vy = dy / dt_diff
+        #     # TODO: still not sure about the error handling here
+        #     x0 = x[0] - dt[0]*vx # np.average(x, weights=x_wt) #
+        #     y0 = y[0] - dt[0]*vy # np.average(y, weights=y_wt) #
+        #     x0e = np.abs(dx) / 2**0.5 # np.sqrt(np.sum(xe**2)/2) #
+        #     y0e = np.abs(dy) / 2**0.5 # np.sqrt(np.sum(ye**2)/2) #
+        #     vxe = 0.0 #np.abs(vx) * np.sqrt(np.sum(xe**2/x**2))
+        #     vye = 0.0 #np.abs(vy) * np.sqrt(np.sum(ye**2/y**2))
             
+        # else:
+        if use_scipy:
+            def linear(t, c0, c1):
+                return c0 + c1*t
+            x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=sigma_x, absolute_sigma=absolute_sigma)
+            y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=sigma_y, absolute_sigma=absolute_sigma)
+            x0, vx = x_opt
+            y0, vy = y_opt
+            x0e, vxe = np.sqrt(x_cov.diagonal())
+            y0e, vye = np.sqrt(y_cov.diagonal())
+
         else:
-            if use_scipy:
-                def linear(t, c0, c1):
-                    return c0 + c1*t
-                x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=1/np.sqrt(x_wt), absolute_sigma=absolute_sigma)
-                y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=1/np.sqrt(y_wt), absolute_sigma=absolute_sigma)
-                x0, vx = x_opt
-                y0, vy = y_opt
-                x0e, vxe = np.sqrt(x_cov.diagonal())
-                y0e, vye = np.sqrt(y_cov.diagonal())
-                x0e, vxe, y0e, vye = self.scale_errors([x0e, vxe, y0e, vye], weighting=weighting)
-            else:
-                # Use  https://en.wikipedia.org/wiki/Weighted_least_squares#Solution scheme
-                x = np.array(x)
-                y = np.array(y)
-                dt = np.array(dt)
-                X_mat_t = np.vander(dt, 2)
-                # x calculation
-                W_mat_x = np.diag(x_wt)
-                XTWX_mat_x = X_mat_t.T @ W_mat_x @ X_mat_t
-                pcov_x = np.linalg.inv(XTWX_mat_x)  # Covariance Matrix
-                popt_x = pcov_x @ X_mat_t.T @ W_mat_x @ x   # Linear Solution
-                perr_x = np.sqrt(np.diag(pcov_x))   # Uncertainty of Linear Solution
-                # y calculation
-                W_mat_y = np.diag(y_wt)
-                XTWX_mat_y = X_mat_t.T @ W_mat_y @ X_mat_t
-                pcov_y = np.linalg.inv(XTWX_mat_y)  # Covariance Matrix
-                popt_y = pcov_y @ X_mat_t.T @ W_mat_y @ y   # Linear Solution
-                perr_y = np.sqrt(np.diag(pcov_y))   # Uncertainty of Linear Solution
-                # prepare values to return
-                x0, vx = popt_x[1], popt_x[0]
-                y0, vy = popt_y[1], popt_y[0]
-                x0e, vxe = perr_x[1], perr_x[0]
-                y0e, vye = perr_y[1], perr_y[0]
-                x0e, vxe, y0e, vye = self.scale_errors([x0e, vxe, y0e, vye], weighting=weighting)
+            # Use  https://en.wikipedia.org/wiki/Weighted_least_squares#Solution scheme
+            x = np.array(x)
+            y = np.array(y)
+            dt = np.array(dt)
+            X_mat_t = np.vander(dt, 2)
+            # x calculation
+            W_mat_x = np.diag(x_wt)
+            XTWX_mat_x = X_mat_t.T @ W_mat_x @ X_mat_t
+            pcov_x = np.linalg.inv(XTWX_mat_x)  # Covariance Matrix
+            popt_x = pcov_x @ X_mat_t.T @ W_mat_x @ x   # Linear Solution
+            perr_x = np.sqrt(np.diag(pcov_x))   # Uncertainty of Linear Solution
+            # y calculation
+            W_mat_y = np.diag(y_wt)
+            XTWX_mat_y = X_mat_t.T @ W_mat_y @ X_mat_t
+            pcov_y = np.linalg.inv(XTWX_mat_y)  # Covariance Matrix
+            popt_y = pcov_y @ X_mat_t.T @ W_mat_y @ y   # Linear Solution
+            perr_y = np.sqrt(np.diag(pcov_y))   # Uncertainty of Linear Solution
+            # prepare values to return
+            x0, vx = popt_x[1], popt_x[0]
+            y0, vy = popt_y[1], popt_y[0]
+            x0e, vxe = perr_x[1], perr_x[0]
+            y0e, vye = perr_y[1], perr_y[0]
         
         params = [x0, vx, y0, vy]
         param_errors = [x0e, vxe, y0e, vye]
@@ -339,15 +342,14 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None,
         if not use_scipy:
             Warning("Acceleration model has no non-scipy fitter option. Running with scipy.")
         dt = t-t0
-        x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting)
         if params_guess is None:
             params_guess = [x.mean(),0.0,0.0,y.mean(),0.0,0.0]
             
         def accel(t, c0,c1,c2):
             return c0 + c1*t + 0.5*c2*t**2
-            
-        x_opt, x_cov = curve_fit(accel, dt, x, p0=np.array(params_guess[:3]), sigma=1/x_wt**0.5, absolute_sigma=True)
-        y_opt, y_cov = curve_fit(accel, dt, y, p0=np.array(params_guess[3:]), sigma=1/y_wt**0.5, absolute_sigma=True)
+        sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting)
+        x_opt, x_cov = curve_fit(accel, dt, x, p0=np.array(params_guess[:3]), sigma=sigma_x, absolute_sigma=True)
+        y_opt, y_cov = curve_fit(accel, dt, y, p0=np.array(params_guess[3:]), sigma=sigma_y, absolute_sigma=True)
         x0 = x_opt[0]
         y0 = y_opt[0]
         vx0 = x_opt[1]
@@ -357,7 +359,6 @@ def accel(t, c0,c1,c2):
         
         x0e, vx0e, axe = np.sqrt(x_cov.diagonal())
         y0e, vy0e, aye = np.sqrt(y_cov.diagonal())
-        x0e, vx0e, axe, y0e, vy0e, aye = self.scale_errors([x0e, vx0e, axe, y0e, vy0e, aye], weighting=weighting)
 
         params = [x0, vx0, ax, y0, vy0, ay]
         param_errors = [x0e, vx0e, axe, y0e, vy0e, aye]
@@ -372,7 +373,7 @@ class Parallax(MotionModel):
     Optional PA is counterclockwise offset of the image y-axis from North.
     Optional obs parameter describes observer location, default is 'earth'.
     """
-    n_pts_req = 4
+    n_pts_req = 3
     n_params=3
     fitter_param_names = ['x0', 'vx', 'y0', 'vy', 'pi']
     fixed_param_names = ['t0']
@@ -451,7 +452,6 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None,
             Warning("Parallax model has no non-scipy fitter option. Running with scipy.")
         t_mjd = Time(t, format='decimalyear', scale='utc').mjd
         pvec = self.get_parallax_vector(t_mjd)
-        x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting)
         def fit_func(use_t, x0,vx, y0,vy, pi):
             x_res = x0 + vx*(use_t-t0) + pi*pvec[0]
             y_res = y0 + vy*(use_t-t0) + pi*pvec[1]
@@ -463,10 +463,12 @@ def fit_func(use_t, x0,vx, y0,vy, pi):
             idx_first, idx_last = np.argmin(t), np.argmax(t)
             params_guess = [x.mean(),(x[idx_last]-x[idx_first])/(t[idx_last]-t[idx_first]),
                             y.mean(),(y[idx_last]-y[idx_first])/(t[idx_last]-t[idx_first]), 0.1]
+        sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting)
+        sigma = np.hstack([sigma_x, sigma_y])
         res = curve_fit(fit_func, t, np.hstack([x,y]),
-                        p0=params_guess, sigma = 1.0/np.hstack([x_wt,y_wt]))
+                        p0=params_guess, sigma=sigma, absolute_sigma=absolute_sigma)
         x0,vx,y0,vy,pi = res[0]
-        x0_err,vx_err,y0_err,vy_err,pi_err = self.scale_errors(np.sqrt(np.diag(res[1])), weighting=weighting)
+        x0_err,vx_err,y0_err,vy_err,pi_err = np.sqrt(res[1].diagonal())
 
         params = [x0, vx, y0, vy, pi]
         param_errors = [x0_err, vx_err, y0_err, vy_err, pi_err]
diff --git a/flystar/parallax.py b/flystar/parallax.py
index 4792ec6..586bda8 100755
--- a/flystar/parallax.py
+++ b/flystar/parallax.py
@@ -36,7 +36,7 @@ def parallax_in_direction(RA, Dec, mjd, obsLocation='earth', PA=0):
     #print('parallax_in_direction: len(t) = ', len(mjd))
 
     # Munge inputs into astropy format.
-    times = Time(mjd + 2400000.5, format='jd', scale='tdb')
+    times = Time(mjd, format='mjd', scale='tdb')
     coord = SkyCoord(RA, Dec, unit=(units.deg, units.deg))
 
     direction = coord.cartesian.xyz.value
diff --git a/flystar/plots.py b/flystar/plots.py
index 2d65b2c..77642d7 100755
--- a/flystar/plots.py
+++ b/flystar/plots.py
@@ -1,4 +1,4 @@
-from flystar import analysis, motion_model, startables
+from flystar import motion_model
 import pylab as py
 import pylab as plt
 import numpy as np
@@ -17,6 +17,89 @@
 from astropy.coordinates import SkyCoord
 from astropy import units as u
 
+
+# Moved here from analysis
+def calc_chi2(ref_mat, starlist_mat, transform, errs='both'):
+    """
+    calculate the chi2 and reduced chi2 of the position 
+    between two matched starlists.
+    Input:
+    ref_mat: astropy table
+        Reference starlist only containing matched stars that were used in the
+        transformation. Standard column headers are assumed.
+        
+    starlist_mat: astropy table
+        Transformed starlist only containing the matched stars used in
+        the transformation. Standard column headers are assumed.
+
+    transform: transformation object
+        Transformation object of final transform. Used in chi-square
+        determination
+
+    errs: string; 'both', 'reference', or 'starlist'
+        If both, add starlist errors in quadrature with reference errors.
+
+        If reference, only consider reference errors. This should be used if the starlist
+        does not have valid errors
+
+        If starlist, only consider starlist errors. This should be used if the reference
+        does not have valid errors
+
+    Output:
+    chi_sq: float
+        chi2 = sum (diff_x**2 / xerr**2 + diff_y**2 /yerr**2)
+    chi_sq_red: float
+        reduced chi2 = chi2/ degree of freedom
+    deg_freedom: int
+        degree of freedom
+
+    """
+    diff_x = ref_mat['x'] - starlist_mat['x']
+    diff_y = ref_mat['y'] - starlist_mat['y']
+
+    # Set errors as per user input
+    if errs == 'both':
+        xerr = np.hypot(ref_mat['xe'], starlist_mat['xe'])
+        yerr = np.hypot(ref_mat['ye'], starlist_mat['ye'])
+    elif errs == 'reference':
+        xerr = ref_mat['xe']
+        yerr = ref_mat['ye']
+    elif errs == 'starlist':
+        xerr = starlist_mat['xe']
+        yerr = starlist_mat['ye']
+          
+
+    # For both X and Y, calculate chi-square. Combine arrays to get combined
+    # chi-square
+    chi_sq_x = diff_x**2. / xerr**2.
+    chi_sq_y = diff_y**2. / yerr**2.
+
+    chi_sq = np.append(chi_sq_x, chi_sq_y)
+    
+    # Calculate degrees of freedom in transformation
+    num_mod_params = calc_nparam(transform)
+    deg_freedom = len(chi_sq) - num_mod_params
+    
+    # Calculate reduced chi-square
+    chi_sq = np.sum(chi_sq)
+    chi_sq_red = chi_sq / deg_freedom
+
+    return chi_sq, chi_sq_red, deg_freedom
+
+
+def calc_nparam(transformation):
+    """
+    calculate the degree of freedom for a transformation
+    """
+    # Read transformation: Extract X, Y coefficients from transform
+    if transformation.__class__.__name__ == 'four_paramNW':
+        nparam = 4
+    elif transformation.__class__.__name__ == 'PolyTransform':
+        order = transformation.order
+        nparam = (order+1) * (order+2) 
+    return nparam
+
+
 ####################################################
 # Code for making diagnostic plots for astrometry
 # alignment
@@ -226,15 +309,15 @@ def pos_diff_err_hist(ref_mat, starlist_mat, transform, nbins=25, bin_width=None
     chi_sq_red = np.sum(chi_sq) / deg_freedom
     """
     # Chi-square analysis for all stars, including outliers
-    chi_sq, chi_sq_red, deg_freedom = analysis.calc_chi2(ref_mat, starlist_mat,
+    chi_sq, chi_sq_red, deg_freedom = calc_chi2(ref_mat, starlist_mat,
                                                         transform, errs=errs)
     # Chi-square analysis for only non-outlier stars
-    chi_sq_good, chi_sq_red_good, deg_freedom_good = analysis.calc_chi2(ref_mat[good],
+    chi_sq_good, chi_sq_red_good, deg_freedom_good = calc_chi2(ref_mat[good],
                                                                    starlist_mat[good],
                                                                    transform,
                                                                    errs=errs)
     
-    num_mod_params = analysis.calc_nparam(transform)
+    num_mod_params = calc_nparam(transform)
 
     #-------------------------------------------#
     # Plotting
@@ -262,7 +345,7 @@ def pos_diff_err_hist(ref_mat, starlist_mat, transform, nbins=25, bin_width=None
     py.plot(x, norm.pdf(x,mean,sigma), 'g-', linewidth=2)
     
     # Annotate reduced chi-sqared values in plot: with outliers
-    xstr = '$\chi^2_r$ = {0}'.format(np.round(chi_sq_red, decimals=3))
+    xstr = r'$\chi^2_r$ = {0}'.format(np.round(chi_sq_red, decimals=3))
     py.annotate(xstr, xy=(0.3, 0.77), xycoords='figure fraction', color='black')
     txt = r'$\nu$ = 2*{0} - {1} = {2}'.format(len(diff_x), num_mod_params,
                                                  deg_freedom)
@@ -273,7 +356,7 @@ def pos_diff_err_hist(ref_mat, starlist_mat, transform, nbins=25, bin_width=None
     py.annotate(xstr3, xy=(0.25, 0.80), xycoords='figure fraction', color='black')
     
     # Annotate reduced chi-sqared values in plot: without outliers
-    xstr = '$\chi^2_r$ = {0}'.format(np.round(chi_sq_red_good, decimals=3))
+    xstr = r'$\chi^2_r$ = {0}'.format(np.round(chi_sq_red_good, decimals=3))
     py.annotate(xstr, xy=(0.7, 0.8), xycoords='figure fraction', color='black')
     txt = r'$\nu$ = 2*{0} - {1} = {2}'.format(len(good[0]), num_mod_params,
                                                  deg_freedom_good)
@@ -2226,8 +2309,8 @@ def plot_chi2_dist(tab, Ndetect, motion_model_dict={}, xlim=40, n_bins=50, boot_
     plt.hist(x[idx], bins=chi2_bins, histtype='step', label='X', density=True)
     plt.hist(y[idx], bins=chi2_bins, histtype='step', label='Y', density=True)
     plt.plot(chi2_xaxis, chi2.pdf(chi2_xaxis, Ndof), 'r-', alpha=0.6,
-             label='$\chi^2$ ' + str(round(Ndof,2)) + ' dof')
-    plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(round(Ndof,2)))
+             label=r'$\chi^2$ ' + str(round(Ndof,2)) + ' dof')
+    plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(round(Ndof,2)))
     plt.xlim(0, xlim)
     plt.legend()
 
@@ -2390,8 +2473,8 @@ def plot_chi2_dist_per_filter(tab, Ndetect, motion_model_dict={}, xlim=40, n_bin
     plt.hist(x[idx], bins=chi2_bins, histtype='stepfilled', label='RA', density=True, color='skyblue', alpha=0.8, edgecolor='k')
     plt.hist(y[idx], bins=chi2_bins, histtype='stepfilled', label='DEC', density=True, color='orange', alpha=0.8, edgecolor='k')
     plt.plot(chi2_xaxis, chi2.pdf(chi2_xaxis, Ndof), 'r-', alpha=0.6,
-             label='$\chi^2$ ' + str(Ndof) + ' dof')
-    #plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof))
+             label=r'$\chi^2$ ' + str(Ndof) + ' dof')
+    #plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof))
     plt.title(str(filter)+' (N = '+str(len(chi2_x_list))+')', fontsize=22)
     plt.xlim(0, xlim)
     plt.ylabel(r'PDF', fontsize=28)
@@ -2677,8 +2760,8 @@ def plot_chi2_dist_mag(tab, Ndetect, xlim=40, n_bins=30, boot_err=False):
     plt.clf()
     plt.hist(chi2_m[idx], bins=np.arange(xlim*10), histtype='step', density=True)
     plt.plot(chi2_maxis, chi2.pdf(chi2_maxis, Ndof), 'r-', alpha=0.6, 
-             label='$\chi^2$ ' + str(Ndof) + ' dof')
-    plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof))
+             label=r'$\chi^2$ ' + str(Ndof) + ' dof')
+    plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof))
     plt.xlim(0, xlim)
     plt.legend()
 
@@ -2726,8 +2809,8 @@ def plot_chi2_dist_mag_per_filter(tab, Ndetect, mlim=40, n_bins=30, xlim=40, fil
     plt.clf()
     plt.hist(chi2_m[idx], bins=np.arange(xlim*10), label='mag', histtype='stepfilled', density=True, color='green', alpha=0.7, edgecolor='k')
     plt.plot(chi2_maxis, chi2.pdf(chi2_maxis, Ndof), 'r-', alpha=0.6, 
-             label='$\chi^2$ ' + str(Ndof) + ' dof')
-    #plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof))
+             label=r'$\chi^2$ ' + str(Ndof) + ' dof')
+    #plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof))
     plt.xlim(0, xlim)
     plt.xlabel(r'$\chi^{2}$', fontsize=28)
     plt.ylabel(r'PDF', fontsize=28)
diff --git a/flystar/startables.py b/flystar/startables.py
index d75fca9..f4a04c1 100644
--- a/flystar/startables.py
+++ b/flystar/startables.py
@@ -9,7 +9,7 @@
 import pdb
 import time
 import copy
-from flystar import motion_model
+from . import motion_model
 import pandas as pd
 
 class StarTable(Table):
@@ -463,9 +463,15 @@ def combine_lists(self, col_name_in, weights_col=None, mask_val=None,
                 if all(isinstance(item, int) for item in mask_lists):
                     val_2d.mask[:, mask_lists] = True
                 
+                use_lists = np.array([i for i in np.arange(self[col_name_in].data.shape[1]) if i not in mask_lists])
+
             # Throw a warning if mask_lists is not a list
             if not isinstance(mask_lists, list):
-                raise RuntimeError('mask_lists needs to be a list.')
+                raise RuntimeError(f'mask_lists needs to be a list., not {type(mask_lists)}')
+
+        else:
+            # Use all indices
+            use_lists = np.arange(self[col_name_in].data.shape[1])
 
         # Decide if we are going to have weights (before we
         # do the expensive sigma clipping routine). Note that
@@ -500,15 +506,18 @@ def combine_lists(self, col_name_in, weights_col=None, mask_val=None,
         # the N_lists direction (axis=1).
         if wgt_2d is not None:
             avg = np.ma.average(val_2d_clip, weights=wgt_2d, axis=1)
-            std = np.sqrt(np.ma.average((val_2d_clip.T - avg).T**2, weights=wgt_2d, axis=1))
+            # std = np.sqrt(np.ma.average((val_2d_clip.T - avg).T**2, weights=wgt_2d, axis=1))
+            std = np.ma.sqrt(1. / np.ma.sum(wgt_2d, axis=1)) # Error propagation
         else:
             avg = np.ma.mean(val_2d_clip, axis=1)
-            std = np.ma.std(val_2d_clip, axis=1)
+            # std = np.ma.std(val_2d_clip, axis=1)
+            std = np.ma.std(val_2d_clip, axis=1) / np.sqrt(len(use_lists)) # Error propagation
         # To Do: bring the previous uncertainties of stars that are detected
         # in only one input frame.
-        if (weights_col and weights_col in self.colnames) and (val_2d.shape[1] > 1):
-            mask_for_singles = ((~np.isnan(val_2d_clip)).sum(axis=1)==1)
-            std[mask_for_singles]=np.nanmean(err_2d[mask_for_singles], axis=1)
+        # This can be removed now as error propagation won't result in std=0 anymore.
+        # if (weights_col and weights_col in self.colnames) and (val_2d.shape[1] > 1):
+        #     mask_for_singles = ((~np.isnan(val_2d_clip)).sum(axis=1)==1)
+        #     std[mask_for_singles]=np.nanmean(err_2d[mask_for_singles], axis=1)
 
         # Save off our new AVG and STD into new columns with shape (N_stars).
         col_name_avg = col_name_in + '0'
@@ -607,9 +616,12 @@ def fit_velocities(self, weighting='var', use_scipy=True, absolute_sigma=True, b
         # Define output arrays for the best-fit parameters.
         for col in new_col_list:
             # Clean/remove up old arrays.
-            if col in self.colnames: self.remove_column(col)
-            # Add column #TODO: is this good for filling???
-            self.add_column(Column(data = np.full(N_stars, np.nan, dtype=float), name = col))
+            # if col in self.colnames: self.remove_column(col)
+            # # Add column #TODO: is this good for filling???
+            # self.add_column(Column(data = np.full(N_stars, np.nan, dtype=float), name=col))
+            # Keep existing values
+            if col not in self.colnames:
+                self.add_column(Column(data = np.full(N_stars, np.nan, dtype=float), name=col))
 
         # Add a column to keep track of the number of points used in a fit.
         self['n_fit'] = 0
diff --git a/flystar/tests/compare_branches.py b/flystar/tests/compare_branches.py
new file mode 100644
index 0000000..ed55b39
--- /dev/null
+++ b/flystar/tests/compare_branches.py
@@ -0,0 +1,53 @@
+import os
+import pickle
+import matplotlib.pyplot as plt
+from flystar import align, transforms, motion_model
+from flystar.plots import plot_stars
+
+branch = 'mm_rework'  # 'mm_compare' or 'mm_rework'
+
+test_data_path = f'{os.path.expanduser("~")}/Software/flystar/flystar/tests/test_data'
+
+with open(f'{test_data_path}/my_gaia.pkl', 'rb') as f:
+    my_gaia = pickle.load(f)
+with open(f'{test_data_path}/list_of_starlists.pkl', 'rb') as f:
+    list_of_starlists = pickle.load(f)
+ra_deg, dec_deg = 18.0, -30.0
+my_gaia.remove_column('motion_model_used')
+# my_gaia['motion_model_input'] = 'Fixed'
+if branch == 'mm_compare':
+    msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=1,
+                        dr_tol=[0.2], dm_tol=[5],
+                        outlier_tol=[None], mag_lim=[6, 20],
+                        trans_class=transforms.PolyTransform,
+                        trans_args=[{'order': 1}], 
+                        motion_models=['Fixed', 'Parallax'],
+                        fixed_params_dict = {'ra':ra_deg, 'dec':dec_deg, 'pa':0.0, 'obsLocation':'earth'},
+                        use_ref_new=True,
+                        update_ref_orig=False, 
+                        mag_trans=True,
+                        trans_weights='both,std',
+                        init_guess_mode='name', verbose=3)
+elif branch == 'mm_rework':
+    msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=1,
+                        dr_tol=[0.2], dm_tol=[5],
+                        outlier_tol=[None], mag_lim=[6, 20],
+                        trans_class=transforms.PolyTransform,
+                        trans_args=[{'order': 1}], 
+                        default_motion_model='Parallax',
+                        motion_model_dict = {'Parallax': motion_model.Parallax(RA=ra_deg, Dec=dec_deg, PA=0.0, obsLocation='earth')},
+                        use_ref_new=True,
+                        update_ref_orig=False, 
+                        mag_trans=True,
+                        trans_weights='both,std',
+                        init_guess_mode='name', verbose=3)
+
+msc.fit()
+
+# with open(f'{test_data_path}/ref_table_old.pkl', 'wb') as f:
+#     pickle.dump(msc.ref_table, f)
+
+# for i in range(msc.ref_table['x'].shape[1]):
+#     plt.scatter(msc.ref_table['x'][:, i], msc.ref_table['y'][:, i])
+# plt.show()
+# plot_stars(msc.ref_table, msc.ref_table['name'][:3])
\ No newline at end of file
diff --git a/flystar/tests/test_align.py b/flystar/tests/test_align.py
index 2d6b0dc..4d4cc73 100644
--- a/flystar/tests/test_align.py
+++ b/flystar/tests/test_align.py
@@ -1,21 +1,16 @@
-from flystar import align
-from flystar import starlists
-from flystar import startables
-from flystar import transforms
-from flystar import analysis
-from flystar import motion_model
-from astropy.table import Table
 import numpy as np
 import pylab as plt
-import pdb
-import datetime
-import pytest
+import flystar
+from astropy.table import Table
+from flystar import align, starlists, transforms, analysis, motion_model
+
+test_data_path = f'{flystar.__path__[0]}/tests/test_data'
 
 def test_MosaicSelfRef():
     """
     Cross-match and align 4 starlists using the OO version of mosaic lists.
     """
-    list_files = ['A.lis', 'B.lis', 'C.lis', 'D.lis']
+    list_files = [f'{test_data_path}/{f}' for f in ['A.lis', 'B.lis', 'C.lis', 'D.lis']]
     lists = [starlists.StarList.from_lis_file(lf) for lf in list_files]
 
     ##########
@@ -24,7 +19,7 @@ def test_MosaicSelfRef():
     msc = align.MosaicSelfRef(lists, ref_index=0, iters=2,
                               dr_tol=[3, 3], dm_tol=[1, 1],
                               trans_class=transforms.PolyTransform,
-                              verbose=False,
+                              verbose=2,
                               trans_args={'order': 2})
 
     msc.fit()
@@ -91,7 +86,7 @@ def test_MosaicSelfRef_vel_tconst():
     The 4 lists are all taken at the same time (so 0 velocities should result).
     
     """
-    list_files = ['A.lis', 'B.lis', 'C.lis', 'D.lis']
+    list_files = [f'{test_data_path}/{f}' for f in ['A.lis', 'B.lis', 'C.lis', 'D.lis']]
     lists = [starlists.StarList.from_lis_file(lf) for lf in list_files]
 
     ##########
@@ -149,7 +144,7 @@ def test_MosaicSelfRef_vel():
     Cross-match and align 4 starlists using the OO version of mosaic lists.
     
     """
-    list_files = ['A.lis', 'B.lis', 'C.lis', 'D.lis']
+    list_files = [f'{test_data_path}/{f}' for f in ['A.lis', 'B.lis', 'C.lis', 'D.lis']]
     lists = [starlists.StarList.from_lis_file(lf) for lf in list_files]
 
     # Modify the times so that we get velocities out.
@@ -173,7 +168,7 @@ def test_MosaicSelfRef_vel():
                               dr_tol=[5, 3, 3], dm_tol=[1, 1, 0.5], outlier_tol=None,
                               trans_class=transforms.PolyTransform,
                               trans_args={'order': 2}, default_motion_model='Linear',
-                              verbose=False)
+                              verbose=2)
 
     msc.fit()
     
@@ -215,15 +210,8 @@ def test_MosaicSelfRef_vel():
 def test_MosaicToRef():
     make_fake_starlists_poly1(seed=42)
     
-    ref_file = 'random_ref.fits'
-    list_files = ['random_0.fits',
-                  'random_1.fits',
-                  'random_2.fits',
-                  'random_3.fits',
-                  'random_4.fits',
-                  'random_5.fits',
-                  'random_6.fits',
-                  'random_7.fits']
+    ref_file = f'{test_data_path}/random_ref.fits'
+    list_files = [f'{test_data_path}/random_{i}.fits' for i in range(8)]
 
     ref_list = Table.read(ref_file)
 
@@ -272,15 +260,8 @@ def test_MosaicToRef():
 def test_MosaicToRef_p0_vel():
     make_fake_starlists_poly0_vel(seed=42)
     
-    ref_file = 'random_vel_ref.fits'
-    list_files = ['random_vel_p0_0.fits',
-                  'random_vel_p0_1.fits',
-                  'random_vel_p0_2.fits',
-                  'random_vel_p0_3.fits']
-                  #'random_vel_4.fits',
-                  #'random_vel_5.fits',
-                  #'random_vel_6.fits',
-                  #'random_vel_7.fits']
+    ref_file = f'{test_data_path}/random_vel_ref.fits'
+    list_files = [f'{test_data_path}/random_vel_p0_{i}.fits' for i in range(4)]
 
     ref_list = Table.read(ref_file)
 
@@ -338,15 +319,8 @@ def test_MosaicToRef_p0_vel():
 def test_MosaicToRef_vel():
     make_fake_starlists_poly1_vel(seed=42)
     
-    ref_file = 'random_vel_ref.fits'
-    list_files = ['random_vel_0.fits',
-                  'random_vel_1.fits',
-                  'random_vel_2.fits',
-                  'random_vel_3.fits']
-                  #'random_vel_4.fits',
-                  #'random_vel_5.fits',
-                  #'random_vel_6.fits',
-                  #'random_vel_7.fits']
+    ref_file = f'{test_data_path}/random_vel_ref.fits'
+    list_files = [f'{test_data_path}/random_vel_{i}.fits' for i in range(4)]
 
     ref_list = Table.read(ref_file)
 
@@ -404,15 +378,8 @@ def test_MosaicToRef_vel():
 def test_MosaicToRef_acc():
     make_fake_starlists_poly1_acc(seed=42)
     
-    ref_file = 'random_acc_ref.fits'
-    list_files = ['random_acc_0.fits',
-                  'random_acc_1.fits',
-                  'random_acc_2.fits',
-                  'random_acc_3.fits',
-                  'random_acc_4.fits',
-                  'random_acc_5.fits',
-                  'random_acc_6.fits',
-                  'random_acc_7.fits']
+    ref_file = f'{test_data_path}/random_acc_ref.fits'
+    list_files = [f'{test_data_path}/random_acc_{i}.fits' for i in range(8)]
 
     ref_list = Table.read(ref_file)
 
@@ -500,7 +467,7 @@ def make_fake_starlists_shifts():
 
     # Save original positions as reference (1st) list.
     fmt = '{0:10s}  {1:5.2f} 2015.0 {2:9.4f}  {3:9.4f} 0 0 0 0\n'
-    _out = open('random_0.lis', 'w')
+    _out = open(f'{test_data_path}/random_0.lis', 'w')
     for ii in range(N_stars):
         _out.write(fmt.format(name[ii], m[ii], x[ii], y[ii]))
     _out.close()
@@ -525,7 +492,7 @@ def make_fake_starlists_shifts():
 
         mnew = m + np.random.randn(N_stars) * 0.05
 
-        _out = open('random_shift_{0:d}.lis'.format(ss+1), 'w')
+        _out = open(f'{test_data_path}/random_shift_{ss+1}.lis', 'w')
         for ii in range(N_stars):
             _out.write(fmt.format(name[ii], mnew[ii], xnew[ii], ynew[ii]))
         _out.close()
@@ -563,7 +530,7 @@ def make_fake_starlists_poly1(seed=-1):
 
     # Save original positions as reference (1st) list
     # in a StarList format (with velocities).
-    lis.write('random_ref.fits', overwrite=True)
+    lis.write(f'{test_data_path}/random_ref.fits', overwrite=True)
 
     ##########
     # Shifts
@@ -614,7 +581,7 @@ def make_fake_starlists_poly1(seed=-1):
         new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t],
                                      names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't'))
 
-        new_lis.write('random_{0:d}.fits'.format(ss), overwrite=True)
+        new_lis.write(f'{test_data_path}/random_{ss}.fits', overwrite=True)
 
     return (xy_trans,mag_trans)
 
@@ -656,7 +623,7 @@ def make_fake_starlists_poly0_vel(seed=-1):
 
     # Save original positions as reference (1st) list
     # in a StarList format (with velocities).
-    lis.write('random_vel_ref.fits', overwrite=True)
+    lis.write(f'{test_data_path}/random_vel_ref.fits', overwrite=True)
     
     ##########
     # Propogate to new times and distort.
@@ -707,7 +674,7 @@ def make_fake_starlists_poly0_vel(seed=-1):
         new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t],
                                      names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't'))
 
-        new_lis.write('random_vel_p0_{0:d}.fits'.format(ss), overwrite=True)
+        new_lis.write(f'{test_data_path}/random_vel_p0_{ss}.fits', overwrite=True)
 
     return (xy_trans, mag_trans)
 
@@ -750,7 +717,7 @@ def make_fake_starlists_poly1_vel(seed=-1):
 
     # Save original positions as reference (1st) list
     # in a StarList format (with velocities).
-    lis.write('random_vel_ref.fits', overwrite=True)
+    lis.write(f'{test_data_path}/random_vel_ref.fits', overwrite=True)
     
     ##########
     # Propogate to new times and distort.
@@ -801,7 +768,7 @@ def make_fake_starlists_poly1_vel(seed=-1):
         new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t],
                                      names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't'))
 
-        new_lis.write('random_vel_{0:d}.fits'.format(ss), overwrite=True)
+        new_lis.write(f'{test_data_path}/random_vel_{ss}.fits', overwrite=True)
 
     return (xy_trans, mag_trans)
 
@@ -856,7 +823,7 @@ def make_fake_starlists_poly1_acc(seed=-1):
 
     # Save original positions as reference (1st) list
     # in a StarList format (with velocities).
-    lis.write('random_acc_ref.fits', overwrite=True)
+    lis.write(f'{test_data_path}/random_acc_ref.fits', overwrite=True)
     
     ##########
     # Propogate to new times and distort.
@@ -907,7 +874,7 @@ def make_fake_starlists_poly1_acc(seed=-1):
         new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t],
                                      names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't'))
 
-        new_lis.write('random_acc_{0:d}.fits'.format(ss), overwrite=True)
+        new_lis.write(f'{test_data_path}/random_acc_{ss}.fits', overwrite=True)
 
     return (xy_trans, mag_trans)
     
@@ -959,7 +926,7 @@ def make_fake_starlists_poly1_par(seed=-1):
 
     # Save original positions as reference (1st) list
     # in a StarList format (with velocities).
-    lis.write('random_par_ref.fits', overwrite=True)
+    lis.write(f'{test_data_path}/random_par_ref.fits', overwrite=True)
     
     ##########
     # Propogate to new times and distort.
@@ -1019,7 +986,7 @@ def make_fake_starlists_poly1_par(seed=-1):
         new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t],
                                      names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't'))
 
-        new_lis.write('random_par_{0:d}.fits'.format(ss), overwrite=True)
+        new_lis.write(f'{test_data_path}/random_par_{ss}.fits', overwrite=True)
 
     return (xy_trans, mag_trans)
 
@@ -1036,15 +1003,15 @@ def test_MosaicToRef_hst_me():
     dec = '-34:27:05.01'
     
     # Load up a Gaia catalog (queried around the RA/Dec above)
-    my_gaia = Table.read('mb10364_data/my_gaia.fits')
+    my_gaia = Table.read(f'{test_data_path}/my_gaia.fits')
     my_gaia['me'] = 0.01
     
     # Gather the list of starlists. For first pass, don't modify the starlists.
     # Loop through the observations and read them in, in prep for alignment with Gaia
     epochs = [2011.83, 2012.73, 2013.81]
-    starlist_names = ['mb10364_data/2011_10_31_F606W_MATCHUP_XYMEEE_final.calib',
-                      'mb10364_data/2012_09_25_F606W_MATCHUP_XYMEEE_final.calib',
-                      'mb10364_data/2013_10_24_F606W_MATCHUP_XYMEEE_final.calib']
+    starlist_names = [f'{test_data_path}/mb10364_data/2011_10_31_F606W_MATCHUP_XYMEEE_final.calib',
+                      f'{test_data_path}/mb10364_data/2012_09_25_F606W_MATCHUP_XYMEEE_final.calib',
+                      f'{test_data_path}/mb10364_data/2013_10_24_F606W_MATCHUP_XYMEEE_final.calib']
         
     list_of_starlists = []
     
@@ -1090,9 +1057,9 @@ def test_bootstrap():
     etc.)
     """
     # Read in starlists for MosaicToRef
-    ref = Table.read('ref_vel.lis', format='ascii')
-    list1 = Table.read('E.lis', format='ascii')
-    list2 = Table.read('F.lis', format='ascii')
+    ref = Table.read(f'{test_data_path}/ref_vel.lis', format='ascii')
+    list1 = Table.read(f'{test_data_path}/E.lis', format='ascii')
+    list2 = Table.read(f'{test_data_path}/F.lis', format='ascii')
 
     list1 = starlists.StarList.from_table(list1)
     list2 = starlists.StarList.from_table(list2)
@@ -1202,10 +1169,10 @@ def test_calc_vel_in_bootstrap():
     import copy
     
     # Define match parameters
-    ref = Table.read('ref_vel.lis', format='ascii')
+    ref = Table.read(f'{test_data_path}/ref_vel.lis', format='ascii')
 
-    list1 = Table.read('E.lis', format='ascii')
-    list2 = Table.read('F.lis', format='ascii')
+    list1 = Table.read(f'{test_data_path}/E.lis', format='ascii')
+    list2 = Table.read(f'{test_data_path}/F.lis', format='ascii')
 
     list1 = starlists.StarList.from_table(list1)
     list2 = starlists.StarList.from_table(list2)
@@ -1271,9 +1238,9 @@ def test_transform_xym():
     otherwise
     """
     #---Align 1: self.mag_Trans = False---#
-    ref = Table.read('ref_vel.lis', format='ascii')
-    list1 = Table.read('E.lis', format='ascii')
-    list2 = Table.read('F.lis', format='ascii')
+    ref = Table.read(f'{test_data_path}/ref_vel.lis', format='ascii')
+    list1 = Table.read(f'{test_data_path}/E.lis', format='ascii')
+    list2 = Table.read(f'{test_data_path}/F.lis', format='ascii')
 
     list1 = starlists.StarList.from_table(list1)
     list2 = starlists.StarList.from_table(list2)
@@ -1369,7 +1336,7 @@ def test_MosaicToRef_mag_bug():
     """
     make_fake_starlists_poly1_vel()
 
-    ref_list = starlists.StarList.read('random_vel_0.fits')
+    ref_list = starlists.StarList.read(f'{test_data_path}/random_vel_0.fits')
     lists = [ref_list]
 
     msc = align.MosaicToRef(ref_list, lists, 
@@ -1432,7 +1399,7 @@ def test_masked_cols():
     list_of_starlists = []
 
     for ee in range(len(epochs)):
-        lis_file = 'mag' + epochs[ee] + '_ob150029_kp_rms_named.lis'
+        lis_file = f'{test_data_path}/mag{epochs[ee]}_ob150029_kp_rms_named.lis'
         lis = starlists.StarList.from_lis_file(lis_file)
     
         list_of_starlists.append(lis)
@@ -1451,3 +1418,6 @@ def test_masked_cols():
     msc.fit()
 
     return
+
+if __name__ == "__main__":
+    test_MosaicSelfRef_vel()
\ No newline at end of file
diff --git a/flystar/tests/A.lis b/flystar/tests/test_data/A.lis
similarity index 100%
rename from flystar/tests/A.lis
rename to flystar/tests/test_data/A.lis
diff --git a/flystar/tests/B.lis b/flystar/tests/test_data/B.lis
similarity index 100%
rename from flystar/tests/B.lis
rename to flystar/tests/test_data/B.lis
diff --git a/flystar/tests/C.lis b/flystar/tests/test_data/C.lis
similarity index 100%
rename from flystar/tests/C.lis
rename to flystar/tests/test_data/C.lis
diff --git a/flystar/tests/D.lis b/flystar/tests/test_data/D.lis
similarity index 100%
rename from flystar/tests/D.lis
rename to flystar/tests/test_data/D.lis
diff --git a/flystar/tests/E.lis b/flystar/tests/test_data/E.lis
similarity index 100%
rename from flystar/tests/E.lis
rename to flystar/tests/test_data/E.lis
diff --git a/flystar/tests/F.lis b/flystar/tests/test_data/F.lis
similarity index 100%
rename from flystar/tests/F.lis
rename to flystar/tests/test_data/F.lis
diff --git a/flystar/tests/test_data/list_of_starlists.pkl b/flystar/tests/test_data/list_of_starlists.pkl
new file mode 100644
index 0000000..3662f0f
Binary files /dev/null and b/flystar/tests/test_data/list_of_starlists.pkl differ
diff --git a/flystar/tests/test_data/my_gaia.pkl b/flystar/tests/test_data/my_gaia.pkl
new file mode 100644
index 0000000..58fa1c8
Binary files /dev/null and b/flystar/tests/test_data/my_gaia.pkl differ
diff --git a/flystar/tests/ref.lis b/flystar/tests/test_data/ref.lis
similarity index 100%
rename from flystar/tests/ref.lis
rename to flystar/tests/test_data/ref.lis
diff --git a/flystar/tests/ref_vel.lis b/flystar/tests/test_data/ref_vel.lis
similarity index 100%
rename from flystar/tests/ref_vel.lis
rename to flystar/tests/test_data/ref_vel.lis
diff --git a/flystar/tests/test_all_detected.fits b/flystar/tests/test_data/test_all_detected.fits
similarity index 100%
rename from flystar/tests/test_all_detected.fits
rename to flystar/tests/test_data/test_all_detected.fits
diff --git a/flystar/tests/test_catalog.fits b/flystar/tests/test_data/test_catalog.fits
similarity index 100%
rename from flystar/tests/test_catalog.fits
rename to flystar/tests/test_data/test_catalog.fits