diff --git a/.travis.yml b/.travis.yml
index 13917ed..36a1c7d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -16,10 +16,10 @@ python:
 
 # Setup anaconda
 before_install:
-  - wget http://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh
+  - wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
   - chmod +x miniconda.sh
   - ./miniconda.sh -b
-  - export PATH=/home/travis/miniconda2/bin:$PATH
+  - export PATH=/home/travis/miniconda3/bin:$PATH
   - ls
   - conda update --yes conda
   # The next couple lines fix a crash with multiprocessing on Travis and are not specific to using Miniconda
diff --git a/aotools/aotools.py b/aotools/aotools.py
index 7d78dc5..596d3c5 100644
--- a/aotools/aotools.py
+++ b/aotools/aotools.py
@@ -40,11 +40,26 @@ def image_contrast(image):
 
     Parameters:
         image (ndarray): Image array
-
+ 
     Returns:
         float: Contrast value
     """
 
     contrast = (image.max() - image.min()) / (image.max() + image.min())
 
-    return contrast
\ No newline at end of file
+    return contrast
+
+def rms_contrast(image):
+    """
+    Calculates the RMS contrast - basically the standard deviation of the image
+
+    Parameters:
+        image (ndarray): Image array
+
+    Returns:
+        float: Contrast value
+    """
+    
+    image /= image.max()
+    
+    return image.std()
diff --git a/aotools/centroiders/centroiders.py b/aotools/centroiders/centroiders.py
index b827412..6d90117 100644
--- a/aotools/centroiders/centroiders.py
+++ b/aotools/centroiders/centroiders.py
@@ -35,7 +35,7 @@ def centreOfGravity(img, threshold=0, **kwargs):
     y_centroid+=0.5
     x_centroid+=0.5
 
-    return numpy.array([y_centroid,x_centroid])
+    return numpy.array([x_centroid, y_centroid])
 
 
 def brightestPxl(img, threshold, **kwargs):
diff --git a/aotools/turbulence/infinitephasescreen.py b/aotools/turbulence/infinitephasescreen.py
index fa97766..fd39c10 100644
--- a/aotools/turbulence/infinitephasescreen.py
+++ b/aotools/turbulence/infinitephasescreen.py
@@ -181,7 +181,7 @@ def makeAMatrix(self):
         # Different inversion methods, not sure which is best
         cf = linalg.cho_factor(self.cov_zz)
         inv_cov_zz = linalg.cho_solve(cf, numpy.identity(self.cov_zz.shape[0]))
-        # inv_cov_zz = numpy.linalg.pinv(self.cov_zz)#, 0.001)    
+        # inv_cov_zz = numpy.linalg.pinv(self.cov_zz)
         
         self.A_mat_forwards = self.cov_xz_forwards.dot(inv_cov_zz)
         self.A_mat_backwards = self.cov_xz_backwards.dot(inv_cov_zz) 
@@ -536,7 +536,7 @@ def phaseCovariance(r, r0, L0):
     
     A = (L0/r0)**(5./3) 
     
-    B1 = (2**(1./6)) * gamma(11./6)/(pi**(8./3))
+    B1 = (2**(-5./6)) * gamma(11./6)/(pi**(8./3))
     B2 = ((24./5) * gamma(6./5))**(5./6)
     
     C = (((2 * pi * r)/L0) ** (5./6)) * kv(5./6, (2 * pi * r)/L0)
diff --git a/aotools/turbulence/phasescreen.py b/aotools/turbulence/phasescreen.py
index 4f8cdec..18690a0 100644
--- a/aotools/turbulence/phasescreen.py
+++ b/aotools/turbulence/phasescreen.py
@@ -138,7 +138,7 @@ def ft_phase_screen(r0, N, delta, L0, l0, FFT=None):
     PSD_phi  = (0.023*r0**(-5./3.) * numpy.exp(-1*((f/fm)**2)) /
                 ( ( (f**2) + (f0**2) )**(11./6) ) )
 
-    PSD_phi[(N/2),(N/2)] = 0
+    PSD_phi[int(N/2), int(N/2)] = 0
 
     cn = ( (numpy.random.normal(size=(N,N)) + 1j* numpy.random.normal(size=(N,N)) )
                 * numpy.sqrt(PSD_phi)*del_f )
diff --git a/aotools/turbulence/temporal_ps.py b/aotools/turbulence/temporal_ps.py
index b6e0d39..a4e5b48 100644
--- a/aotools/turbulence/temporal_ps.py
+++ b/aotools/turbulence/temporal_ps.py
@@ -7,6 +7,7 @@
 
 import numpy
 from matplotlib import pyplot
+from scipy import optimize
 
 
 def calc_slope_temporalps(slope_data):
@@ -34,9 +35,11 @@ def calc_slope_temporalps(slope_data):
     tps = abs(numpy.fft.fft(slope_data, axis=-2)[..., :n_frames/2, :])**2
 
     # Find mean across all sub-aps
-    mean_tps = (abs(tps)**2).mean(-1)
+    tps = (abs(tps)**2)
+    mean_tps = tps.mean(-1)
+    tps_err = tps.std(-1)/numpy.sqrt(tps.shape[-1])
 
-    return mean_tps
+    return mean_tps, tps_err
 
 
 def get_tps_time_axis(frame_rate, n_frames):
@@ -64,11 +67,12 @@ def plot_tps(slope_data, frame_rate):
         frame_rate (float): Frame rate of detector observing slope gradients (Hz)
 
     Returns:
+        tuple: The computed temporal power spectrum/a, and the time axis data
 
     """
     n_frames = slope_data.shape[-2]
 
-    tps = calc_slope_temporalps(slope_data)
+    tps, tps_err = calc_slope_temporalps(slope_data)
 
     t_axis_data = get_tps_time_axis(frame_rate, n_frames)
 
@@ -77,7 +81,7 @@ def plot_tps(slope_data, frame_rate):
 
     # plot each power spectrum
     for i_ps, ps in enumerate(tps):
-        ax.semilogy(t_axis_data, ps, label="Spectrum {}".format(i_ps))
+        ax.loglog(t_axis_data, ps, label="Spectrum {}".format(i_ps))
 
     ax.legend()
 
@@ -86,13 +90,104 @@ def plot_tps(slope_data, frame_rate):
 
     pyplot.show()
 
-def fit_tps(tps, axis):
+    return tps, tps_err, t_axis_data
+
+def fit_tps(tps, t_axis, D, V_guess=10, f_noise_guess=20, A_guess=9, tps_err=None, plot=False):
     """
+    Runs minimization routines to get t0.
+
     Parameters:
-        tps: The temporal power spectrum to fit
-        axis: fit parallel or perpendicular to wind direction
+        tps (ndarray): The temporal power spectrum to fit
+        axis (str): fit parallel ('par') or perpendicular ('per') to wind direction
+        D (float): (Sub-)Aperture diameter
 
     Returns:
+    """
+    if plot:
+        fig = pyplot.figure()
+
+    opt_result = optimize.minimize(
+            test_tps_fit_minimize_func,
+            x0=(V_guess, f_noise_guess, A_guess),
+            args=(tps, t_axis, D, tps_err, plot),
+            method="COBYLA")
+
+    print(opt_result)
+
+
+def test_tps_fit(tps, t_axis_data, D, V, f_noise, A=1, tps_err=None, plot=False):
+    """
+    Tests the validaty of a fit to the temporal power spectrum.
+
+    Uses the temporal power spectrum and time-axis data to test the validity of a coherence time. A frequency above which fitting is not performaed should also be given, as noise will be the dominant contributor above this.
 
+    Parameters:
+        tps (ndarray): Temporal power spectrum to fit
+        t_axis_data (ndarray): Time axis data
+        D (float): (sub-) Aperture diameter
+        V (float): Integrated wind speed
+        f_noise (float): Frequency above which noise dominates.
+        A (float): Initial Guess of
     """
+    f0 = 0.3 * V/D
+
+    if f0<t_axis_data[0] or f0>f_noise or f_noise>t_axis_data.max():
+        return 10**99
+
+    tps_tt_indices = numpy.where((t_axis_data<f0) & (t_axis_data>0))[0]
+    tt_t_data = t_axis_data[tps_tt_indices]
+    tt_fit = 10**A * tt_t_data**(-2./3)
+
+    # get scaling for next part of distribution so it matches up at cutof freq.
+    tps_ho_indices = numpy.where((t_axis_data>f0) & (t_axis_data<f_noise))
+    ho_t_data = t_axis_data[tps_ho_indices]
+    B = tt_fit[-1]/(ho_t_data[0] ** (-11./3))
+    ho_fit = B * ho_t_data ** (-11./3)
+
+    ps_fit = numpy.append(tt_fit, ho_fit)
+    fit_t_data = numpy.append(tt_t_data, ho_t_data)
+    fit_t_coords = numpy.append(tps_tt_indices, tps_ho_indices)
+
+
+    if tps_err is None:
+        err = numpy.mean((numpy.sqrt(numpy.square(ps_fit - tps[fit_t_coords]))))
+    else:
+        chi2 = numpy.square((ps_fit - tps[fit_t_coords])/tps_err[fit_t_coords]).sum()
+        err = chi2/ps_fit.shape[0]
+
+
+
+    if plot:
+        print("V: {}, f_noise: {}, A: {}".format(V, f_noise, A))
+        pyplot.cla()
+        ax = pyplot.gca()
+        ax.set_xscale('log')
+        ax.set_yscale('log')
+        if tps_err is None:
+            ax.loglog(t_axis_data, tps)
+        else:
+            ax.errorbar(t_axis_data, tps, tps_err)
+        ax.plot(fit_t_data, ps_fit, linewidth=2)
+
+        ax.plot([f0]*2, [0.1, tps.max()], color='k')
+        pyplot.draw()
+        pyplot.pause(0.01)
+        print("Error Func: {}".format(err))
+
+    return err
+
+def test_tps_fit_minimize_func(args, tps, t_axis, D, tps_err=None, plot=False):
+
+    V, f_noise, A = args
+
+    return test_tps_fit(tps, t_axis, D=D, V=V, f_noise=f_noise, A=A, tps_err=None, plot=plot)
+
+
+if __name__ == "__main__":
+
+    from astropy.io import fits
+
+    data = fits.getdata("t0_5ms_canary_slopes.fits")
 
+    tps, tps_err, t_data = plot_tps(data, 150)
+    fit_tps(tps[0], t_data, D=4.2/7, V_guess=20, f_noise_guess=50, A_guess=9, tps_err=tps_err[0], plot=True)
\ No newline at end of file
diff --git a/aotools/wfs/wfslib.py b/aotools/wfs/wfslib.py
index 9d7ae7d..c6159a4 100644
--- a/aotools/wfs/wfslib.py
+++ b/aotools/wfs/wfslib.py
@@ -4,6 +4,12 @@
 
 import numpy
 
+# Best range for python2 and 3
+try:
+    range = xrange
+except:
+    pass
+
 def r0fromSlopes(slopes, wavelength, subapDiam):
     """
     Measures the value of R0 from a set of WFS slopes.
@@ -92,4 +98,27 @@ def computeFillFactor(mask, subapPos, subapSpacing):
         y2 = int(round(y + subapSpacing))
         fills[i] = mask[x1:x2, y1:y2].mean()
 
-    return fills
\ No newline at end of file
+    return fills
+
+def make_subaps_2d(data, mask):
+    """
+    Fills in a pupil shape with 2-d sub-apertures
+
+    Parameters:
+        data (ndarray): slope data, of shape, (frames, 2, nSubaps)
+        mask (ndarray): 2-d array of shape (nxSubaps, nxSubaps), where 1 indicates a valid subap position and 0 a masked subap 
+    """ 
+    n_frames = data.shape[0]
+    n_subaps = data.shape[-1]
+    nx_subaps = mask.shape[0]
+
+    subaps_2d = numpy.zeros((data.shape[0], 2, mask.shape[0], mask.shape[1]), dtype=data.dtype)
+
+    n_subap = 0
+    for x in range(nx_subaps):
+        for y in range(nx_subaps):
+            if mask[x, y] == 1:
+                subaps_2d[:, :, x, y] = data[:, :, n_subap]
+                n_subap += 1
+
+    return subaps_2d