diff --git a/flystar/align.py b/flystar/align.py
index f5f0742..994a3b1 100755
--- a/flystar/align.py
+++ b/flystar/align.py
@@ -137,10 +137,12 @@ def = None. If not None, then this should contain an array or list of transform
             Dict of motion model name keys (strings) and corresponding MotionModel object values
 
         use_scipy : bool, optional
-            If True, use scipy.optimize.curve_fit for velocity fitting. If False, use linear algebra fitting, by default True.
+            If True, use scipy.optimize.curve_fit for velocity fitting. If False, use linear
+            algebra fitting of posible, by default True.
 
         absolute_sigma : bool, optional
-            If True, the velocity fit will use absolute errors in the data. If False, relative errors will be used, by default False.
+            If True, the velocity fit will use absolute errors in the data. If False, relative
+            errors will be used, by default False.
 
         save_path : str, optional
             Path to save the MosaicSelfRef object as a pickle file.
diff --git a/flystar/motion_model.py b/flystar/motion_model.py
index 66acce6..e24067b 100644
--- a/flystar/motion_model.py
+++ b/flystar/motion_model.py
@@ -430,20 +430,19 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None,
         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(t, x0,vx, y0,vy, pi):
-            use_t = t[:int(len(t)/2)]
+        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]
-            return np.append(x_res, y_res)
+            return np.hstack([x_res, y_res])
         # Initial guesses, x0,y0 as x,y averages;
         #     vx,vy as average velocity if first and last points are perfectly measured;
         #     pi for 10 pc disance
         if params_guess is None:
             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]), 1]
-        res = curve_fit(fit_func, np.append(t,t), np.append(x,y),
-                        p0=params_guess, sigma = 1.0/np.append(x_wt,y_wt))
+                            y.mean(),(y[idx_last]-y[idx_first])/(t[idx_last]-t[idx_first]), 0.1]
+        res = curve_fit(fit_func, t, np.hstack([x,y]),
+                        p0=params_guess, sigma = 1.0/np.hstack([x_wt,y_wt]))
         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)
 
diff --git a/flystar/parallax.py b/flystar/parallax.py
index 27cd632..b4b3a1b 100755
--- a/flystar/parallax.py
+++ b/flystar/parallax.py
@@ -1,3 +1,6 @@
+# Parallax calculation module for motion models involving parallax
+# Adapted from BAGLE's parallax.py
+
 import math
 
 import numpy as np