diff --git a/README.md b/README.md
index bb1b3b4..93f3b37 100644
--- a/README.md
+++ b/README.md
@@ -22,7 +22,17 @@ pip install astrosource
 
 ### Install development version
 
-If you need to install the development branch, download from [GitHub](https://github.com/zemogle/astrosource) and from the root of the repo, run:
+The development version tends to be more cutting edge with hotfixes applied to latest problems and with new features that are new enough that we are still monitoring for bugs. Having said that, it is probably the better choice to use, just be aware that there may be bugs and if you find one, let us know! 
+
+If you need to install the development branch, there are two methods:
+
+(1) Use pip to install the development branch
+
+```bash
+pip install -U git+https://github.com/zemogle/astrosource@dev
+```
+
+(2) download from [GitHub](https://github.com/zemogle/astrosource/tree/dev) and from the root of the repo, run:
 
 ```bash
 cd astrosource
diff --git a/astrosource/analyse.py b/astrosource/analyse.py
index 670baef..a4e2700 100644
--- a/astrosource/analyse.py
+++ b/astrosource/analyse.py
@@ -13,19 +13,23 @@
 from tqdm import tqdm
 #import traceback
 import logging
-from multiprocessing import Pool
-
+import multiprocessing as mp
+from multiprocessing import Pool, cpu_count, shared_memory
+import traceback
 #from astrosource.utils import photometry_files_to_array, AstrosourceException
 from astrosource.utils import AstrosourceException
 from astrosource.plots import plot_variability
+from astropy.stats import sigma_clip
 
+from scipy.optimize import curve_fit
 import matplotlib
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import matplotlib.colors as colors
-
+from functools import partial
 logger = logging.getLogger('astrosource')
 
+import platform
 
 def get_total_counts(photFileArray, compFile, loopLength, photCoords):
 
@@ -37,7 +41,6 @@ def get_total_counts(photFileArray, compFile, loopLength, photCoords):
     for photFile in photFileArray:
         allCounts = 0.0
         allCountsErr = 0.0
-        #fileRaDec = SkyCoord(ra=photFile[:, 0]*degree, dec=photFile[:, 1]*degree)
         fileRaDec = photCoords[counter]
         counter=counter+1
         #Array of comp measurements
@@ -53,10 +56,14 @@ def get_total_counts(photFileArray, compFile, loopLength, photCoords):
             if (compFile.shape[0] == 5 and compFile.size == 5) or (compFile.shape[0] == 3 and compFile.size == 3):
                 break
         allCountsArray.append([allCounts, allCountsErr])
-    #logger.debug(allCountsArray)
     return allCountsArray
 
-def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs):
+#def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs):
+def process_varsearch_target(target, photFileArray_shape, photFileArray_dtype, shm_name, allCountsArray, matchRadius, minimumNoOfObs):
+    
+    # Attach to the shared memory for photFileArray
+    existing_shm = shared_memory.SharedMemory(name=shm_name)
+    photFileArray = np.ndarray(photFileArray_shape, dtype=photFileArray_dtype, buffer=existing_shm.buf)    
     diffMagHolder = []
     
     for allcountscount, photFile in enumerate(photFileArray):
@@ -82,30 +89,239 @@ def process_varsearch_target(target, photFileArray, allCountsArray, matchRadius,
         if diffMagHolder.size == sizeBefore:
             break
 
+    existing_shm.close()
+
     # Append to output if sufficient observations are available
     if diffMagHolder.size > minimumNoOfObs:
         return [target[0], target[1], np.median(diffMagHolder), np.std(diffMagHolder), diffMagHolder.size]
     return None
 
+def process_varsearch_target_singleprocess(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs):    
+        
+    diffMagHolder = []
+    
+    for allcountscount, photFile in enumerate(photFileArray):
+        # Efficiently calculate the closest match using numpy
+        idx = (np.abs(photFile[:, 0] - target[0]) + np.abs(photFile[:, 1] - target[1])).argmin()
+        d2d = np.sqrt((photFile[idx, 0] - target[0])**2 + (photFile[idx, 1] - target[1])**2) * 3600
+        multTemp = -2.5 * np.log10(photFile[idx, 11] / allCountsArray[allcountscount][0])
+        
+        if d2d < matchRadius and not np.isinf(multTemp):
+            diffMagHolder.append(multTemp)
+    
+    # Remove major outliers
+    diffMagHolder = np.array(diffMagHolder)
+    while True:
+        stdVar = np.std(diffMagHolder)
+        avgVar = np.mean(diffMagHolder)
+        sizeBefore = diffMagHolder.size
+
+        # Mask outliers
+        mask = (diffMagHolder <= avgVar + 4 * stdVar) & (diffMagHolder >= avgVar - 4 * stdVar)
+        diffMagHolder = diffMagHolder[mask]
+
+        if diffMagHolder.size == sizeBefore:
+            break
+
+    # Append to output if sufficient observations are available
+    if diffMagHolder.size > minimumNoOfObs:
+        return [target[0], target[1], np.median(diffMagHolder), np.std(diffMagHolder), diffMagHolder.size]
+    
+    return None
+
 
 def process_varsearch_targets_multiprocessing(targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs):
+    # Create shared memory for photFileArray
+    # As this can lead to gigantic RAM use for large datasets
+    shm = shared_memory.SharedMemory(create=True, size=photFileArray.nbytes)
+    shared_photFileArray = np.ndarray(photFileArray.shape, dtype=photFileArray.dtype, buffer=shm.buf)
+    shared_photFileArray[:] = photFileArray[:]
+    
+    
     # Partial function to pass shared arguments
-    from functools import partial
+    
     worker = partial(
         process_varsearch_target,
-        photFileArray=photFileArray,
+        photFileArray_shape=photFileArray.shape,
+        photFileArray_dtype=photFileArray.dtype,
+        shm_name=shm.name,               
         allCountsArray=allCountsArray,
         matchRadius=matchRadius,
         minimumNoOfObs=minimumNoOfObs,
     )
     
     # Multiprocessing with Pool
-    with Pool() as pool:
+    with Pool(processes=max([cpu_count()-1,1])) as pool:
         results = pool.map(worker, targetFile)
     
+    # Clean up shared memory
+    shm.close()
+    shm.unlink()
+    
     # Filter out None results
     return [res for res in results if res is not None]
 
+
+def fit_sigma_clipped_poly(x, y, order=3, sigma=2, parentPath=''):
+    """
+    Fits a sigma-clipped polynomial to the data and identifies outliers.
+
+    Parameters:
+    x (array-like): The x-values of the data points.
+    y (array-like): The y-values of the data points.
+    order (int): The order of the polynomial to fit (default: 3).
+    sigma (float): The number of standard deviations for sigma clipping (default: 2).
+
+    Returns:
+    None (plots the data, fitted polynomial, and identified outliers).
+    """
+    # Ensure inputs are numpy arrays
+    x = np.array(x)
+    y = np.array(y)
+
+    # Sigma clipping
+    clipped_data = sigma_clip(y, sigma=sigma, maxiters=5)
+    mask = clipped_data.mask
+
+    # Fit polynomial to non-masked data
+    poly_coeff = np.polyfit(x[~mask], y[~mask], order)
+    poly_func = np.poly1d(poly_coeff)
+
+    # Calculate residuals and standard deviation
+    y_fit = poly_func(x)
+    residuals = y - y_fit
+    std_dev = np.std(residuals[~mask])
+
+    # Identify outliers
+    outliers = residuals > (2 * std_dev)
+
+    # Plot results
+    plt.figure(figsize=(10, 6))
+    plt.scatter(x, y, color='blue', label='Data')
+    plt.plot(x, y_fit, color='red', label=f'{order}-order Polynomial Fit')
+    plt.scatter(x[outliers], y[outliers], color='orange', label='Outliers (>2 std dev)', zorder=5)
+    plt.xlabel('x')
+    plt.ylabel('y')
+    plt.title('Sigma-Clipped Polynomial Fit and Outliers')
+    plt.legend()
+    plt.grid(True)
+    #plt.show()
+    
+    plt.savefig(parentPath / "results/polifitdiagram.png")
+    
+def sigmoid_func(x, L, k, x0):
+    """Sigmoid function: y = L / (1 + exp(-k * (x - x0))) + 0.01"""
+    return L / (1 + np.exp(-k * (x - x0))) + 0.01
+
+def fit_sigma_clipped_sigmoid(x, y, sigma=2, parentPath=''):
+    """
+    Fits a sigma-clipped sigmoid function to the data with weighted fitting and identifies outliers.
+
+    Parameters:
+    x (array-like): The x-values of the data points.
+    y (array-like): The y-values of the data points.
+    sigma (float): The number of standard deviations for sigma clipping (default: 2).
+
+    Returns:
+    outlier_indices (array): Indices of the identified outliers.
+    """
+    # Ensure inputs are numpy arrays
+    x = np.array(x)
+    y = np.array(y)
+
+    # Sigma clipping
+    clipped_data = sigma_clip(y, sigma=sigma, maxiters=5)
+    mask = clipped_data.mask
+
+    # Create weights that emphasize the middle of the x-range
+    weights = np.exp(-((x - np.median(x))**2) / (2 * (np.std(x) / 2)**2))
+
+    # Fit sigmoid function to non-masked data with weights
+    p0 = [max(y), 1, np.median(x)]  # Initial guesses for L, k, x0
+    try:
+        popt, _ = curve_fit(sigmoid_func, x[~mask], y[~mask], p0=p0, sigma=1/weights[~mask])
+    except:
+        logger.info("Sigmoid curve fit failed")
+        return []
+
+    # Calculate residuals and standard deviation
+    y_fit = sigmoid_func(x, *popt)
+    residuals = y - y_fit
+    std_dev = np.std(residuals[~mask])
+
+    # Identify outliers
+    outliers = residuals > (2 * std_dev)
+    outlier_indices = np.where(outliers)[0]
+
+    # Plot results
+    plt.figure(figsize=(10, 6))
+    plt.scatter(x, y, color='blue', label='Data')
+    plt.plot(x, y_fit, color='red', label='Sigmoid Fit')
+    plt.scatter(x[outliers], y[outliers], color='orange', label='Outliers (>2 std dev)', zorder=5)
+    plt.xlabel('x')
+    plt.ylabel('y')
+    plt.title('Weighted Sigma-Clipped Sigmoid Fit and Outliers')
+    plt.legend()
+    plt.grid(True)
+    plt.show()
+    
+    plt.savefig(parentPath / "results/polifitdiagram.png")
+
+    return outlier_indices
+
+def exponential_func(x, a, b, c):
+    """Exponential function: y = a * exp(b * x) + c"""
+    return a * np.exp(b * x) + c
+
+def fit_sigma_clipped_exp(x, y, sigma=2, parentPath=''):
+    """
+    Fits a sigma-clipped exponential function to the data and identifies outliers.
+
+    Parameters:
+    x (array-like): The x-values of the data points.
+    y (array-like): The y-values of the data points.
+    sigma (float): The number of standard deviations for sigma clipping (default: 2).
+
+    Returns:
+    outlier_indices (array): Indices of the identified outliers.
+    """
+    # Ensure inputs are numpy arrays
+    x = np.array(x)
+    y = np.array(y)
+
+    # Sigma clipping
+    clipped_data = sigma_clip(y, sigma=sigma, maxiters=5)
+    mask = clipped_data.mask
+
+    # Fit exponential function to non-masked data
+    popt, _ = curve_fit(exponential_func, x[~mask], y[~mask], p0=(1, 0.1, 1))
+
+    # Calculate residuals and standard deviation
+    y_fit = exponential_func(x, *popt)
+    residuals = y - y_fit
+    std_dev = np.std(residuals[~mask])
+
+    # Identify outliers
+    outliers = residuals > (2 * std_dev)
+    outlier_indices = np.where(outliers)[0]
+
+    # Plot results
+    plt.figure(figsize=(10, 6))
+    plt.scatter(x, y, color='blue', label='Data')
+    plt.plot(x, y_fit, color='red', label='Exponential Fit')
+    plt.scatter(x[outliers], y[outliers], color='orange', label='Outliers (>2 std dev)', zorder=5)
+    plt.xlabel('x')
+    plt.ylabel('y')
+    plt.title('Sigma-Clipped Exponential Fit and Outliers')
+    plt.legend()
+    plt.grid(True)
+    plt.show()
+    
+    plt.savefig(parentPath / "results/polifitdiagram.png")
+
+    return outlier_indices
+
+
 def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None, varsearchglobalstdev=-99.9, varsearchthresh=10000, varsearchstdev=2.0, varsearchmagwidth=0.25, varsearchminimages=0.3, photCoords=None, photFileHolder=None, fileList=None):
     '''
     Find stable comparison stars for the target photometry and remove variables
@@ -126,7 +342,7 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None,
 
 
     photFileArray=photFileHolder
-    photFileCoords=photCoords
+    #photFileCoords=photCoords
 
     # allocate minimum images to detect
     minimumNoOfObs=int(varsearchminimages*len(fileList))
@@ -149,7 +365,7 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None,
 
     compFile = genfromtxt(parentPath / "results/compsUsed.csv", dtype=float, delimiter=',')
     logger.debug("Stable Comparison Candidates below variability threshold")
-    outputPhot = []
+    #outputPhot = []
 
     # Get total counts for each file
     allCountsArray = get_total_counts(photFileArray, compFile, loopLength=compFile.shape[0], photCoords=photCoords)
@@ -168,51 +384,28 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None,
 
     ## NEED TO REMOVE COMPARISON STARS FROM TARGETLIST
 
-    allcountscount=0
+    #allcountscount=0
     # For each variable calculate the variability
     outputVariableHolder=[]
 
     logger.info("Measuring variability of stars...... ")
     taketime=time.time()
 
-    # for target in targetFile:
-    #     diffMagHolder=[]
-    #     allcountscount=0
-
+    # Hack to get windows to not multiprocess until a good solution is found.    
+    if platform.system() == "Windows":       
         
-    #     for photFile in photFileArray:
-    #         # A bit rougher than using SkyCoord, but way faster
-    #         # The amount of calculations is too slow for SkyCoord
-    #         idx=(np.abs(photFile[:,0] - target[0]) + np.abs(photFile[:,1] - target[1])).argmin()
-    #         d2d=pow(pow(photFile[idx,0] - target[0],2) + pow(photFile[idx,1] - target[1],2),0.5) * 3600
-    #         multTemp=(multiply(-2.5,log10(divide(photFile[idx][11],allCountsArray[allcountscount][0]))))
-    #         if less(d2d, matchRadius) and (multTemp != inf) :
-    #             diffMagHolder=append(diffMagHolder,multTemp)
-    #         allcountscount=add(allcountscount,1)
-            
-    #     ## REMOVE MAJOR OUTLIERS FROM CONSIDERATION
-    #     diffMagHolder=np.array(diffMagHolder)
-    #     while True:
-    #         stdVar=std(diffMagHolder)
-    #         avgVar=average(diffMagHolder)
-
-    #         sizeBefore=diffMagHolder.shape[0]
-    #         #print (sizeBefore)
-
-    #         diffMagHolder[diffMagHolder > avgVar+(4*stdVar) ] = np.nan
-    #         diffMagHolder[diffMagHolder < avgVar-(4*stdVar) ] = np.nan
-    #         diffMagHolder=diffMagHolder[~np.isnan(diffMagHolder)]
-
-    #         if diffMagHolder.shape[0] == sizeBefore:
-    #             break
-
-    #     if (diffMagHolder.shape[0] > minimumNoOfObs):
-    #         outputVariableHolder.append( [target[0],target[1],median(diffMagHolder), std(diffMagHolder), diffMagHolder.shape[0]])
-
-# process_varsearch_target
-    outputVariableHolder = process_varsearch_targets_multiprocessing(
-        targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs
-    )
+        outputVariableHolder=[]
+        for target in targetFile:
+            result= process_varsearch_target_singleprocess(target, photFileArray, allCountsArray, matchRadius, minimumNoOfObs)
+            if result == None:
+                pass
+            else:
+                outputVariableHolder.append(result)
+    else:
+        
+        outputVariableHolder = process_varsearch_targets_multiprocessing(
+            targetFile, photFileArray, allCountsArray, matchRadius, minimumNoOfObs
+        )
 
     print ("Star Variability done in " + str(time.time()-taketime))
 
@@ -222,6 +415,10 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None,
     starVar = np.asarray(outputVariableHolder)
 
     
+    outliers=fit_sigma_clipped_sigmoid(starVar[:,2],starVar[:,3], parentPath=parentPath)
+
+    potentialVariables=starVar[outliers]
+
     meanMags = starVar[:,2]
     variations = starVar[:,3]
 
@@ -230,51 +427,29 @@ def find_variable_stars(targets, matchRadius, errorReject=0.05, parentPath=None,
     xbins = np.arange(np.min(meanMags), np.max(meanMags), xStepSize)
     ybins = np.arange(np.min(variations), np.max(variations), yStepSize)
 
-    #split it into one array and identify variables in bins with centre
-    variationsByMag=[]
-    potentialVariables=[]
-    for xbinner in range(len (xbins)):
-
-        starsWithin=[]
-        for q in range(len(meanMags)):
-            if meanMags[q] >= xbins[xbinner] and meanMags[q] < xbins[xbinner]+xStepSize:
-                starsWithin.append(variations[q])
-
-        meanStarsWithin= (np.mean(starsWithin))
-        stdStarsWithin= (np.std(starsWithin))
-        variationsByMag.append([xbins[xbinner]+0.5*xStepSize,meanStarsWithin,stdStarsWithin])
-
-        # At this point extract RA and Dec of stars that may be variable
-        for q in range(len(starVar[:,2])):
-            if starVar[q,2] >= xbins[xbinner] and starVar[q,2] < xbins[xbinner]+xStepSize:
-                if varsearchglobalstdev != -99.9:
-                    if starVar[q,3] > varsearchglobalstdev :
-                        potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]])
-                elif starVar[q,3] > (meanStarsWithin + varsearchstdev*stdStarsWithin):
-                    #print (starVar[q,3])
-                    potentialVariables.append([starVar[q,0],starVar[q,1],starVar[q,2],starVar[q,3]])
-
-    potentialVariables=np.array(potentialVariables)
-    logger.debug("Potential Variables Identified: " + str(potentialVariables.shape[0]))
-
     if potentialVariables.shape[0] == 0:
         logger.info("No Potential Variables identified in this set of data using the parameters requested.")
     else:
         savetxt(parentPath / "results/potentialVariables.csv", potentialVariables , delimiter=",", fmt='%0.8f', header='RA,DEC,DiffMag,Variability')
 
+    try:
         plot_variability(outputVariableHolder, potentialVariables, parentPath, compFile)
-
-        plt.cla()
-        fig, ax = plt.subplots(figsize =(10, 7))
-        plt.hist2d(meanMags, variations,  bins =[xbins, ybins], norm=colors.LogNorm(), cmap = plt.cm.YlOrRd)
-        plt.colorbar()
-        plt.title("Variation Histogram")
-        ax.set_xlabel('Mean Differential Magnitude')
-        ax.set_ylabel('Variation (Standard Deviation)')
-        plt.plot(potentialVariables[:,2],potentialVariables[:,3],'bo')
-        plt.tight_layout()
-
-        plt.savefig(parentPath / "results/Variation2DHistogram.png")
+    except:
+        print ("MTF hunting this bug")
+        logger.error(traceback.print_exc())
+        breakpoint()
+
+    plt.cla()
+    fig, ax = plt.subplots(figsize =(10, 7))
+    plt.hist2d(meanMags, variations,  bins =[xbins, ybins], norm=colors.LogNorm(), cmap = plt.cm.YlOrRd)
+    plt.colorbar()
+    plt.title("Variation Histogram")
+    ax.set_xlabel('Mean Differential Magnitude')
+    ax.set_ylabel('Variation (Standard Deviation)')
+    plt.plot(potentialVariables[:,2],potentialVariables[:,3],'bo')
+    plt.tight_layout()
+
+    plt.savefig(parentPath / "results/Variation2DHistogram.png")
 
     return outputVariableHolder
 
@@ -324,15 +499,15 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file
         starDistanceRejCount=0
         logger.debug("****************************")
         logger.debug("Processing Variable {}".format(q+1))
-        if int(len(targets)) == 4 and targets.size==4:
+        if int(len(targets)) == 5 and targets.size==5:
             logger.debug("RA {}".format(targets[0]))
         else:
             logger.debug("RA {}".format(targets[q][0]))
-        if int(len(targets)) == 4 and targets.size==4:
+        if int(len(targets)) == 5 and targets.size==5:
             logger.debug("Dec {}".format(targets[1]))
         else:
             logger.debug("Dec {}".format(targets[q][1]))
-        if int(len(targets)) == 4 and targets.size==4:
+        if int(len(targets)) == 5 and targets.size==5:
             #varCoord = SkyCoord(targets[0],(targets[1]), frame='icrs', unit=degree) # Need to remove target stars from consideration
             singleCoordRA=targets[0]
             singleCoordDEC=targets[1]
@@ -347,20 +522,11 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file
         allcountscount=0
 
         for imgs, photFile in enumerate(tqdm(photFileArray)):
-            #fileRaDec = photCoordsFile[imgs]
-            #idx, d2d, _ = varCoord.match_to_catalog_sky(fileRaDec)
-            #breakpoint()
             idx=(np.abs(photFile[:,0] - singleCoordRA) + np.abs(photFile[:,1] - singleCoordDEC)).argmin()
             d2d=pow(pow(photFile[idx,0] - singleCoordRA,2) + pow(photFile[idx,1] - singleCoordDEC,2),0.5) * 3600
             
-            #print (d2d)
-            
             starRejected=0
-            if (less(d2d, targetRadius)):
-                #magErrVar = 1.0857 * (photFile[idx][5]/photFile[idx][4])
-                
-                
-                
+            if (less(d2d, targetRadius)):                
                 # If the file hasn't been calibrated, then it still contains the countrate in it.
                 # So convert these to mags, otherwise use the calibrated error.
                 if photFile[idx][4] > 50:
@@ -456,8 +622,6 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file
                     fileCount.append(allCountsArray[allcountscount][0])
                     allcountscount=allcountscount+1
 
-        #breakpoint()
-
         # Check for dud images
         imageReject=[]
         for j in range(asarray(outputPhot).shape[0]):
@@ -488,7 +652,6 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file
                     break
 
             # Reject by outsized error elimination
-
             while True:
                 errorsArray=[]
                 for j in range(asarray(outputPhot).shape[0]):
@@ -518,12 +681,8 @@ def photometric_calculations(targets, paths, targetRadius, errorReject=0.1, file
             outputPhot=delete(outputPhot, starReject, axis=0)
 
         except ValueError:
-            #raise AstrosourceException("No target stars were detected in your dataset. Check your input target(s) RA/Dec")
-            import traceback; logger.error(traceback.print_exc())
+            logger.error(traceback.print_exc())
             logger.error("This target star was not detected in your dataset. Check your input target(s) RA/Dec")
-            #logger.info("Rejected Stdev Measurements: : {}".format(stdevReject))
-            #logger.error("Rejected Error Measurements: : {}".format(starErrorRejCount))
-            #logger.error("Rejected Distance Measurements: : {}".format(starDistanceRejCount))
 
         # Add calibration columns
         outputPhot= np.c_[outputPhot, np.ones(outputPhot.shape[0]),np.ones(outputPhot.shape[0])]
@@ -606,4 +765,8 @@ def calibrated_photometry(paths, photometrydata, colourterm, colourerror, colour
         logger.info("then consider providing an appropriate colour for this filter using the --targetcolour option")
         logger.info("as well as an appropriate colour term for this filter (using --colourdetect or --colourterm).")
 
-    return pdata
\ No newline at end of file
+    return pdata
+
+# Needed for windows to multiprocess appropriately
+if __name__ == "__main__":
+    pass
\ No newline at end of file
diff --git a/astrosource/astrosource.py b/astrosource/astrosource.py
index e4b0a83..80c97fe 100644
--- a/astrosource/astrosource.py
+++ b/astrosource/astrosource.py
@@ -4,6 +4,7 @@
 import ssl
 import sys
 from pathlib import Path
+import csv
 
 if (not os.environ.get('PYTHONHTTPSVERIFY', '') and getattr(ssl, '_create_unverified_context', None)):
     ssl._create_default_https_context = ssl._create_unverified_context
@@ -240,6 +241,22 @@ def plot(self, detrend=False, eebls=False, period=False, phaseShift=0.0, filesav
                 phased_plots(filterCode=self.filtercode, paths=self.paths, targets=self.targets, period=self.period, phaseShift=phaseShift)
         if eebls:
             plot_bls(paths=self.paths)
+        
+        # Output the list of targets that was used in this run.
+        #results_dir = parentPath / 'results'
+        #results_dir.mkdir(exist_ok=True)
+        
+        csv_path = self.paths['results'] / 'TargetsUsed.csv'
+        with open(csv_path, 'w', newline='') as csvfile:
+            writer = csv.writer(csvfile)
+            num_cols = len(self.targets[0]) + 1
+            header = [f'col{i}' for i in range(1, num_cols + 1)]
+            writer.writerow(header)
+        
+            for idx, targ in enumerate(self.targets, start=1):
+                writer.writerow([f'V{idx}'] + list(targ))
+        
+        
 
     def output(self, mode, data):
         output_files(paths=self.paths, photometrydata=data, mode=mode)
diff --git a/astrosource/comparison.py b/astrosource/comparison.py
index 3684601..4b1d7f6 100644
--- a/astrosource/comparison.py
+++ b/astrosource/comparison.py
@@ -13,7 +13,7 @@
 import matplotlib
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
-
+from multiprocessing import Pool, cpu_count
 from numpy import min, max, median, std, isnan, delete, genfromtxt, savetxt, load, \
     asarray, add, append, log10, average, array, where, cos
 from astropy.units import degree, arcsecond
@@ -22,19 +22,17 @@
 from astroquery.vizier import Vizier
 from tqdm import tqdm
 from prettytable import PrettyTable
-
+import traceback
 import requests
 import http
 import urllib3
-
+from collections import OrderedDict
 from astrosource.utils import AstrosourceException
 
 import logging
 
 logger = logging.getLogger('astrosource')
 
-
-
 def check_comparisons_files(parentPath=None, fileList=None, photFileArray=None, photSkyCoord=None, matchRadius=1.45):
     '''
     Find stable comparison stars for the target photometry
@@ -54,7 +52,6 @@ def check_comparisons_files(parentPath=None, fileList=None, photFileArray=None,
     if fileList==None:
         fileList= genfromtxt(parentPath / 'results/usedImages.txt', dtype=str)
 
-    #print (fileList)
     usedImages=[]
     q=0
     imageRemove=[]
@@ -73,7 +70,6 @@ def check_comparisons_files(parentPath=None, fileList=None, photFileArray=None,
             logger.debug('The Following file does not contain one or more of the provided comparison stars. It has been removed')
             logger.debug(file)
             logger.debug('**********************')
-            #fileList.remove(file)
             imageRemove.append(q)
         else:
             usedImages.append(file)
@@ -282,77 +278,88 @@ def read_data_files(parentPath, fileList):
 
 def ensemble_comparisons(photFileArray, compFile, parentPath, photSkyCoord):
 
-    # get rid of dumb for loop
-    fileCount = []
-    q=0
-    for photFile in photFileArray:
-        allCounts = 0.0
+    fileCount = np.zeros(len(photFileArray))  # Pre-allocate output array
 
-        if compFile.size ==2 and compFile.shape[0]==2:
-            matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree)
-        else:
-            matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree)
-
-        idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q])
-
-        allCounts = add(allCounts,sum(photFile[idx,4]))
+    # Prepare the `matchCoord` SkyCoord object
+    if compFile.size == 2 and compFile.shape[0] == 2:
+        matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree)
+    else:
+        matchCoord = SkyCoord(ra=compFile[:,0] * degree, dec=compFile[:,1] * degree)
 
-        fileCount.append(allCounts)
-        q=q+1
+    # Process all photSkyCoord entries in a vectorized way
+    for q, (photFile, skyCoord) in enumerate(zip(photFileArray, photSkyCoord)):
+        # Match coordinates
+        idx, _, _ = matchCoord.match_to_catalog_sky(skyCoord)
+        # Accumulate counts
+        fileCount[q] = np.sum(photFile[idx, 4])
 
-    logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(np.sum(np.array(fileCount))))
+    total_counts = np.sum(fileCount)
+    logger.debug("Total Ensemble Star Counts in Reference Frame {}".format(total_counts))
     return fileCount
 
+def process_phot_file_for_variation(args):
+    #q, matchCoord, photSkyCoord, photFileArray, fileCount = args
+    q, matchCoord, photSkyCoord, photFile, fileCount = args
+    idx, _, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q])
+    #photFile = photFileArray[q]
+    compDiffMags = 2.5 * np.log10(photFile[idx, 4] / fileCount[q])
+    instrMags = -2.5 * np.log10(photFile[idx, 4])
+    return compDiffMags, instrMags
+
 def calculate_comparison_variation(compFile, photFileArray, fileCount, parentPath, photSkyCoord):
-    stdCompStar=[]
-    sortStars=[]
+        
+    stdCompStar = []
+    sortStars = []
     logger.debug("Calculating Variation in Individual Comparisons")
 
-    if compFile.size ==2 and compFile.shape[0]==2:
-        compDiffMags = []
-
-        matchCoord = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree)
-        for q, photFile in enumerate(photFileArray):
-            idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q])
-            compDiffMags = append(compDiffMags,2.5 * log10(photFile[idx][4]/fileCount[q]))
-            instrMags = -2.5 * log10(photFile[idx][4])
-
+    if compFile.size == 2 and compFile.shape[0] == 2:
+        # Case for single comparison star
+        matchCoord = SkyCoord(ra=compFile[0] * degree, dec=compFile[1] * degree)
 
-        stdCompDiffMags=std(compDiffMags)
-        medCompDiffMags=np.nanmedian(compDiffMags)
-        medInstrMags= np.nanmedian(instrMags)
+        args = [(q, matchCoord, photSkyCoord, photFileArray[q], fileCount) for q in range(len(photFileArray))]
+        
+        with Pool(processes=max([cpu_count()-1,1])) as pool:
+            results = pool.map(process_phot_file_for_variation, args)
+        
+        compDiffMags = np.array([result[0] for result in results])
+        instrMags = np.array([result[1] for result in results])
+        
+        stdCompDiffMags = np.std(compDiffMags)
+        medCompDiffMags = np.nanmedian(compDiffMags)
+        medInstrMags = np.nanmedian(instrMags)
 
-        if np.isnan(stdCompDiffMags) :
+        if np.isnan(stdCompDiffMags):
             logger.error("Star Variability non rejected")
-            stdCompDiffMags=99
-        stdCompStar.append(stdCompDiffMags)
-        sortStars.append([compFile[0],compFile[1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])
-
+            stdCompDiffMags = 99
 
+        stdCompStar.append(stdCompDiffMags)
+        sortStars.append([compFile[0], compFile[1], stdCompDiffMags, medCompDiffMags, medInstrMags] + [0.0] * 8)
 
     else:
-        sortStars=[]
+        # Case for multiple comparison stars
+        matchCoord = SkyCoord(ra=compFile[:, 0] * degree, dec=compFile[:, 1] * degree)
         compDiffMags = []
-        instrMags=[]
-        matchCoord = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree)
+        instrMags = []
+
         for q, photFile in enumerate(photFileArray):
-            idx, d2d, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q])
-            compDiffMags.append(2.5 * log10(photFile[idx,4]/fileCount[q]))
-            instrMags.append(-2.5 * log10(photFile[idx,4]))
-
-        compDiffMags=array(compDiffMags)
-        instrMags=array(instrMags)
-
-        sortStars=[]
-        for z in range(len(compDiffMags[0])):
-            stdCompDiffMags=std(compDiffMags[:,z])
-            medCompDiffMags=np.nanmedian(compDiffMags[:,z])
-            medInstrMags=np.nanmedian(instrMags[:,z])
-            if np.isnan(stdCompDiffMags) :
+            idx, _, _ = matchCoord.match_to_catalog_sky(photSkyCoord[q])
+            compDiffMags.append(2.5 * np.log10(photFile[idx, 4] / fileCount[q]))
+            instrMags.append(-2.5 * np.log10(photFile[idx, 4]))
+
+        compDiffMags = np.array(compDiffMags)
+        instrMags = np.array(instrMags)
+
+        for z in range(compDiffMags.shape[1]):
+            stdCompDiffMags = np.std(compDiffMags[:, z])
+            medCompDiffMags = np.nanmedian(compDiffMags[:, z])
+            medInstrMags = np.nanmedian(instrMags[:, z])
+
+            if np.isnan(stdCompDiffMags):
                 logger.error("Star Variability non rejected")
-                stdCompDiffMags=99
+                stdCompDiffMags = 99
+
             stdCompStar.append(stdCompDiffMags)
-            sortStars.append([compFile[z,0],compFile[z,1],stdCompDiffMags,medCompDiffMags,medInstrMags,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])
+            sortStars.append([compFile[z, 0], compFile[z, 1], stdCompDiffMags, medCompDiffMags, medInstrMags] + [0.0] * 8)
 
     return stdCompStar, sortStars
 
@@ -361,10 +368,10 @@ def remove_stars_targets(parentPath, compFile, acceptDistance, targetFile, remov
     logger.info("Removing Target Stars from potential Comparisons")
     tableFound=False
 
-    if not (compFile.shape[0] == 2 and compFile.size ==2):
-        fileRaDec = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree)
-    else:
-        fileRaDec = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree)
+    # if not (compFile.shape[0] == 2 and compFile.size ==2):
+    #     fileRaDec = SkyCoord(ra=compFile[:,0]*degree, dec=compFile[:,1]*degree)
+    # else:
+    #     fileRaDec = SkyCoord(ra=compFile[0]*degree, dec=compFile[1]*degree)
 
     # Remove any nan rows from targetFile
     if targetFile is not None:
@@ -413,20 +420,13 @@ def remove_stars_targets(parentPath, compFile, acceptDistance, targetFile, remov
     tableFound=False
     # Check VSX for any known variable stars and remove them from the list
     logger.info("Searching for known variable stars in VSX......")
-    vServers = ['vizier.cfa.harvard.edu','vizier.u-strasbg.fr']
+    vServers = ['vizier.cds.unistra.fr','vizier.cfa.harvard.edu','vizier.u-strasbg.fr']
     vS=0
     try:
         while vS < len(vServers):
             v=Vizier(columns=['RAJ2000', 'DEJ2000']) # Skymapper by default does not report the error columns
             v.ROW_LIMIT=-1
-            
-                 #'vizier.nao.ac.jp',
-                 #]
-                 #'vizier.iucaa.in',
-                 #'vizier.china-vo.org',
-                 #'vizier.inasan.ru',
-                 #'vizier.idia.ac.za']
-            
+                        
             v.VIZIER_SERVER=vServers[vS]
             vS=vS + 1
             variableResult=v.query_region(avgCoord, radius=1.5*radius*degree, catalog='VSX')
@@ -517,16 +517,16 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd):
     TABLES = {'APASS':'II/336/apass9',
               'SDSS' :'V/154/sdss16',
               'PanSTARRS' : 'II/349/ps1',
-              'SkyMapper' : 'II/358/smss'
+              'SkyMapper' : 'II/379/smssdr4'
               }
 
     tbname = TABLES.get(cat_name, None)
     kwargs = {'radius': str(1.5*radius)+' deg'}
-    kwargs['catalog'] = cat_name
+    kwargs['catalog'] = cat_name 
 
-    vServers = ['vizier.u-strasbg.fr',
+    vServers = ['vizier.cds.unistra.fr','vizier.u-strasbg.fr']#,
          #'vizier.nao.ac.jp',
-         'vizier.cfa.harvard.edu']
+         #'vizier.cfa.harvard.edu']
          #'vizier.iucaa.in',
          #'vizier.china-vo.org',
          #'vizier.inasan.ru',
@@ -546,17 +546,15 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd):
         radecname = {'ra' :'raj2000', 'dec': 'dej2000'}
 
     if cat_name in ['APASS']:
-        searchColumns=[]
-        #searchColumns=[radecname['ra'], radecname['dec'], 'i_mag', opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')]
-        #print (searchColumns)
+
+        searchColumns=[radecname['ra'], radecname['dec'], opt['filter'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['error'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag'), opt['colmatch'].replace('i_mag','i\'mag').replace('r_mag','r\'mag').replace('g_mag','g\'mag'), opt['colerr'].replace('e_i_mag','e_i\'mag').replace('e_r_mag','e_r\'mag').replace('e_g_mag','e_g\'mag')]
+    
     else:
         searchColumns=[radecname['ra'], radecname['dec'], opt['filter'], opt['error'], opt['colmatch'], opt['colerr']]
-    #print (searchColumns)
-
+    
     v=Vizier(columns=searchColumns) # Skymapper by default does not report the error columns
     v.ROW_LIMIT=-1
 
-    #queryConstraint=
     # Filter out bad data from catalogues
     if cat_name == 'PanSTARRS':
         queryConstraint={'Qual' : '52 || 60 || 61'}
@@ -567,40 +565,28 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd):
     elif cat_name == 'APASS':
         queryConstraint={}
 
-
-
-    # vServers = ['vizier.cfa.harvard.edu','vizier.u-strasbg.fr']
-    # vS=0
-    # try:
-    #     while vS < len(vServers):
-    #         v=Vizier(columns=['RAJ2000', 'DEJ2000']) # Skymapper by default does not report the error columns
-    #         v.ROW_LIMIT=-1
-            
-    #              #'vizier.nao.ac.jp',
-    #              #]
-    #              #'vizier.iucaa.in',
-    #              #'vizier.china-vo.org',
-    #              #'vizier.inasan.ru',
-    #              #'vizier.idia.ac.za']
-            
-    #         v.VIZIER_SERVER=vServers[vS]
-    #         vS=vS + 1
-    #         variableResult=v.query_region(avgCoord, radius=1.5*radius*degree, catalog='VSX')
-    #         if str(variableResult)=="Empty TableList":
-    #             logger.info("VSX Returned an Empty Table from " + str(vServers[vS]) + ".")
-    #         else:
-    #             variableResult=variableResult['B/vsx/vsx']
-    #             tableFound=True
-
     cycler=0
+    vS=0
     try:
+        vtimeout=180
         while vS < len(vServers):
+            
             v.VIZIER_SERVER=vServers[vS]
-            query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs)
-            if str(query)=="Empty TableList":
+            v.TIMEOUT = vtimeout 
+            try:
+                query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs)
+                if str(query)=="Empty TableList":
+                    vS=vS+1
+                else:
+                    break
+            except requests.exceptions.ReadTimeout:
+                
+                logger.info("Read Timeout for Vizier, might be a very dense field or slow server today")
                 vS=vS+1
-            else:
-                break
+            except:
+                logger.info("Failed vizier query")
+                print(traceback.print_exc())
+            
             
     except VOSError:
         raise AstrosourceException("Could not find RA {} Dec {} in {}".format(avgCoord.ra.value,avgCoord.dec.value, cat_name))
@@ -655,13 +641,6 @@ def catalogue_call(avgCoord, radius, opt, cat_name, targets, closerejectd):
     if len(resp) != 0:
         catReject=[]
         
-        # v=Vizier() # Skymapper by default does not report the error columns
-        # v.ROW_LIMIT=-1
-        # query = v.query_region(avgCoord, column_filters=queryConstraint, **kwargs)
-        # resp = query[tbname]
-
-        #breakpoint()
-
         for q in range(len(resp)):
             if np.asarray(resp[opt['filter']][q]) == 0.0 or np.asarray(resp[opt['error']][q]) == 0.0 or np.isnan(resp[opt['filter']][q]) or np.isnan(resp[opt['error']][q]):
                 catReject.append(q)
@@ -741,8 +720,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
     if not calibPath.exists():
         os.makedirs(calibPath)
 
-    #Vizier.ROW_LIMIT = -1
-
     # Get List of Files Used
     fileList=[]
     for line in (parentPath / "results/usedImages.txt").read_text().strip().split('\n'):
@@ -804,7 +781,15 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         raise AstrosourceException(f"{filterCode} is not accepted at present")
 
     # Look up in online catalogues and make sure there are sufficient comparison stars
-
+    # If in the realm of skymapper and a skymapper filter, put skymapper at the front
+    if avgCoord.dec.deg < 15 and 'SkyMapper' in catalogues:
+        ordered = OrderedDict()
+        ordered['SkyMapper'] = catalogues['SkyMapper']
+        for key, val in catalogues.items():
+            if key != 'SkyMapper':
+                ordered[key] = val
+        catalogues = ordered
+        
     coords=[]
     for cat_name, opt in catalogues.items():
         try:
@@ -817,8 +802,18 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
                 elif cat_name == 'SkyMapper' and noskymapper==True:
                     logger.info("Skipping Skymapper")
                 else:
-
-                    coords = catalogue_call(avgCoord, 1.5*radius, opt, cat_name, targets=targets, closerejectd=closerejectd)
+                    try:
+                        coords = catalogue_call(avgCoord, 1.5*radius, opt, cat_name, targets=targets, closerejectd=closerejectd)
+                    except Exception as e:
+                        
+                        if 'Could not find RA' in str(e):
+                            raise AstrosourceException("Empty catalogue produced from catalogue call")
+                            #pass
+                        else:
+                            print (e)
+                            print ("MTF hunting this bug")
+                            print(traceback.print_exc())
+                            breakpoint()
                     # If no results try next catalogue
                     if len(coords.ra) == 0:
                         coords=[]
@@ -829,111 +824,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
                     else:
                         max_sep=1.5 * arcsecond
 
-                    # The following commented out code is temporary code to manually access Skymapper DR3 and APASS10 prior to them becoming publically available through vizier.
-                    # It will be removed in later versions.
-
-                    # SKYMAPPER OVERRIDE
-                    # max_sep=1.5 * arcsecond
-                    # smapdr3=genfromtxt(parentPath / 'SkymapperDR3.csv', dtype=float, delimiter=',')
-                    # coords = namedtuple(typename='data',field_names=['ra','dec','mag','emag','cat_name', 'colmatch', 'colerr'])
-                    # coords.ra=smapdr3[:,0]
-                    # coords.dec=smapdr3[:,1]
-                    # print (opt['filter'])
-                    # if opt['filter'] == 'uPSF':
-                    #     coords.mag = smapdr3[:,2]
-                    #     coords.emag = smapdr3[:,3]
-                    #     coords.colmatch = smapdr3[:,4]
-                    #     coords.colerr = smapdr3[:,5]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'gPSF':
-                    #     coords.mag = smapdr3[:,4]
-                    #     coords.emag = smapdr3[:,5]
-                    #     coords.colmatch = smapdr3[:,6]
-                    #     coords.colerr = smapdr3[:,7]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'rPSF':
-                    #     coords.mag = smapdr3[:,6]
-                    #     coords.emag = smapdr3[:,7]
-                    #     coords.colmatch = smapdr3[:,8]
-                    #     coords.colerr = smapdr3[:,9]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'iPSF':
-                    #     coords.mag = smapdr3[:,8]
-                    #     coords.emag = smapdr3[:,9]
-                    #     coords.colmatch = smapdr3[:,6]
-                    #     coords.colerr = smapdr3[:,7]
-                    #     coords.colrev = 1
-
-                    # if opt['filter'] == 'zPSF':
-                    #     coords.mag = smapdr3[:,10]
-                    #     coords.emag = smapdr3[:,11]
-                    #     coords.colmatch = smapdr3[:,8]
-                    #     coords.colerr = smapdr3[:,9]
-                    #     coords.colrev = 1
-
-                    # APASS OVERRIDE
-                    # max_sep=2.5 * arcsecond
-                    # smapdr3=genfromtxt(parentPath / 'ap10.csv', dtype=float, delimiter=',')
-                    # coords = namedtuple(typename='data',field_names=['ra','dec','mag','emag','cat_name', 'colmatch', 'colerr'])
-                    # coords.ra=smapdr3[:,0]
-                    # coords.dec=smapdr3[:,1]
-                    # print (opt['filter'])
-
-                    # if opt['filter'] == 'Bmag':
-                    #     coords.mag = smapdr3[:,2]
-                    #     coords.emag = smapdr3[:,3]
-                    #     coords.colmatch = smapdr3[:,4]
-                    #     coords.colerr = smapdr3[:,5]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'Vmag':
-                    #     coords.mag = smapdr3[:,4]
-                    #     coords.emag = smapdr3[:,5]
-                    #     coords.colmatch = smapdr3[:,2]
-                    #     coords.colerr = smapdr3[:,3]
-                    #     coords.colrev = 1
-
-                    # if opt['filter'] == 'umag':
-                    #     coords.mag = smapdr3[:,6]
-                    #     coords.emag = smapdr3[:,7]
-                    #     coords.colmatch = smapdr3[:,8]
-                    #     coords.colerr = smapdr3[:,9]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'gmag':
-                    #     coords.mag = smapdr3[:,8]
-                    #     coords.emag = smapdr3[:,9]
-                    #     coords.colmatch = smapdr3[:,10]
-                    #     coords.colerr = smapdr3[:,11]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'rmag':
-                    #     coords.mag = smapdr3[:,10]
-                    #     coords.emag = smapdr3[:,11]
-                    #     coords.colmatch = smapdr3[:,12]
-                    #     coords.colerr = smapdr3[:,13]
-                    #     coords.colrev = 0
-
-                    # if opt['filter'] == 'imag':
-                    #     coords.mag = smapdr3[:,12]
-                    #     coords.emag = smapdr3[:,13]
-                    #     coords.colmatch = smapdr3[:,10]
-                    #     coords.colerr = smapdr3[:,11]
-                    #     coords.colrev = 1
-
-                    # if opt['filter'] == 'zmag':
-                    #     coords.mag = smapdr3[:,14]
-                    #     coords.emag = smapdr3[:,15]
-                    #     coords.colmatch = smapdr3[:,12]
-                    #     coords.colerr = smapdr3[:,13]
-                    #     coords.colrev = 1
-
                     # Save catalogue search
                     catalogueOut=np.transpose(np.vstack(np.array([[coords.ra],[coords.dec],[coords.mag],[coords.emag],[coords.colmatch],[coords.colerr]])))
-                    #catalogueOut=np.array([[coords.ra],[coords.dec],[coords.mag],[coords.emag],[coords.colmatch],[coords.colerr]])
                     
                     savetxt(parentPath / "results/catalogueSearch.csv", np.asarray(catalogueOut), newline='\n' , delimiter=",", fmt='%0.8f',header='RA,Dec,Mag,err_Mag,colourmatch_mag,colourmatch_mag_err')
                     del catalogueOut
@@ -990,13 +882,11 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
 
 
                     ### If looking for colour, remove those without matching colour information
-
                     calibStandsReject=[]
 
                     if (asarray(calibStands).shape[0] != 9 and asarray(calibStands).size !=9) and len(calibStands) != 0:
                         for q in range(len(asarray(calibStands)[:,0])):
                             reject=0
-                            #if colourdetect == True:
                             # There needs to be a colour in here to make various plots.
                             if np.isnan(calibStands[q][6]): # if no matching colour
                                 reject=1
@@ -1025,10 +915,14 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
                             logger.debug('Calibration Stars Identified below')
                             tabl = PrettyTable()
                             tabl.field_names = ["RA","Dec","MagDiff","Mag"]
-                            tabl.add_rows(calibStands[:,0:4])
+                            #abl.add_rows(calibStands[:,0:4])
+                            rows = calibStands[:,0:4].tolist()
+                            if rows:                 
+                                tabl.add_rows(rows)
                             logger.debug(tabl)
                     except:
                         print ("something fishy in the calib mags table")
+                        logger.error(traceback.print_exc())
 
                     # Get the set of least variable stars to use as a comparison to calibrate the files (to eventually get the *ACTUAL* standards
                     if asarray(calibStands).shape[0] == 0:
@@ -1065,8 +959,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
 
 
     #Colour Term Estimation routine
-    
-
     # use a temporary calibStands array for colour terms
     arrayCalibStands=np.asarray(calibStands)
 
@@ -1110,7 +1002,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         colTemp=delete(colTemp, calibStandsReject, axis=0)
 
     # Pre-colour Reference plot
-    
+    try:
         plt.cla()
         fig = plt.gcf()
         outplotx=colTemp[:,4]
@@ -1127,9 +1019,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
             x, residuals, rank, s = sqsol
             plt.plot(outplotx,m*outplotx+c,'r')
         except:
-            logger.info("Couldn't fit a line to pre-colour reference plot")
-
-        
+            logger.info("Couldn't fit a line to pre-colour reference plot")        
         
         plt.ylim(max(outploty)+0.05,min(outploty)-0.05)
         plt.xlim(min(outplotx)-0.05,max(outplotx)+0.05)
@@ -1141,9 +1031,11 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         plt.savefig(parentPath / str("results/CalibrationSanityPlot_PreColour_Reference.eps"))
     
         logger.info('Estimated Colour Slope in Reference Frame: ' + str(m))
-    
+    except:
+        print ("Couldn't make colour plot")
+        logger.error(traceback.print_exc())
+        
     # MAKE ALL PRE-COLOUR PLOTS
-
     if colourdetect == True and colourTerm == 0.0:
 
         logger.info('Estimating Colour Slope from all Frames...')
@@ -1161,14 +1053,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
             os.makedirs(colourPath)
         for qert in range(len(photFileHolder)):
 
-            
-            #breakpoint()
-            #photFrame = load(parentPath / filen)
             photFrame = photFileHolder[z]
             z=z+1
-            #photFrame[:,5] = 1.0857 * (photFrame[:,5]/photFrame[:,4])
-            #photFrame[:,4]=-2.5 * np.log10(photFrame[:,4])
-
             photCoords=SkyCoord(ra=photFrame[:,0]*degree, dec=photFrame[:,1]*degree)
             colTemp=[]
             for q in range(len(arrayCalibStands[:,0])):
@@ -1344,14 +1230,9 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         # Get colour information into photFile
         # adding in colour columns to photfile
         
-        idx,d2d,_=catCoords.match_to_catalog_sky(photCoords)
-        
-        
-
+        idx,d2d,_=catCoords.match_to_catalog_sky(photCoords)      
         rejectStars=where(d2d > max_sep)[0]
-        idx=delete(idx, rejectStars, axis=0)
-
-        
+        idx=delete(idx, rejectStars, axis=0)       
 
         for q in range(len(idx)):
             if colrev == 1:
@@ -1361,16 +1242,10 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
             photFile[idx[q],9]=pow(pow(coords.colerr[q],2)+pow(coords.emag[q],2),0.5)
             photFile[idx[q],10]=1
 
-        #breakpoint()
-
         #Replace undetected colours with average colour of the field
         photFile[:,8]=np.nan_to_num(photFile[:,8],nan=np.nanmedian(photFile[:,8]))
         photFile[:,9]=np.nan_to_num(photFile[:,9],nan=np.nanmedian(photFile[:,9]))
-        # If the colour detection routine was run, this was calculated already
-        #if not colourdetect:
-        #    photFile[:,5]=1.0857 * (photFile[:,5]/photFile[:,4])
-
-
+        
         notNanTemp=photFile[np.isnan(photFile).any(axis=1)] #rows where any value is nan
         if colourdetect:
             notNanTemp[:,4]=notNanTemp[:,4] -(colourTerm*notNanTemp[:,8])
@@ -1379,9 +1254,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
 
         emptyCalibFlag=0
         NanTemp=photFile[~np.isnan(photFile).any(axis=1)] #rows where no value is nan
-        if len(NanTemp) != 0:
-            #if not colourdetect:
-            #    NanTemp[:,4]=-2.5*log10(photFile[:,4])
+        if len(NanTemp) != 0:            
             NanTemp[:,10]=2 # 2 means that there was no colour to use to embed the colour.
         else:
             emptyCalibFlag=1
@@ -1398,9 +1271,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         else:
             calibCoord=SkyCoord(ra=calibStands[:,0]*degree,dec=calibStands[:,1]*degree)
 
-
-
-
         if calibStands.size == 13 and calibStands.shape[0]== 13:
             for q in range(len(calibStands[:,0])):
                 idx,d2d,_=calibCoord.match_to_catalog_sky(photCoords)
@@ -1450,10 +1320,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
             weights=1/(calibOut[:,1])
             linA = np.vstack([outplotx,np.ones(len(outplotx))]).T * np.sqrt(weights[:,np.newaxis])
             linB = outploty * np.sqrt(weights)
-
-            #breakpoint()
-            
-
             
             try:
                 sqsol = np.linalg.lstsq(linA,linB, rcond=None)
@@ -1471,8 +1337,7 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
                 plt.subplots_adjust(left=0.15, right=0.98, top=0.98, bottom=0.17, wspace=0.3, hspace=0.4)
                 fig.set_size_inches(6,3)
                 plt.savefig(colourPath / str("CalibrationSanityPlot_Colour_" + str(z) + "_Post.png"))
-                plt.savefig(colourPath  / str("CalibrationSanityPlot_Colour_" + str(z) + "_Post.eps"))
-                
+                plt.savefig(colourPath  / str("CalibrationSanityPlot_Colour_" + str(z) + "_Post.eps"))                
             except:
                 logger.info("\nFailed linear fit on a postcolour plot: ")
 
@@ -1497,7 +1362,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
 
     # PRINT FINAL REFERENCE IMAGE COLOUR PLOT
     
-    #referenceFrame = genfromtxt(parentPath / 'results/referenceFrame.csv', dtype=float, delimiter=',')
     POSTreferenceFrame=np.copy(referenceFrame)
     # add the colour error in here
     POSTreferenceFrame[:,5] = pow(pow(referenceFrame[:,5],2) + pow(colourError,2),0.5)
@@ -1630,7 +1494,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         fig = plt.gcf()
         outplotx=calibOverlord[:,0]
         outploty=calibOverlord[:,5]
-        #breakpoint()
         sqsol = np.linalg.lstsq(np.vstack([calibOverlord[:,0],np.ones(len(calibOverlord[:,0]))]).T,calibOverlord[:,5], rcond=None)
         m, c = sqsol[0]
         x, residuals, rank, s = sqsol
@@ -1668,7 +1531,6 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
     except:
         logger.info("Could not create Difference versus Magnitude calibration plot")
 
-
     # Non-linearity correction routine
     if linearise == False:
         logger.info("Skipping Non-linearity correction routine")
@@ -1741,14 +1603,10 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
         logger.debug("CORRECTING EACH FILE FOR NONLINEARITY")
         correctFileList = calibPath.glob("*.calibrated.csv")
         for file in correctFileList:
-
             # NEED TO FIX UP COMPARED!
-
             photFile = genfromtxt(file, dtype=float, delimiter=',')
-
             for q in range(len(photFile[:,0])):
                 photFile[q,4]=photFile[q,4] - nonlinearSlope*photFile[q,4]+nonlinearZero
-
             savetxt(file, photFile, delimiter=",", fmt='%0.8f')
 
     # Difference vs time calibration plot
@@ -1820,4 +1678,8 @@ def find_comparisons_calibrated(targets, paths, filterCode, nopanstarrs=False, n
     savetxt(parentPath / "results/calibCompsUsed.csv", compFile, delimiter=",", fmt='%0.8f',header='RA,Dec,Variability')
     sys.stdout.write('\n')
 
-    return colourTerm, colourError, compFile
\ No newline at end of file
+    return colourTerm, colourError, compFile
+
+# Needed for windows to multiprocess appropriately
+if __name__ == "__main__":
+    pass
diff --git a/astrosource/detrend.py b/astrosource/detrend.py
index 0407aef..dc6acda 100644
--- a/astrosource/detrend.py
+++ b/astrosource/detrend.py
@@ -29,74 +29,35 @@ def detrend_data(paths, filterCode, detrendfraction=0.1):
 
     fileList = paths['outcatPath'].glob('*diffExcel*csv')
     r=0
-    #logger.debug(fileList)
     for file in fileList:
-        #photFile = load(paths['parent'] / file)
-        #print (file)
         photFile = genfromtxt(file, dtype=float, delimiter=',')
-        #print (photFile)
         exists=os.path.isfile(str(file).replace('diff','calib'))
         if exists:
-            calibFile = genfromtxt(str(file).replace('diff','calib'), dtype=float, delimiter=',')
-            #logger.debug("Calibration difference")
-            #logger.debug(-(photFile[:,1]-calibFile[:,1])[0])
+            calibFile = genfromtxt(str(file).replace('diff','calib'), dtype=float, delimiter=',')            
             calibDiff=-((photFile[:,1]-calibFile[:,1])[0])
-        #logger.debug(photFile[:,1])
-        #logger.debug(photFile[:,0])
-        
-        #print (calibFile)
-        #logger.debug(file)
-        #logger.debug(photFile[:,1])
-
+       
         baseSubDate=min(photFile[:,0])
-        #logger.debug(baseSubDate)
-        #logger.debug(math.floor(baseSubDate))
-
-        photFile[:,0]=photFile[:,0]-baseSubDate
-
-
-        
-
-        #leftMost = click.prompt("Enter left side most valid date:")
-        #leftFlat = click.prompt("Enter left side end of flat region:")
-
-        #rightFlat = click.prompt("Enter right side start of flat region:")
-        #rightMost = click.prompt("Enter right side most valid date:")
-
-        
+        photFile[:,0]=photFile[:,0]-baseSubDate        
         detrendPercent=0.1
         leftMost=min(photFile[:,0])
-        leftFlat=min(photFile[:,0])+ (detrendPercent * max(photFile[:,0])-min(photFile[:,0]))
-        
+        leftFlat=min(photFile[:,0])+ (detrendPercent * max(photFile[:,0])-min(photFile[:,0]))        
         rightFlat=max(photFile[:,0])- (detrendPercent * max(photFile[:,0])-min(photFile[:,0]))
         rightMost=max(photFile[:,0])
-
-
-        #print (leftMost)
-        #print (leftFlat)
-        #print (rightFlat)
-        #print (rightMost)
-
         # Clip off edges
         clipReject=[]
         for i in range(photFile.shape[0]):
             if photFile[i,0] < float(leftMost) or photFile[i,0] > float(rightMost):
                 clipReject.append(i)
-                #logger.debug(photFile[i,1])
-                #logger.debug("REJECT")
         logger.debug(photFile.shape[0])
         photFile=delete(photFile, clipReject, axis=0)
         logger.debug(photFile.shape[0])
 
-
         # Get an array holding only the flat bits
         transitReject=[]
         flatFile=asarray(photFile)
         for i in range(flatFile.shape[0]):
             if not (flatFile[i,0] > float(leftMost) and flatFile[i,0] < float(leftFlat)) or (flatFile[i,0] > float(rightFlat) and flatFile[i,0] < float(rightMost)):
                 transitReject.append(i)
-                #logger.debug(flatFile[i,0])
-                #logger.debug("REJECT")
         logger.debug(flatFile.shape[0])
         flatFile=delete(flatFile, transitReject, axis=0)
         logger.debug(flatFile.shape[0])
@@ -113,7 +74,6 @@ def detrend_data(paths, filterCode, detrendfraction=0.1):
             for i in range(flatFile.shape[0]):
                 flatFile[i,1]=flatFile[i,1]-(polyFit[1]+(polyFit[0]*flatFile[i,0]))
 
-
         #Remove trend from actual data
         if polyFitRequest == 2:
             for i in range(photFile.shape[0]):
@@ -122,7 +82,6 @@ def detrend_data(paths, filterCode, detrendfraction=0.1):
             for i in range(photFile.shape[0]):
                 photFile[i,1]=photFile[i,1]-(polyFit[1]+(polyFit[0]*photFile[i,0]))
 
-
         #return basedate to the data
         photFile[:,0]=photFile[:,0]+baseSubDate
 
diff --git a/astrosource/identify.py b/astrosource/identify.py
index bfc65d1..dcd0571 100644
--- a/astrosource/identify.py
+++ b/astrosource/identify.py
@@ -4,8 +4,11 @@
 import os
 import logging
 import pickle
+from concurrent.futures import ProcessPoolExecutor
 
-from numpy import genfromtxt, delete, asarray, save, savetxt, load, transpose, isnan, zeros, max, min, nan, where, average, cos, hstack, array, column_stack, copy, c_, sqrt, ptp
+#import copy
+
+from numpy import genfromtxt, delete, argmax, asarray, save, savetxt, load, transpose, isnan, zeros, max, min, nan, where, average, cos, hstack, array, column_stack, copy, c_, sqrt, ptp
 from astropy import units as u
 from astropy.units import degree, arcsecond
 from astropy import wcs
@@ -15,13 +18,61 @@
 from barycorrpy import utc_tdb
 from tqdm import tqdm
 from prettytable import PrettyTable
-from multiprocessing import Pool
+from multiprocessing import Pool,cpu_count
 
+import platform
 from astrosource.utils import AstrosourceException
 from astrosource.comparison import catalogue_call
 
 logger = logging.getLogger('astrosource')
 
+def process_fits_file(f, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut):
+    try:
+        # Determine if the file is from S3
+        s3 = False
+        try:
+            fitsobj = Path(f)
+        except TypeError:
+            fitsobj = f.open()
+            s3 = True
+
+        # Extract photometry
+        filepath, photFile = extract_photometry(
+            fitsobj, indir, bjd=bjd, ignoreedgefraction=ignoreedgefraction, 
+            lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut
+        )
+
+        if filepath is None or photFile is None or len(photFile) == 0:
+            return None  # Skip if extraction failed
+
+        # Close S3 file if applicable
+        if s3:
+            f.close()
+            filename = f.name
+        else:
+            filename = fitsobj.name
+
+        # Return results if photFile is valid
+        if photFile.size > 100:
+            photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree)
+            return filepath, photFile, photSkyCoord
+        return None
+    except Exception as e:
+        print(f"Error processing file {f}: {e}")
+        return None
+
+def process_files_multiprocessing(filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut):
+    # Wrap arguments for multiprocessing
+    args = (indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut)
+    with Pool(processes=max([cpu_count()-1,1])) as pool:
+        results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(list(filelist))))
+
+    # Filter valid results
+    results = [res for res in results if res is not None]
+    phot_dict, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], [])
+    return list(phot_dict), list(photFileHolder), list(photSkyCoord)    
+
+
 def rename_data_file(prihdr, bjd=False):
 
     prihdrkeys = prihdr.keys()
@@ -47,8 +98,6 @@ def rename_data_file(prihdr, bjd=False):
     airMass=(str(prihdr['AIRMASS']).replace('.','a'))
     instruMe=(prihdr['INSTRUME']).replace(' ','').replace('/','').replace('-','')
 
-
-
     if (prihdr['MJD-OBS'] == 'UNKNOWN'):
         timeobs = 'UNKNOWN'
     elif bjd:
@@ -63,86 +112,11 @@ def export_photometry_files(filelist, indir, filetype='csv', bjd=False, ignoreed
     phot_dict = []
     new_files = []
     photFileHolder=[]
-    photSkyCoord=[]
-    # filelist = list(filelist)
-    # for f in tqdm(filelist):
-    #     s3 = False
-    #     try:
-    #         fitsobj = Path(f)
-    #     except TypeError:
-    #         fitsobj = f.open()
-    #         s3 = True
-
-
-    #     filepath, photFile = extract_photometry(fitsobj, indir, bjd=bjd, ignoreedgefraction=ignoreedgefraction, lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut)
-    #     if not filepath and not photFile :
-    #         continue
-    #     if s3:
-    #         f.close()
-    #         filename = f.name
-    #     else:
-    #         filename = fitsobj.name
-
-    #     if photFile.size > 100:
-    #         phot_dict.append(filepath)
-
-    #         photFileHolder.append(photFile)
-    #         photSkyCoord.append(SkyCoord(ra=photFile[:,0]*u.degree, dec=photFile[:,1]*u.degree))
-    
-    def process_fits_file(f, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut):
-        try:
-            # Determine if the file is from S3
-            s3 = False
-            try:
-                fitsobj = Path(f)
-            except TypeError:
-                fitsobj = f.open()
-                s3 = True
-    
-            # Extract photometry
-            filepath, photFile = extract_photometry(
-                fitsobj, indir, bjd=bjd, ignoreedgefraction=ignoreedgefraction, 
-                lowestcounts=lowestcounts, racut=racut, deccut=deccut, radiuscut=radiuscut
-            )
-    
-            if not filepath or not photFile:
-                return None  # Skip if extraction failed
-    
-            # Close S3 file if applicable
-            if s3:
-                f.close()
-                filename = f.name
-            else:
-                filename = fitsobj.name
-    
-            # Return results if photFile is valid
-            if photFile.size > 100:
-                photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree)
-                return filepath, photFile, photSkyCoord
-            return None
-        except Exception as e:
-            print(f"Error processing file {f}: {e}")
-            return None
-    
-    
-    def process_files_multiprocessing(filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut):
-        # Wrap arguments for multiprocessing
-        args = (indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut)
-        with Pool() as pool:
-            results = list(tqdm(pool.starmap(process_fits_file, [(f, *args) for f in filelist]), total=len(filelist)))
-    
-        # Filter valid results
-        results = [res for res in results if res is not None]
-        phot_dict, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], [])
-        return list(phot_dict), list(photFileHolder), list(photSkyCoord)
-    
-    
+    photSkyCoord=[]    
     # Example usage
     phot_dict, photFileHolder, photSkyCoord = process_files_multiprocessing(
             filelist, indir, bjd, ignoreedgefraction, lowestcounts, racut, deccut, radiuscut
-        )
-
-    
+        )    
 
     return phot_dict, photFileHolder, photSkyCoord
 
@@ -208,14 +182,20 @@ def extract_photometry(infile, parentPath, outfile=None, bjd=False, ignoreedgefr
                         photFile[:,0] = photFile[:,0] - 180.0
                         for entry in range(len(photFile[:,0])):
                             if photFile[entry,0] < 0:
-                                photFile[entry,0] = photFile[entry,0] + 360
+                                photFile[entry,0] = photFile[entry,0] + 360     
+                                
                 #remove odd zero entries
                 photFile[:,0][photFile[:,0] == 0.0 ] = nan
                 photFile[:,0][photFile[:,0] == 0.0 ] = nan
                 photFile[:,1][photFile[:,1] == 0.0 ] = nan
                 photFile[:,1][photFile[:,1] == 0.0 ] = nan
 
-                photFile=photFile[~isnan(photFile).any(axis=1)]
+                num_rows_with_nans = sum(isnan(photFile).any(axis=1))
+                
+                if len(photFile) == num_rows_with_nans:
+                    print(f"REJECT {infile}: Very likely no wcs fit - RA and Dec columns are zero.")
+
+                photFile=photFile[~isnan(photFile).any(axis=1)]                               
 
                 #remove lowcounts
                 rejectStars=where(photFile[:,4] < lowestcounts)[0]
@@ -229,183 +209,92 @@ def extract_photometry(infile, parentPath, outfile=None, bjd=False, ignoreedgefr
 
     return outfile, photFile
 
-def convert_photometry_files(filelist, ignoreedgefraction=0.05, lowestcounts=1800,  racut=-99.9, deccut=-99.9, radiuscut=-99.9):
 
-    new_files = []
-    photFileHolder=[]
-    photSkyCoord=[]
-    #breakpoint()
-    # for fn in filelist:
-    #     photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1)
-    #     logger.info(fn)
-    #     # reject nan entries in file
-    #     if photFile.size > 16 and photFile.shape[1] == 8: #ignore zero sized files and files with only one or two entries
-    #         if max(photFile[:,0]) < 360 and max(photFile[:,1]) < 90:
-
-    #             if (photFile.size > 50):
-    #                 if not (( asarray(photFile[:,0]) > 360).sum() > 0) :
-    #                     if not(( asarray(photFile[:,1]) > 90).sum() > 0) :
-
-    #                         photFile=photFile[~isnan(photFile).any(axis=1)]
-
-
-
-    #                         # if radial cut do that otherwise chop off image edges
-    #                         if racut != -99.9 and deccut !=-99.9 and radiuscut !=-99.9:
-    #                             distanceArray=pow((pow(photFile[:,0]-float(racut),2)+pow(photFile[:,1]-float(deccut),2)),0.5)
-    #                             photFile=delete(photFile, where(distanceArray > float(radiuscut)/60), axis=0)
-    #                         else:
-    #                             #Routine to deal with RA 0 crossovers
-    #                             crossover=0
-    #                             if (max(photFile[:,0])-min(photFile[:,0])) > 180:
-    #                                 for entry in range(len(photFile[:,0])):
-    #                                     if photFile[entry,0] > 180:
-    #                                         photFile[entry,0] = photFile[entry,0] - 180
-    #                                     else:
-    #                                         photFile[entry,0] = photFile[entry,0] + 180
-
-    #                                 crossover=1
-
-    #                             # Remove edge detections
-    #                             raRange=(max(photFile[:,0])-min(photFile[:,0]))
-    #                             decRange=(max(photFile[:,1])-min(photFile[:,1]))
-    #                             raMin=min(photFile[:,0])
-    #                             raMax=max(photFile[:,0])
-    #                             decMin=min(photFile[:,1])
-    #                             decMax=max(photFile[:,1])
-    #                             raClip=raRange*ignoreedgefraction
-    #                             decClip=decRange*ignoreedgefraction
-    #                             photFile[:,0][photFile[:,0] < raMin + raClip ] = nan
-    #                             photFile[:,0][photFile[:,0] > raMax - raClip ] = nan
-    #                             photFile[:,1][photFile[:,1] > decMax - decClip ] = nan
-    #                             photFile[:,1][photFile[:,1] < decMin + decClip ] = nan
-    #                             if crossover == 1:
-    #                                 photFile[:,0] = photFile[:,0] - 180.0
-    #                                 for entry in range(len(photFile[:,0])):
-    #                                     if photFile[entry,0] < 0:
-    #                                         photFile[entry,0] = photFile[entry,0] + 360
-    #                         #remove odd zero entries
-    #                         photFile[:,0][photFile[:,0] == 0.0 ] = nan
-    #                         photFile[:,0][photFile[:,0] == 0.0 ] = nan
-    #                         photFile[:,1][photFile[:,1] == 0.0 ] = nan
-    #                         photFile[:,1][photFile[:,1] == 0.0 ] = nan
-    #                         photFile=photFile[~isnan(photFile).any(axis=1)]
-
-    #                         #remove lowcounts
-    #                         rejectStars=where(photFile[:,4] < lowestcounts)[0]
-    #                         photFile=delete(photFile, rejectStars, axis=0)
-
-    #                         filepath = Path(fn).with_suffix('.npy')
-                            
-    #                         photFile=c_[photFile,zeros(len(photFile[:,0])),zeros(len(photFile[:,0])),zeros(len(photFile[:,0])),zeros(len(photFile[:,0]))]
-    #                         photFile[:,8]=nan
-    #                         photFile[:,9]=nan
-    #                         photFile[:,10]=0
-    #                         photFile[:,11]=copy(photFile[:,4])
-
-    #                         if photFile.size > 16:
-    #                             new_files.append(filepath.name)
-    #                             photFileHolder.append(photFile)
-    #                             photSkyCoord.append(SkyCoord(ra=photFile[:,0]*u.degree, dec=photFile[:,1]*u.degree))
-                                
-                            
-
-                            
-    #                     else:
-    #                         logger.debug("REJECT")
-    #                         logger.debug(fn)
-    #                 else:
-    #                     logger.debug("REJECT")
-    #                     logger.debug(fn)
-    #             else:
-    #                 logger.debug("REJECT")
-    #                 logger.debug(fn)
-    #         else:
-    #             logger.debug("REJECT")
-    #             logger.debug(fn)
-    #     else:
-    #         logger.debug("REJECT")
-    #         logger.debug(fn)
-    
-    
-    def process_file(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger):
-        try:
-            # Read the file
-            photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1)
-            logger.info(fn)
-    
-            # Validate file shape and range
-            if photFile.size <= 16 or photFile.shape[1] != 8:
-                logger.debug(f"REJECT {fn}: Invalid size or shape.")
-                return None
-    
-            if max(photFile[:, 0]) >= 360 or max(photFile[:, 1]) >= 90:
-                logger.debug(f"REJECT {fn}: Invalid RA/Dec range.")
-                return None
-    
-            # Reject NaN entries
-            photFile = photFile[~isnan(photFile).any(axis=1)]
-    
-            # Apply radial cut or image edge trimming
-            if racut != -99.9 and deccut != -99.9 and radiuscut != -99.9:
-                distanceArray = sqrt((photFile[:, 0] - racut)**2 + (photFile[:, 1] - deccut)**2)
-                photFile = delete(photFile, where(distanceArray > radiuscut / 60), axis=0)
-            else:
-                # Handle RA crossover
-                if (max(photFile[:, 0]) - min(photFile[:, 0])) > 180:
-                    photFile[:, 0] = where(photFile[:, 0] > 180, photFile[:, 0] - 360, photFile[:, 0] + 360)
-    
-                # Clip edges
-                raRange, decRange = ptp(photFile[:, 0]), ptp(photFile[:, 1])
-                raClip, decClip = raRange * ignoreedgefraction, decRange * ignoreedgefraction
-                photFile = photFile[
-                    (photFile[:, 0] > min(photFile[:, 0]) + raClip) &
-                    (photFile[:, 0] < max(photFile[:, 0]) - raClip) &
-                    (photFile[:, 1] > min(photFile[:, 1]) + decClip) &
-                    (photFile[:, 1] < max(photFile[:, 1]) - decClip)
-                ]
-    
-            # Remove zero and low-count entries
-            photFile = photFile[(photFile[:, 4] > lowestcounts) & (photFile[:, 0] != 0) & (photFile[:, 1] != 0)]
-    
-            # Prepare output
-            photFile = c_[photFile, zeros((len(photFile), 4))]
-            photFile[:, 8:10] = nan
-            photFile[:, 11] = photFile[:, 4].copy()
-    
-            # Return processed data
-            if photFile.size > 16:
-                filepath = Path(fn).with_suffix('.npy')
-                photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree)
-                return filepath.name, photFile, photSkyCoord
-    
-            logger.debug(f"REJECT {fn}: Insufficient valid entries.")
+def process_convert_photometry_file(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger):
+    try:
+        # Read the file
+        photFile = genfromtxt(fn, dtype=float, delimiter=',', skip_header=1, invalid_raise=False)
+        logger.info(fn)
+
+        # Validate file shape and range
+        if photFile.size <= 16 or photFile.shape[1] != 8:
+            logger.debug(f"REJECT {fn}: Invalid size or shape.")
             return None
-        except Exception as e:
-            logger.error(f"Error processing file {fn}: {e}")
+
+        if max(photFile[:, 0]) >= 360 or max(photFile[:, 1]) >= 90:
+            logger.debug(f"REJECT {fn}: Invalid RA/Dec range.")
             return None
-    
-    
-    def process_files_multiprocessing(filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger):
-        # Pool for multiprocessing
-        with Pool() as pool:
-            results = pool.starmap(
-                process_file,
-                [(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger) for fn in filelist]
-            )
+
+        # Reject NaN entries
+        photFile = photFile[~isnan(photFile).any(axis=1)]
+
+        # Apply radial cut or image edge trimming
+        if racut != -99.9 and deccut != -99.9 and radiuscut != -99.9:
+            distanceArray = sqrt((photFile[:, 0] - racut)**2 + (photFile[:, 1] - deccut)**2)
+            photFile = delete(photFile, where(distanceArray > radiuscut / 60), axis=0)
+        else:
+            # Handle RA crossover
+            if (max(photFile[:, 0]) - min(photFile[:, 0])) > 180:
+                photFile[:, 0] = where(photFile[:, 0] > 180, photFile[:, 0] - 360, photFile[:, 0] + 360)
+
+            # Clip edges
+            raRange, decRange = ptp(photFile[:, 0]), ptp(photFile[:, 1])
+            raClip, decClip = raRange * ignoreedgefraction, decRange * ignoreedgefraction
+            photFile = photFile[
+                (photFile[:, 0] > min(photFile[:, 0]) + raClip) &
+                (photFile[:, 0] < max(photFile[:, 0]) - raClip) &
+                (photFile[:, 1] > min(photFile[:, 1]) + decClip) &
+                (photFile[:, 1] < max(photFile[:, 1]) - decClip)
+            ]
+
+        # Remove zero and low-count entries        
+        if len((photFile[:, 0] != 0) & (photFile[:, 1] != 0)) == 0:
+            print(f"REJECTED {fn}: Very likely no wcs fit - RA and Dec columns are zero.")
         
-        # Filter results
-        results = [res for res in results if res is not None]
-        new_files, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], [])
-        return list(new_files), list(photFileHolder), list(photSkyCoord)
-    
-    
-    # Example usage
-    if __name__ == "__main__":
-        new_files, photFileHolder, photSkyCoord = process_files_multiprocessing(
-            filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger
+        photFile = photFile[(photFile[:, 4] > lowestcounts) & (photFile[:, 0] != 0) & (photFile[:, 1] != 0)]
+
+        # Prepare output
+        photFile = c_[photFile, zeros((len(photFile), 4))]
+        photFile[:, 8:10] = nan
+        photFile[:, 11] = photFile[:, 4].copy()
+
+        # Return processed data
+        if photFile.size > 16:
+            filepath = Path(fn).with_suffix('.npy')
+            photSkyCoord = SkyCoord(ra=photFile[:, 0] * u.degree, dec=photFile[:, 1] * u.degree)
+            return filepath.name, photFile, photSkyCoord
+
+        logger.debug(f"REJECT {fn}: Insufficient valid entries.")
+        return None
+    except Exception as e:
+        logger.error(f"Error processing file {fn}: {e}")
+        return None
+
+
+def process_convert_files_multiprocessing(filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger):
+    # Pool for multiprocessing
+    with Pool(processes=max([cpu_count()-1,1])) as pool:
+        results = pool.starmap(
+            process_convert_photometry_file,
+            [(fn, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger) for fn in filelist]
         )
     
+    # Filter results
+    results = [res for res in results if res is not None]
+    new_files, photFileHolder, photSkyCoord = zip(*results) if results else ([], [], [])
+    return list(new_files), list(photFileHolder), list(photSkyCoord)
+
+
+
+def convert_photometry_files(filelist, ignoreedgefraction=0.05, lowestcounts=1800,  racut=-99.9, deccut=-99.9, radiuscut=-99.9):
+
+    new_files = []
+    photFileHolder=[]
+    photSkyCoord=[]
+    
+    new_files, photFileHolder, photSkyCoord = process_convert_files_multiprocessing(
+        filelist, racut, deccut, radiuscut, ignoreedgefraction, lowestcounts, logger
+    )    
 
     if new_files ==[] or photFileHolder==[] :
         raise AstrosourceException("Either there are no files of this photometry type or radial Cut seems to be outside of the range of your images... check your racut, deccut and radiuscut")
@@ -466,6 +355,100 @@ def gather_files(paths, filelist=None, filetype="fz", bjd=False, ignoreedgefract
 
     return phot_list, filterCode, photFileHolder, photSkyCoord
 
+def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger):
+    try:
+        imgRejFlag = 0
+        photReject = []
+
+        #breakpoint()
+
+        # Check minimum stars requirement
+        if referenceFrame.shape[0] < mincompstars:
+            return None, referenceFrame, photReject
+
+        # Validate photometry file
+        if photFile.size > imgsize and photFile.size > 7:
+            photRAandDec = photCoords[index]
+            rejectStars = []
+
+            if (max(photFile[:, 0]) <= 360) and (photFile[0][0] != 'null') and (photFile[0][0] != 0.0):
+                testStars = SkyCoord(ra=referenceFrame[:, 0] * u.degree, dec=referenceFrame[:, 1] * u.degree)
+                idx, d2d, _ = testStars.match_to_catalog_sky(photRAandDec)
+                rejectStars = where(d2d.arcsecond > acceptDistance)[0]
+            else:
+                logger.debug(f"Image {index} rejected due to problematic entries.")
+                imgRejFlag = 1
+                photReject.append(index)
+
+            # Remove rejected stars from the reference frame
+            if rejectStars.size > 0:
+                if not ((len(rejectStars) / referenceFrame.shape[0]) > starreject and index > rejectStart):
+                    referenceFrame = delete(referenceFrame, rejectStars, axis=0)
+                    logger.debug(f"Image {index}: Removed {len(rejectStars)} stars.")
+                else:
+                    logger.debug(f"Image {index} rejected due to too high a fraction of rejected stars.")
+                    imgRejFlag = 1
+                    photReject.append(index)
+
+            if imgRejFlag == 0:
+                logger.debug(f"Image {index}: All stars present.")
+        elif photFile.size <= 7:
+            logger.error(f"Image {index}: WCS coordinates broken or too few stars.")
+            photReject.append(index)
+        else:
+            logger.debug(f"Image {index}: Contains too few stars.")
+            photReject.append(index)
+
+        return imgRejFlag, referenceFrame, photReject
+    except Exception as e:
+        logger.error(f"Error processing image {index}: {e}")
+        return None, referenceFrame, [index]
+
+
+def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger):
+        
+    # Hack to get windows to not multiprocess until I figure out how to do it.    
+    if platform.system() == "Windows":
+        results=[]
+        for index in range(len(photFileHolder)):
+            results.append(process_phot_file(index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger))
+    
+    else:    
+        with Pool(processes=max([cpu_count()-1,1])) as pool:
+            # Prepare arguments
+            args = [
+                (index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger)
+                for index in range(len(photFileHolder))
+            ]
+    
+            # Process files in parallel
+            results = pool.starmap(process_phot_file, args)
+    
+    # Collect results
+    updated_referenceFrame = referenceFrame
+    photReject = []
+
+    for imgRejFlag, refFrameUpdate, rejects in results:
+        if refFrameUpdate is not None:
+            updated_referenceFrame = refFrameUpdate
+        photReject.extend(rejects)
+        
+    return updated_referenceFrame, photReject
+
+
+def evaluate_file_size(file_index, file_list, phot_file_holder, phot_coords):
+    """Evaluate a file to determine its size and suitability."""
+    file = file_list[file_index]
+    phot_file = phot_file_holder[file_index]
+    file_size = phot_file.size
+    file_ra_dec = phot_coords[file_index]
+    return file_size, file_index, file_ra_dec
+
+
+def compute_size(phot_file):
+    """Compute the size for a specific file."""
+    return phot_file.size
+
 def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskymapper=False,closerejectd=5.0, photCoords=None, photFileHolder=None, mincompstars=0.1, mincompstarstotal=-99, starreject=0.1 , acceptDistance=1.0, lowcounts=2000, hicounts=3000000, imageFracReject=0.0,  rejectStart=3, maxcandidatestars=10000, restrictcompcolourcentre=-99.0, restrictcompcolourrange=-99.0, filterCode=None, restrictmagbrightest=-99.0, restrictmagdimmest=99.0, minfractionimages=0.5):
     """
     Finds stars useful for photometry in each photometry/data file
@@ -510,22 +493,22 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym
     if os.path.exists(paths['parent'] / "photFileHolder"):
         os.remove(paths['parent'] / "photFileHolder")
 
-    fileSizer=0
-    logger.info("Finding image with most stars detected and reject ones with bad WCS")
+       
+    logger.info("Finding image with most stars detected and rejecting ones with bad WCS")
     referenceFrame = None
-
-    counter=0
-    for file in fileList:
-        photFile = photFileHolder[counter]
-
-        # Sort through and find the largest file and use that as the reference file
-        if photFile.size > fileSizer:            
-            referenceFrame = photFile
-            fileSizer = photFile.size
-            fileRaDec = photCoords[counter]
-            logger.debug("{} - {}".format(photFile.size, file))
-        counter=counter+1
-
+    #file_sizer = 0
+    fileRaDec = None
+    
+    list_of_sizes=[]
+    for filething in photFileHolder:
+        list_of_sizes.append(filething.size)
+    
+    largest_file_index= argmax(list_of_sizes)
+    
+    referenceFrame=photFileHolder[largest_file_index]
+    fileRaDec = photCoords[largest_file_index]
+    fileSizer=referenceFrame.size
+    
     if not referenceFrame.size:
         raise AstrosourceException("No suitable reference files found")
 
@@ -558,7 +541,6 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym
             lowcounts=lowcounts+(0.05*lowcounts)
             hicounts=hicounts-(0.05*hicounts)
 
-
     logger.debug("Number of stars post")
     referenceFrame = delete(referenceFrame, rejectStars, axis=0)
     logger.debug("Final count range, Low: " +str(int(lowcounts))+ " High: "+ str(int(hicounts)))
@@ -587,169 +569,18 @@ def find_stars(targets, paths, fileList, nopanstarrs=False, nosdss=False, noskym
     while (compchecker < mincompstars): # Keep going until you get the minimum number of Comp Stars
         imgsize=imageFracReject * fileSizer # set threshold size
         rejStartCounter = 0
-        #usedImages=[] # Set up used images array
         imgReject = 0 # Number of images rejected due to high rejection rate
         loFileReject = 0 # Number of images rejected due to too few stars in the photometry file
         wcsFileReject=0
         imgOverride=0
-        # q=0
-        # photReject=[]
-
-        # for Nholder in range(len(photFileHolder)):
-
-        #     if ( not referenceFrame.shape[0] < mincompstars):
-        #         rejStartCounter = rejStartCounter +1
-        #         #photFile = load(paths['parent'] / file)
-        #         photFile = photFileHolder[q]
-        #         logger.debug('Image Number: ' + str(rejStartCounter))
-        #         logger.debug(file)
-        #         logger.debug("Image threshold size: "+str(imgsize))
-        #         logger.debug("Image catalogue size: "+str(photFile.size))
-        #         imgRejFlag=0
-        #         if photFile.size > imgsize and photFile.size > 7 :
-        #             phottmparr = asarray(photFile)
-        #             photRAandDec = photCoords[q]
-
-        #             # Checking existance of stars in all photometry files
-        #             rejectStars=[] # A list to hold what stars are to be rejected
-
-        #             if (( phottmparr[:,0] > 360).sum() == 0) and ( phottmparr[0][0] != 'null') and ( phottmparr[0][0] != 0.0) :
-        #                 # Faster way than loop
-        #                 testStars=SkyCoord(ra = referenceFrame[:,0]*u.degree, dec = referenceFrame[:,1]*u.degree)
-        #                 idx, d2d, _ = testStars.match_to_catalog_sky(photRAandDec)
-        #                 rejectStars=where(d2d.arcsecond > acceptDistance)[0]
-
-        #             else:
-        #                 logger.debug('**********************')
-        #                 logger.debug('Image Rejected due to problematic entries in Photometry File')
-        #                 logger.debug('**********************')
-        #                 imgReject=imgReject+1
-        #                 photReject.append(q)
-        #                 imgRejFlag=1
-
-        #             # if the rejectstar list is not empty, remove the stars from the reference List
-        #             if len(rejectStars) != 0:
-
-        #                 if not (((len(rejectStars) / referenceFrame.shape[0]) > starreject) and rejStartCounter > rejectStart):
-        #                     referenceFrame = delete(referenceFrame, rejectStars, axis=0)
-        #                     logger.debug('**********************')
-        #                     logger.debug('Stars Removed  : ' +str(len(rejectStars)))
-        #                     logger.debug('Remaining Stars: ' +str(referenceFrame.shape[0]))
-        #                     logger.debug('**********************')
-
-        #                 else:
-        #                     logger.debug('**********************')
-        #                     logger.debug('Image Rejected due to too high a fraction of rejected stars')
-        #                     logger.debug(len(rejectStars) / referenceFrame.shape[0])
-        #                     logger.debug('**********************')
-        #                     imgReject=imgReject+1
-        #                     photReject.append(q)
-
-        #             elif imgRejFlag==0:
-        #                 logger.debug('**********************')
-        #                 logger.debug('All Stars Present')
-        #                 logger.debug('**********************')
-
-        #             # If we have removed all stars, we have failed!
-        #             if (referenceFrame.shape[0]==0):
-        #                 logger.error("Problem file - {}".format(file))
-        #                 logger.error("Running Loop again")
-
-        #         elif photFile.size < 7:
-        #             logger.error('**********************')
-        #             logger.error("WCS Coordinates broken")
-        #             logger.error('**********************')
-        #             wcsFileReject=wcsFileReject+1
-        #             photReject.append(q)
-        #         else:
-        #             logger.debug('**********************')
-        #             logger.debug("CONTAINS TOO FEW STARS")
-        #             logger.debug('**********************')
-        #             loFileReject=loFileReject+1
-        #             photReject.append(q)
-
-        #     q=q+1
-        
-        def process_phot_file(index, photFile, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger):
-            try:
-                imgRejFlag = 0
-                photReject = []
-        
-                # Check minimum stars requirement
-                if referenceFrame.shape[0] < mincompstars:
-                    return None, referenceFrame, photReject
-        
-                # Validate photometry file
-                if photFile.size > imgsize and photFile.size > 7:
-                    photRAandDec = photCoords[index]
-                    rejectStars = []
-        
-                    if (max(photFile[:, 0]) <= 360) and (photFile[0][0] != 'null') and (photFile[0][0] != 0.0):
-                        testStars = SkyCoord(ra=referenceFrame[:, 0] * u.degree, dec=referenceFrame[:, 1] * u.degree)
-                        idx, d2d, _ = testStars.match_to_catalog_sky(photRAandDec)
-                        rejectStars = where(d2d.arcsecond > acceptDistance)[0]
-                    else:
-                        logger.debug(f"Image {index} rejected due to problematic entries.")
-                        imgRejFlag = 1
-                        photReject.append(index)
-        
-                    # Remove rejected stars from the reference frame
-                    if rejectStars.size > 0:
-                        if not ((len(rejectStars) / referenceFrame.shape[0]) > starreject and index > rejectStart):
-                            referenceFrame = delete(referenceFrame, rejectStars, axis=0)
-                            logger.debug(f"Image {index}: Removed {len(rejectStars)} stars.")
-                        else:
-                            logger.debug(f"Image {index} rejected due to too high a fraction of rejected stars.")
-                            imgRejFlag = 1
-                            photReject.append(index)
         
-                    if imgRejFlag == 0:
-                        logger.debug(f"Image {index}: All stars present.")
-                elif photFile.size <= 7:
-                    logger.error(f"Image {index}: WCS coordinates broken or too few stars.")
-                    photReject.append(index)
-                else:
-                    logger.debug(f"Image {index}: Contains too few stars.")
-                    photReject.append(index)
-        
-                return imgRejFlag, referenceFrame, photReject
-            except Exception as e:
-                logger.error(f"Error processing image {index}: {e}")
-                return None, referenceFrame, [index]
-        
-        
-        def process_photometry_files_multiprocessing(photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger):
-            with Pool() as pool:
-                # Prepare arguments
-                args = [
-                    (index, photFileHolder[index], photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger)
-                    for index in range(len(photFileHolder))
-                ]
-        
-                # Process files in parallel
-                results = pool.starmap(process_phot_file, args)
-        
-            # Collect results
-            updated_referenceFrame = referenceFrame
-            photReject = []
-        
-            for imgRejFlag, refFrameUpdate, rejects in results:
-                if refFrameUpdate is not None:
-                    updated_referenceFrame = refFrameUpdate
-                photReject.extend(rejects)
-        
-            return updated_referenceFrame, photReject
-        
-        
-        # Example usage
         updated_referenceFrame, photReject = process_photometry_files_multiprocessing(
             photFileHolder, photCoords, referenceFrame, acceptDistance, starreject, rejectStart, mincompstars, imgsize, logger
         )
 
-        # Remove files and Hold the photSkyCoords in memory
 
+        # Remove files and Hold the photSkyCoords in memory
         photCoords=asarray(photCoords, dtype=object)
-
         photCoords=delete(photCoords, photReject, axis=0)
         photFileHolder=asarray(photFileHolder, dtype=object)
         photFileHolder=delete(photFileHolder, photReject, axis=0)
@@ -838,7 +669,6 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen
             photCoords=originalphotSkyCoord
             photFileHolder=originalphotFileHolder
 
-
         elif starreject == 0.15 and imageFracReject == 0.8 and mincompstars !=1:
             logger.error("Maximum number of Candidate Comparison Stars found this cycle: " + str(compchecker))
             logger.error("Failed to find sufficient comparison candidates with the maximum restrictions, trying with a lower value for mincompstars")
@@ -990,9 +820,9 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen
                             'SkyMapper' : {'filter' : 'iPSF', 'error' : 'e_iPSF', 'colmatch' : 'rPSF', 'colerr' : 'e_rPSF', 'colname' : 'r-i', 'colrev' : '1'},
                             'PanSTARRS': {'filter' : 'imag', 'error' : 'e_imag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-i', 'colrev' : '1'},
                             'APASS' : {'filter' : 'i_mag', 'error' : 'e_i_mag', 'colmatch' : 'r_mag', 'colerr' : 'e_r_mag', 'colname' : 'r-i', 'colrev' : '1'}},
-                    'zs' : {'PanSTARRS': {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'},
+                    'zs' : {'SDSS' : {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'},
                             'SkyMapper' : {'filter' : 'zPSF', 'error' : 'e_zPSF', 'colmatch' : 'rPSF', 'colerr' : 'e_rPSF', 'colname' : 'r-zs', 'colrev' : '1'},
-                            'SDSS' : {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'}},
+                            'PanSTARRS': {'filter' : 'zmag', 'error' : 'e_zmag', 'colmatch' : 'rmag', 'colerr' : 'e_rmag', 'colname' : 'r-zs', 'colrev' : '1'}},
                     }
 
         try:
@@ -1013,7 +843,6 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen
 
                         coords = catalogue_call(avgCoord, 1.5*radius, opt, cat_name, targets=targets, closerejectd=closerejectd)
                         # If no results try next catalogue
-                        #print (len(coords.ra))
                         if len(coords.ra) == 0:
                             coords=[]
                             raise AstrosourceException("Empty catalogue produced from catalogue call")
@@ -1095,4 +924,8 @@ def process_photometry_files_multiprocessing(photFileHolder, photCoords, referen
             f.write(str(filename) +"\n")
 
 
-    return fileList, outputComps, photFileHolder, photCoords
\ No newline at end of file
+    return fileList, outputComps, photFileHolder, photCoords
+
+# Needed for windows to multiprocess appropriately
+if __name__ == "__main__":
+    pass
\ No newline at end of file
diff --git a/astrosource/main.py b/astrosource/main.py
index 0b48688..18df835 100644
--- a/astrosource/main.py
+++ b/astrosource/main.py
@@ -2,7 +2,8 @@
 import click
 import sys
 import logging
-
+import csv
+import os
 from numpy import array, all
 from astropy.utils.exceptions import AstropyWarning, AstropyDeprecationWarning
 
@@ -55,7 +56,7 @@
 
 
 
-@click.option('--period', is_flag=True, type=float, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.')
+@click.option('--period', is_flag=True, default=False, help='Search for periodicity in the data, currently with PDM and String methods. This will autoselect a reasonable search range if not provided a range.')
 @click.option('--periodlower', '-pl', type=float, default=-99.9, help='Shortest period to trial in days. Default is 0.05 days')
 @click.option('--periodupper', '-pu', type=float, default=-99.9, help='Longest period to trial in days. Default is one-third the observational baseline of your dataset.')
 @click.option('--periodtests', '-pt', type=int, default=10000, help='Number of different trial periods to run')
@@ -126,12 +127,8 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree
             logger.info('All output files removed')
             return
         if not (ra and dec) and not target_file and not variablehunt:
-            #logger.error("Either RA and Dec or a targetfile must be specified")
             logger.error("No specified RA or Dec nor targetfile nor request for a variable hunt. It is assumed you have no target to analyse.")
             notarget=True
-            #return
-
-
 
         if ra and dec:
             ra, dec = convert_coords(ra, dec)
@@ -217,16 +214,18 @@ def main(full, stars, comparison, variablehunt, notarget, lowestcounts, usescree
             ts.find_variables()
 
         if variablehunt == True:
-            targets = get_targets(parentPath / 'results/potentialVariables.csv')
-
-
-
+            if os.path.exists(parentPath / 'results/potentialVariables.csv'):
+                targets = get_targets(parentPath / 'results/potentialVariables.csv')
+            else:
+                targets=None
+                
         if targets is not None:
             if full or phot:
                 ts.photometry(filesave=True, targets=targets)
             if full or plot:
                 ts.plot(detrend=detrend, period=period, eebls=eebls, filesave=True)
 
+        
         sys.stdout.write("✅ AstroSource analysis complete\n")
 
     except AstrosourceException as e:
diff --git a/astrosource/periodic.py b/astrosource/periodic.py
index 10263c4..410a496 100644
--- a/astrosource/periodic.py
+++ b/astrosource/periodic.py
@@ -5,10 +5,12 @@
 import sys
 import os
 
+import traceback
+
 import matplotlib
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
-from multiprocessing import Pool
+from multiprocessing import Pool, cpu_count
 from functools import partial
 import logging
 
@@ -17,7 +19,7 @@
 from astropy.timeseries import LombScargle
 
 logger = logging.getLogger('astrosource')
-NCPUS = 1
+#NCPUS = 1
 
 # Note that the functions that calculate the ANOVA periodograms have been adapted from the astrobase codeset
 # These are aov_theta, resort_by_time, get_frequency_grid, sigclip_magseries, phase_magseries, aov_periodfind, phase_magseries_with_errs, aovhm_theta, aovhm_periodfind
@@ -36,47 +38,279 @@
 # SOFTWARE.
 
 def process_file(file, periodsteps, minperiod, maxperiod, numBins, periodPath, paths, filterCode, logger):
+        
     try:
-        logger.debug(file)
-        variableName = file.stem.split('_')[0]
+        trialRange=[minperiod, maxperiod]        
+        variableName=file.stem.split('_')[0]
         logger.debug("Variable Name: {}".format(variableName))
-
-        # Load data
         varData = genfromtxt(file, dtype=float, delimiter=',')
-        calibFile = file.parent / "{}{}".format(file.stem.replace('diff', 'calib'), file.suffix)
+        calibFile = file.parent / "{}{}".format(file.stem.replace('diff','calib'), file.suffix)
         logger.debug(calibFile)
-
-        calibData = None
         if calibFile.exists():
-            calibData = genfromtxt(calibFile, dtype=float, delimiter=',')
+            calibData=genfromtxt(calibFile, dtype=float, delimiter=',')
 
-        # Determine which dataset to use for analysis
-        if calibFile.exists() and calibData is not None and calibData.size > 3:
-            dataToUse = calibData
-            logger.info("Using calibrated data")
+        if calibFile.exists():
+            if (calibData.size > 3):
+                pdm=phase_dispersion_minimization(calibData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
+            else:
+                logger.info("Calibration File not large enough to run period methods on Calibrated data, running on differential data")
+                pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
         else:
-            dataToUse = varData
-            logger.info("Using differential data")
+            pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
 
-        # Phase dispersion minimization
-        pdm = phase_dispersion_minimization(dataToUse, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
-
-        # Save plots and results
         plt.figure(figsize=(15, 5))
-        plt.plot(pdm["periodguess_array"], pdm["distance_results"])
-        plt.gca().invert_yaxis()
-        plt.title("Range {0} d  Steps: {1}".format(maxperiod - minperiod, periodsteps))
-        plt.xlabel(r"Trial Period")
-        plt.ylabel(r"Likelihood of Period")
-        plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png")
-        plt.clf()
+
 
         with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:
-            f.write(f"Variable: {variableName}\n")
-            f.write(f"Distance Method Estimate (days): {pdm['distance_minperiod']}\n")
-            f.write(f"Distance method error: {pdm['distance_error']}\n")
+            f.write("Variable : "+str(variableName) +"\n")
+
+
+        try:
+            logger.debug("Distance Method Estimate (days): " + str(pdm["distance_minperiod"]))
+            logger.debug("Distance method error: " + str(pdm["distance_error"]))
+            pdmfailed=False
+        except:
+            logger.debug("Distance Method Failed")
+            pdmfailed=True
+            
+        if not pdmfailed:
+            with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:
+                f.write("Variable : "+str(variableName) +"\n")
+                f.write("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])+"\n")
+                f.write("Distance method error: " + str(pdm["distance_error"])+"\n")
+            
+            plt.plot(pdm["periodguess_array"], pdm["distance_results"])
+            plt.gca().invert_yaxis()
+            plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
+            plt.xlabel(r"Trial Period")
+            plt.ylabel(r"Likelihood of Period")
+            plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png")
+            plt.clf()
+            
+            if (varData.size > 3):
+                phaseTest=(varData[:,0] / (pdm["distance_minperiod"])) % 1        
+        
+                plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None')
+                plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None')
+                plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None')
+                plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None')
+                plt.gca().invert_yaxis()
+                plt.title("Period: {0} d  Steps: {1}".format(pdm["distance_minperiod"], periodsteps))
+                plt.xlabel(r"Phase ($\phi$)")
+                plt.ylabel(f"Differential {filterCode} Magnitude")
+                plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot.png")
+                plt.clf()
+    
+            if calibFile.exists():
+                if (calibData.size > 3):
+                    phaseTestCalib=(calibData[:,0] / (pdm["distance_minperiod"])) % 1
+                    plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None')
+                    plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None')
+                    plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None')
+                    plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None')
+                    plt.gca().invert_yaxis()
+                    plt.title("Period: {0} d  Steps: {1}".format(pdm["distance_minperiod"], periodsteps))
+                    plt.xlabel(r"Phase ($\phi$)")
+                    plt.ylabel(f"Calibrated {filterCode} Magnitude")
+                    plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot_Calibrated.png")
+                    plt.clf()
+        
+                    tempPeriodCatOut=[]
+                    for g in range(len(calibData[:,0])):
+                        tempPeriodCatOut.append([(calibData[g,0]/(pdm["distance_minperiod"]) % 1), calibData[g,1], calibData[g,2]])
+                    tempPeriodCatOut=asarray(tempPeriodCatOut)
+                    savetxt(periodPath / f"{variableName}_String_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
+    
+            if (varData.size > 3):
+    
+                tempPeriodCatOut=[]
+                for g in range(len(phaseTest)):
+                    tempPeriodCatOut.append([phaseTest[g],varData[g,1]])
+                tempPeriodCatOut=asarray(tempPeriodCatOut)
+                savetxt(periodPath / f"{variableName}_StringTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
+        
+                tempPeriodCatOut=[]
+                for g in range(len(varData[:,0])):
+                    tempPeriodCatOut.append([(varData[g,0]/(pdm["distance_minperiod"]) % 1), varData[g,1], varData[g,2]])
+                tempPeriodCatOut=asarray(tempPeriodCatOut)
+                savetxt(periodPath / f"{variableName}_String_PhasedDiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
+    
+            if np.isnan(pdm["stdev_results"][0]) or pdm["stdev_results"][0] == 0.0 or (varData.size < 4):
+                logger.info("No PDM results due to lack of datapoint coverage")
+            else:
+                logger.debug("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"]))
+                phaseTest=(varData[:,0] / (pdm["stdev_minperiod"])) % 1
+                logger.debug("PDM method error: " + str(pdm["stdev_error"]))
+    
+                with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:
+                    f.write("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])+"\n")
+                    f.write("PDM method error: " + str(pdm["stdev_error"])+"\n\n")
+    
+    
+            plt.plot(pdm["periodguess_array"], pdm["stdev_results"])
+            plt.gca().invert_yaxis()
+            plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
+            plt.xlabel(r"Trial Period")
+            plt.ylabel(r"Likelihood of Period")
+            plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot.png")
+    
+            plt.clf()
+    
+            if (varData.size > 3):
+                plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None')
+                plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None')
+                plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None')
+                plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None')
+                plt.gca().invert_yaxis()
+                plt.title("Period: {0} d  Steps: {1}".format(pdm["stdev_minperiod"], periodsteps))
+                plt.xlabel(r"Phase ($\phi$)")
+                plt.ylabel(r"Differential " + str(filterCode) + " Magnitude")
+                plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot.png")
+                plt.clf()
+    
+            if calibFile.exists():
+                if (calibData.size > 3):
+                    phaseTestCalib=(calibData[:,0] / (pdm["stdev_minperiod"])) % 1
+                    plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None')
+                    plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None')
+                    plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None')
+                    plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None')
+                    plt.gca().invert_yaxis()
+                    plt.title("Period: {0} d  Steps: {1}".format(pdm["stdev_minperiod"], periodsteps))
+                    plt.xlabel(r"Phase ($\phi$)")
+                    plt.ylabel(r"Calibrated " + str(filterCode) + " Magnitude")
+                    plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot_Calibrated.png")
+                    plt.clf()
+        
+                    tempPeriodCatOut=[]
+                    for g in range(len(calibData[:,0])):
+                        tempPeriodCatOut.append([(calibData[g,0]/(pdm["stdev_minperiod"])) % 1, calibData[g,1], calibData[g,2]])
+                    tempPeriodCatOut=asarray(tempPeriodCatOut)
+                    savetxt(periodPath / f"{variableName}_PDM_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
+    
+            if (varData.size > 3):
+    
+                tempPeriodCatOut=[]
+                for g in range(len(phaseTest)):
+                    tempPeriodCatOut.append([phaseTest[g],varData[g,1]])
+                tempPeriodCatOut=asarray(tempPeriodCatOut)
+                savetxt(periodPath / f"{variableName}_PDMTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
+        
+                tempPeriodCatOut=[]
+                for g in range(len(varData[:,0])):
+                    tempPeriodCatOut.append([(varData[g,0]/(pdm["stdev_minperiod"])) % 1, varData[g,1], varData[g,2]])
+                tempPeriodCatOut=asarray(tempPeriodCatOut)
+                savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
+    
+            # Plot publication plots    
+            plt.figure(figsize=(5, 3))
+    
+            plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5)
+            plt.gca().invert_yaxis()
+            plt.xlabel(r"Trial Period")
+            plt.ylabel(r"Likelihood of Period")
+            plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4)
+            plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.png", dpi=300)
+            plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.eps")
+    
+            plt.clf()
+    
+            plt.figure(figsize=(5, 3))
+    
+            plt.plot(pdm["periodguess_array"], pdm["distance_results"], linewidth=0.5)
+            plt.gca().invert_yaxis()
+            plt.xlabel(r"Trial Period")
+            plt.ylabel(r"Likelihood of Period")
+            plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4)
+            plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.png", dpi=300)
+            plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.eps")
+    
+            plt.clf()
+
+        # ANOVA        
+        if calibFile.exists():
+            if (calibData.size > 3):
+                if len(calibData[:,0]) < 75:
+                    binsize=0.1
+                else:
+                    binsize=0.05
+                minperbin=int((len(calibData[:,0])/10))
+            elif (varData.size > 3):
+                if len(varData[:,0]) < 75:
+                    binsize=0.1
+                else:
+                    binsize=0.05
+                minperbin=int((len(varData[:,0])/10))
+                
+        else:
+            if (varData.size > 3):
+                if len(varData[:,0]) < 75:
+                    binsize=0.1
+                else:
+                    binsize=0.05
+                minperbin=int((len(varData[:,0])/10))
+        
+        if 'minperbin' in locals():
+            if minperbin > 10:
+                minperbin=10
+        else:
+            minperbin = 3
+            
+        
+        # Theta Anova Method off for the moment until I put in a command-line option        
+        # if calibFile.exists():
+        #     if (calibData.size > 3):
+        #         aovoutput=aov_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
+        # else:
+        #     if (varData.size > 3):
+        #         aovoutput=aov_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
 
+        # logger.debug("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"]))
+        #with open(paths['parent'] / "periodEstimates.txt", "a+") as f:            
+        #    f.write("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])+"\n")
+        
+        if calibFile.exists():
+            if (calibData.size > 3):
+                aovhmoutput=aovhm_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
+            elif (varData.size > 3):                
+                aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
+                
+        else:
+            if (varData.size > 3):
+                aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
+            
+        try:
+            logger.debug("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"]))
+            with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:            
+                f.write("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])+"\n")
+        except:
+            logger.debug("Harmonic Anova Method Estimate Failed")
+
+        # LOMB SCARGLE
+        for nts in range(2):
+            if calibFile.exists():
+                if (calibData.size > 3):
+                    lscargoutput = LombScargleMultiterm('periodifile', (calibData[:, 0]), (calibData[:, 1]), (calibData[:, 2]),
+                                                        nterms=nts+1,
+                                                        periodlower=minperiod, periodupper=maxperiod, samples=20,
+                                                        disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
+        
+                    logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput))
+                    with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:            
+                        f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n")
+            else:
+                if (varData.size > 3):
+                    lscargoutput = LombScargleMultiterm('periodifile', (varData[:, 0]), (varData[:, 1]), (varData[:, 2]),
+                                                        nterms=nts+1,
+                                                        periodlower=minperiod, periodupper=maxperiod, samples=20,
+                                                        disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
+        
+                    logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput))
+                    with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:            
+                        f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n")
+        
     except Exception as e:
+        print(traceback.format_exc())
         logger.error(f"Error processing file {file}: {e}")
 
 def aov_theta(times, mags, errs, frequency,
@@ -302,7 +536,6 @@ def sigclip_magseries(times, mags, errs,
     # fake the errors if they don't exist
     # this is inconsequential to sigma-clipping
     # we don't return these dummy values if the input errs are None
-    #print (errs)
     if errs is None:
         # assume 0.1% errors if not given
         # this should work for mags and fluxes
@@ -783,16 +1016,6 @@ def aov_periodfind(times,
             # default end period is length of time series
             startf = 1.0/(stimes.max() - stimes.min())
 
-        # # if we're not using autofreq, then use the provided frequencies
-        # if not autofreq:
-        #     frequencies = nparange(startf, endf, stepsize)
-
-        # else:
-        #     # this gets an automatic grid of frequencies to use
-        #     frequencies = get_frequency_grid(stimes,
-        #                                      minfreq=startf,
-        #                                      maxfreq=endf)
-
         stepsize = (endf-startf)/periodsteps
         frequencies = nparange(startf, endf, stepsize)
 
@@ -814,8 +1037,6 @@ def aov_periodfind(times,
         periods = 1.0/frequencies
 
         plt.plot(periods, lsp)
-        #plt.gca().invert_yaxis()
-        #plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
         plt.xlabel(r"Trial Period")
         plt.ylabel(r"Likelihood of Period")
         plt.savefig(periodPath / f"{variableName}_ANOVAThetaLikelihoodPlot.png")
@@ -882,10 +1103,6 @@ def aov_periodfind(times,
             perioddiff = abs(period - prevperiod)
             bestperiodsdiff = [abs(period - x) for x in nbestperiods]
 
-            # print('prevperiod = %s, thisperiod = %s, '
-            #       'perioddiff = %s, peakcount = %s' %
-            #       (prevperiod, period, perioddiff, peakcount))
-
             # this ensures that this period is different from the last
             # period and from all the other existing best periods by
             # periodepsilon to make sure we jump to an entire different peak
@@ -1298,26 +1515,9 @@ def aovhm_periodfind(times,
             # default end period is length of time series
             startf = 1.0/(stimes.max() - stimes.min())
 
-
-    
-
-        # # if we're not using autofreq, then use the provided frequencies
-        # if not autofreq:
-        #     frequencies = nparange(startf, endf, stepsize)
-
-        # else:
-        #     # this gets an automatic grid of frequencies to use
-        #     frequencies = get_frequency_grid(stimes,
-        #                                      minfreq=startf,
-        #                                      maxfreq=endf)
-
         stepsize = (endf-startf)/periodsteps
         frequencies = nparange(startf, endf, stepsize)
 
-        # map to parallel workers
-        if (not nworkers) or (nworkers > NCPUS):
-            nworkers = NCPUS
-
         # renormalize the working mags to zero and scale them so that the
         # variance = 1 for use with our LSP functions
         if normalize:
@@ -1331,9 +1531,6 @@ def aovhm_periodfind(times,
         magvariance_bot = (nmags.size - 1)*npsum(1.0/(serrs*serrs)) / nmags.size
         magvariance = magvariance_top/magvariance_bot
 
-        #tasks = [(stimes, nmags, serrs, x, nharmonics, magvariance)
-        #         for x in frequencies]
-
         lsp=[]
         for x in frequencies:
             lsp.append(aovhm_theta(times, mags, errs, x, nharmonics, magvariance))
@@ -1342,8 +1539,6 @@ def aovhm_periodfind(times,
         periods = 1.0/frequencies
 
         plt.plot(periods, lsp)
-        #plt.gca().invert_yaxis()
-        #plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
         plt.xlabel(r"Trial Period")
         plt.ylabel(r"Likelihood of Period")
         plt.savefig(periodPath / f"{variableName}_ANOVAharmonic_LikelihoodPlot.png")
@@ -1410,10 +1605,6 @@ def aovhm_periodfind(times,
             perioddiff = abs(period - prevperiod)
             bestperiodsdiff = [abs(period - x) for x in nbestperiods]
 
-            # print('prevperiod = %s, thisperiod = %s, '
-            #       'perioddiff = %s, peakcount = %s' %
-            #       (prevperiod, period, perioddiff, peakcount))
-
             # this ensures that this period is different from the last
             # period and from all the other existing best periods by
             # periodepsilon to make sure we jump to an entire different peak
@@ -1434,8 +1625,6 @@ def aovhm_periodfind(times,
         plt.gca().invert_yaxis()
         plt.title("Period: " + str(finperiods[bestperiodind]))
         plt.xlim(-0.01,2.01)
-        #plt.gca().invert_yaxis()
-        #plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
         plt.xlabel(r"Phase")
         plt.ylabel(r"Magnitude")
         plt.savefig(periodPath / f"{variableName}_ANOVAHarmonicLightcurve.png")
@@ -1578,7 +1767,13 @@ def phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, nu
     distance_results = []
     stdev_results = []
 
+
+    if len(varData) < 5:
+        return
+
+    #try:
     (julian_dates, fluxes) = (varData[:,0],varData[:,1])
+    
     normalizedFluxes = normalize(fluxes)
 
     for r in range(periodsteps):
@@ -1646,8 +1841,6 @@ def phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, nu
                 break
             stepper=stepper+1
 
-
-        #print ("Stdev method error: " + str((righthandP - lefthandP)/2))
         pdm["stdev_error"] = (righthandP - lefthandP)/2
 
 
@@ -1689,10 +1882,7 @@ def phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, nu
 
 def LombScargleMultiterm(infile, t, m, d, periodlower=0.2, periodupper=2.5, nterms=1, multisearch=False, samples=5,
                          disablelightcurve=False, periodPath=False, variableName="NoName", periodsteps=10000):
-    #print(
-    #    'using ' + str(samples) + ' samples per peak, start P = ' + str(periodlower) + ', end P = ' + str(periodupper))
-    # Calculate the Lomb-Scargle periodogram values
-        
+    
     ls = LombScargle(t, m, d, nterms=nterms, fit_mean=True)
 
     # create equally spaced test in time space, then swap to unequal spaces in frequency    
@@ -1700,10 +1890,7 @@ def LombScargleMultiterm(infile, t, m, d, periodlower=0.2, periodupper=2.5, nter
     freq = 1/freq
 
     
-    try:
-        
-        #freq, power = ls.autopower(samples_per_peak=samples, minimum_frequency=1 / periodupper,
-        #                           maximum_frequency=1 / periodlower)
+    try:        
         power = ls.power(freq)
     except:
         print ("Lomb Scargle failed.") # Need to hunt down a very rare but existant memory problem
@@ -1793,7 +1980,7 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1
     logger.info("Maximum Period Tested  : " +str(maxperiod))
     logger.info("Number of Period Trials: " +str(periodsteps))
 
-    trialRange=[minperiod, maxperiod]
+    #trialRange=[minperiod, maxperiod]
 
     # Get list of phot files
     periodPath = paths['periods']
@@ -1820,274 +2007,11 @@ def plot_with_period(paths, filterCode, numBins = 10, minperiod=0.2, maxperiod=1
     )
 
     # Use multiprocessing
-    with Pool(processes=os.cpu_count()) as pool:
+    with Pool(processes=max([cpu_count()-1,1])) as pool:
         pool.map(worker, fileList)
-
-    # Load in the files
-    # for file in fileList:
-    #     logger.debug(file)
-    #     variableName=file.stem.split('_')[0]
-    #     #logger.debug(str(outcatPath).replace('//',''))
-    #     logger.debug("Variable Name: {}".format(variableName))
-    #     varData = genfromtxt(file, dtype=float, delimiter=',')
-    #     calibFile = file.parent / "{}{}".format(file.stem.replace('diff','calib'), file.suffix)
-    #     logger.debug(calibFile)
-    #     if calibFile.exists():
-    #         calibData=genfromtxt(calibFile, dtype=float, delimiter=',')
-
-    #         #print (calibData.size)
-    #     if calibFile.exists():
-    #         if (calibData.size > 3):
-    #             pdm=phase_dispersion_minimization(calibData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
-    #         else:
-    #             logger.info("Calibration File not large enough to run period methods on Calibrated data, running on differential data")
-    #             pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
-    #     else:
-    #         pdm=phase_dispersion_minimization(varData, periodsteps, minperiod, maxperiod, numBins, periodPath, variableName)
-
-    #     plt.figure(figsize=(15, 5))
-
-    #     logger.debug("Distance Method Estimate (days): " + str(pdm["distance_minperiod"]))
-    #     logger.debug("Distance method error: " + str(pdm["distance_error"]))
-        
-    #     with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:
-    #         f.write("Variable : "+str(variableName) +"\n")
-    #         f.write("Distance Method Estimate (days): " + str(pdm["distance_minperiod"])+"\n")
-    #         f.write("Distance method error: " + str(pdm["distance_error"])+"\n")
         
-    #     plt.plot(pdm["periodguess_array"], pdm["distance_results"])
-    #     plt.gca().invert_yaxis()
-    #     plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
-    #     plt.xlabel(r"Trial Period")
-    #     plt.ylabel(r"Likelihood of Period")
-    #     plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot.png")
-    #     plt.clf()
-        
-    #     if (varData.size > 3):
-    #         phaseTest=(varData[:,0] / (pdm["distance_minperiod"])) % 1        
-    
-    #         plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None')
-    #         plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None')
-    #         plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None')
-    #         plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None')
-    #         plt.gca().invert_yaxis()
-    #         plt.title("Period: {0} d  Steps: {1}".format(pdm["distance_minperiod"], periodsteps))
-    #         plt.xlabel(r"Phase ($\phi$)")
-    #         plt.ylabel(f"Differential {filterCode} Magnitude")
-    #         plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot.png")
-    #         plt.clf()
-
-    #     if calibFile.exists():
-    #         if (calibData.size > 3):
-    #             phaseTestCalib=(calibData[:,0] / (pdm["distance_minperiod"])) % 1
-    #             plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None')
-    #             plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None')
-    #             plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None')
-    #             plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None')
-    #             plt.gca().invert_yaxis()
-    #             plt.title("Period: {0} d  Steps: {1}".format(pdm["distance_minperiod"], periodsteps))
-    #             plt.xlabel(r"Phase ($\phi$)")
-    #             plt.ylabel(f"Calibrated {filterCode} Magnitude")
-    #             plt.savefig(periodPath / f"{variableName}_StringTestPeriodPlot_Calibrated.png")
-    #             plt.clf()
-    
-    #             tempPeriodCatOut=[]
-    #             for g in range(len(calibData[:,0])):
-    #                 tempPeriodCatOut.append([(calibData[g,0]/(pdm["distance_minperiod"]) % 1), calibData[g,1], calibData[g,2]])
-    #             tempPeriodCatOut=asarray(tempPeriodCatOut)
-    #             savetxt(periodPath / f"{variableName}_String_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
-
-    #     if (varData.size > 3):
-
-    #         tempPeriodCatOut=[]
-    #         for g in range(len(phaseTest)):
-    #             tempPeriodCatOut.append([phaseTest[g],varData[g,1]])
-    #         tempPeriodCatOut=asarray(tempPeriodCatOut)
-    #         savetxt(periodPath / f"{variableName}_StringTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
-    
-    #         tempPeriodCatOut=[]
-    #         for g in range(len(varData[:,0])):
-    #             tempPeriodCatOut.append([(varData[g,0]/(pdm["distance_minperiod"]) % 1), varData[g,1], varData[g,2]])
-    #         tempPeriodCatOut=asarray(tempPeriodCatOut)
-    #         savetxt(periodPath / f"{variableName}_String_PhasedDiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
-
-    #     if np.isnan(pdm["stdev_results"][0]) or pdm["stdev_results"][0] == 0.0 or (varData.size < 4):
-    #         logger.info("No PDM results due to lack of datapoint coverage")
-    #     else:
-    #         logger.debug("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"]))
-    #         phaseTest=(varData[:,0] / (pdm["stdev_minperiod"])) % 1
-    #         logger.debug("PDM method error: " + str(pdm["stdev_error"]))
-
-    #         with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:
-    #             f.write("PDM Method Estimate (days): "+ str(pdm["stdev_minperiod"])+"\n")
-    #             f.write("PDM method error: " + str(pdm["stdev_error"])+"\n\n")
-
-
-    #     plt.plot(pdm["periodguess_array"], pdm["stdev_results"])
-    #     plt.gca().invert_yaxis()
-    #     plt.title("Range {0} d  Steps: {1}".format(trialRange, periodsteps))
-    #     plt.xlabel(r"Trial Period")
-    #     plt.ylabel(r"Likelihood of Period")
-    #     plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot.png")
-
-    #     plt.clf()
-
-    #     if (varData.size > 3):
-    #         plt.plot(phaseTest, varData[:,1], 'bo', linestyle='None')
-    #         plt.plot(phaseTest+1, varData[:,1], 'ro', linestyle='None')
-    #         plt.errorbar(phaseTest, varData[:,1], yerr=varData[:,2], linestyle='None')
-    #         plt.errorbar(phaseTest+1, varData[:,1], yerr=varData[:,2], linestyle='None')
-    #         plt.gca().invert_yaxis()
-    #         plt.title("Period: {0} d  Steps: {1}".format(pdm["stdev_minperiod"], periodsteps))
-    #         plt.xlabel(r"Phase ($\phi$)")
-    #         plt.ylabel(r"Differential " + str(filterCode) + " Magnitude")
-    #         plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot.png")
-    #         plt.clf()
-
-    #     if calibFile.exists():
-    #         if (calibData.size > 3):
-    #             phaseTestCalib=(calibData[:,0] / (pdm["stdev_minperiod"])) % 1
-    #             plt.plot(phaseTestCalib, calibData[:,1], 'bo', linestyle='None')
-    #             plt.plot(phaseTestCalib+1, calibData[:,1], 'ro', linestyle='None')
-    #             plt.errorbar(phaseTestCalib, calibData[:,1], yerr=calibData[:,2], linestyle='None')
-    #             plt.errorbar(phaseTestCalib+1, calibData[:,1], yerr=calibData[:,2], linestyle='None')
-    #             plt.gca().invert_yaxis()
-    #             plt.title("Period: {0} d  Steps: {1}".format(pdm["stdev_minperiod"], periodsteps))
-    #             plt.xlabel(r"Phase ($\phi$)")
-    #             plt.ylabel(r"Calibrated " + str(filterCode) + " Magnitude")
-    #             plt.savefig(periodPath / f"{variableName}_PDMTestPeriodPlot_Calibrated.png")
-    #             plt.clf()
-    
-    #             tempPeriodCatOut=[]
-    #             for g in range(len(calibData[:,0])):
-    #                 tempPeriodCatOut.append([(calibData[g,0]/(pdm["stdev_minperiod"])) % 1, calibData[g,1], calibData[g,2]])
-    #             tempPeriodCatOut=asarray(tempPeriodCatOut)
-    #             savetxt(periodPath / f"{variableName}_PDM_PhasedCalibMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
-
-    #     if (varData.size > 3):
-
-    #         tempPeriodCatOut=[]
-    #         for g in range(len(phaseTest)):
-    #             tempPeriodCatOut.append([phaseTest[g],varData[g,1]])
-    #         tempPeriodCatOut=asarray(tempPeriodCatOut)
-    #         savetxt(periodPath / f"{variableName}_PDMTrial.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
-    
-    #         tempPeriodCatOut=[]
-    #         for g in range(len(varData[:,0])):
-    #             tempPeriodCatOut.append([(varData[g,0]/(pdm["stdev_minperiod"])) % 1, varData[g,1], varData[g,2]])
-    #         tempPeriodCatOut=asarray(tempPeriodCatOut)
-    #         savetxt(periodPath / f"{variableName}_PDM_PhaseddiffMags.csv", tempPeriodCatOut, delimiter=",", fmt='%0.8f')
-
-    #     # Plot publication plots
-
-    #     plt.figure(figsize=(5, 3))
+    return 0.0
 
-    #     plt.plot(pdm["periodguess_array"], pdm["stdev_results"], linewidth=0.5)
-    #     plt.gca().invert_yaxis()
-    #     plt.xlabel(r"Trial Period")
-    #     plt.ylabel(r"Likelihood of Period")
-    #     plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4)
-    #     plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.png", dpi=300)
-    #     plt.savefig(periodPath / f"{variableName}_PDMLikelihoodPlot_Publication.eps")
-
-    #     plt.clf()
-
-    #     plt.figure(figsize=(5, 3))
-
-    #     plt.plot(pdm["periodguess_array"], pdm["distance_results"], linewidth=0.5)
-    #     plt.gca().invert_yaxis()
-    #     plt.xlabel(r"Trial Period")
-    #     plt.ylabel(r"Likelihood of Period")
-    #     plt.subplots_adjust(left=0.15, right=0.99, top=0.98, bottom=0.15, wspace=0.3, hspace=0.4)
-    #     plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.png", dpi=300)
-    #     plt.savefig(periodPath / f"{variableName}_StringLikelihoodPlot_Publication.eps")
-
-    #     plt.clf()
-
-
-    #     # ANOVA
-        
-    #     if calibFile.exists():
-    #         if (calibData.size > 3):
-    #             if len(calibData[:,0]) < 75:
-    #                 binsize=0.1
-    #             else:
-    #                 binsize=0.05
-    #             minperbin=int((len(calibData[:,0])/10))
-    #         elif (varData.size > 3):
-    #             if len(varData[:,0]) < 75:
-    #                 binsize=0.1
-    #             else:
-    #                 binsize=0.05
-    #             minperbin=int((len(varData[:,0])/10))
-                
-    #     else:
-    #         if (varData.size > 3):
-    #             if len(varData[:,0]) < 75:
-    #                 binsize=0.1
-    #             else:
-    #                 binsize=0.05
-    #             minperbin=int((len(varData[:,0])/10))
-        
-    #     if 'minperbin' in locals():
-    #         if minperbin > 10:
-    #             minperbin=10
-    #     else:
-    #         minperbin = 3
-            
-        
-    #     # Theta Anova Method off for the moment until I put in a command-line option
-        
-    #     # if calibFile.exists():
-    #     #     if (calibData.size > 3):
-    #     #         aovoutput=aov_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-    #     # else:
-    #     #     if (varData.size > 3):
-    #     #         aovoutput=aov_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, phasebinsize=binsize, mindetperbin=minperbin, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-
-
-    #     # logger.debug("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"]))
-    #     #with open(paths['parent'] / "periodEstimates.txt", "a+") as f:            
-    #     #    f.write("Theta Anova Method Estimate (days): " + str(aovoutput["bestperiod"])+"\n")
-        
-    #     if calibFile.exists():
-    #         if (calibData.size > 3):
-    #             aovhmoutput=aovhm_periodfind((calibData[:,0]),(calibData[:,1]),(calibData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-    #         elif (varData.size > 3):                
-    #             aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-                
-    #     else:
-    #         if (varData.size > 3):
-    #             aovhmoutput=aovhm_periodfind((varData[:,0]),(varData[:,1]),(varData[:,2]), sigclip=False, autofreq=False, startp=minperiod, endp=maxperiod, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-            
-    #     logger.debug("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"]))
-    #     with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:            
-    #         f.write("Harmonic Anova Method Estimate (days): " + str(aovhmoutput["bestperiod"])+"\n")
-
-    #     # LOMB SCARGLE
-    #     for nts in range(2):
-    #         if calibFile.exists():
-    #             if (calibData.size > 3):
-    #                 lscargoutput = LombScargleMultiterm('periodifile', (calibData[:, 0]), (calibData[:, 1]), (calibData[:, 2]),
-    #                                                     nterms=nts+1,
-    #                                                     periodlower=minperiod, periodupper=maxperiod, samples=20,
-    #                                                     disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-        
-    #                 logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput))
-    #                 with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:            
-    #                     f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n")
-    #         else:
-    #             if (varData.size > 3):
-    #                 lscargoutput = LombScargleMultiterm('periodifile', (varData[:, 0]), (varData[:, 1]), (varData[:, 2]),
-    #                                                     nterms=nts+1,
-    #                                                     periodlower=minperiod, periodupper=maxperiod, samples=20,
-    #                                                     disablelightcurve=False, periodPath=periodPath, variableName=variableName, periodsteps=periodsteps)
-        
-    #                 logger.debug('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput))
-    #                 with open(paths['parent'] / "results/periodEstimates.txt", "a+") as f:            
-    #                     f.write('Lomb-Scargle N=' + str(nts+1) + ' Period Best Estimate: ' + str(lscargoutput)+"\n")
-
-    if 'pdm' in locals():
-        return pdm["distance_minperiod"]
-    else:
-        return 0.0
+# Needed for windows to multiprocess appropriately
+if __name__ == "__main__":
+    pass
\ No newline at end of file
diff --git a/astrosource/plots.py b/astrosource/plots.py
index 686c494..f55011f 100644
--- a/astrosource/plots.py
+++ b/astrosource/plots.py
@@ -9,8 +9,8 @@
 from astropy.units import degree, arcsecond
 import logging
 import numpy as np
+import traceback
 
-#from astrosource.utils import photometry_files_to_array, AstrosourceException
 from astrosource.utils import AstrosourceException
 
 logger = logging.getLogger('astrosource')
@@ -94,17 +94,20 @@ def plot_variability(output, variableID, parentPath, compFile):
         if (parentPath / 'results/calibCompsUsed.csv').exists():
             logger.debug("Calibrated")
             calibCompFile=genfromtxt(parentPath / 'results/calibCompsUsed.csv', dtype=float, delimiter=',')
-            
-            if len(calibCompFile) == 5:
-                calibCompSkyCoord = SkyCoord(calibCompFile[0],calibCompFile[1], frame='icrs', unit=degree)
-                calibnumber=1
-            else:
-                calibCompSkyCoord = SkyCoord(calibCompFile[:,0],calibCompFile[:,1], frame='icrs', unit=degree)
-                calibnumber=len(calibCompSkyCoord)
+            try:
+                if len(calibCompFile) == 5 and calibCompFile.size == 5:
+                    calibCompSkyCoord = SkyCoord(calibCompFile[0],calibCompFile[1], frame='icrs', unit=degree)
+                    calibnumber=1
+                else:
+                    calibCompSkyCoord = SkyCoord(calibCompFile[:,0],calibCompFile[:,1], frame='icrs', unit=degree)
+                    calibnumber=len(calibCompSkyCoord)
+            except:
+                print ("MTF FIXING THIS ERROR")
+                logger.error(traceback.print_exc())
+                breakpoint()
             calibCompExist=True
             
             calibCompStarPlot = []        
-            #for q in range(len(calibCompSkyCoord)):
             if calibnumber==1:
                 idx, d2d, _ = calibCompSkyCoord.match_to_catalog_sky(outputSkyCoord)
                 calibCompStarPlot.append([output[idx][2],output[idx][3]])
@@ -135,7 +138,6 @@ def plot_variability(output, variableID, parentPath, compFile):
         plt.plot(np.asarray(compStarPlot)[:,0],np.asarray(compStarPlot)[:,1], 'yo')
         if calibCompExist:
             plt.plot(np.asarray(calibCompStarPlot)[:,0],np.asarray(calibCompStarPlot)[:,1], 'ro', mfc='none')
-        # plt.plot(linex, liney)
 
         plt.ylim(max([min(outploty)-0.04,0.0]), max(outploty)+0.04)
         plt.xlim(min(outplotx)-0.1, max(outplotx)+0.1)
@@ -190,7 +192,6 @@ def plot_variability(output, variableID, parentPath, compFile):
         plt.plot(np.asarray(compStarPlot)[:,0],np.asarray(compStarPlot)[:,1], 'yo')
         if calibCompExist:
             plt.plot(np.asarray(calibCompStarPlot)[:,0],np.asarray(calibCompStarPlot)[:,1], 'ro', mfc='none')
-        # plt.plot(linex, liney)
         plt.ylim(max([min(outploty)-0.04,0.0]), max(outploty)+0.04)
         plt.xlim(min(outplotx)-0.1, max(outplotx)+0.1)
         plt.grid(True)
@@ -210,7 +211,6 @@ def plot_variability(output, variableID, parentPath, compFile):
         if calibCompExist:
             plt.plot(np.asarray(calibCompStarPlot)[:,0],np.asarray(calibCompStarPlot)[:,1], 'ro', mfc='none')
         fig.set_size_inches(16,9)
-        # plt.plot(linex, liney)
         plt.ylim(max([min(outploty)-0.04,0.0]), max(outploty)+0.04)
         plt.xlim(min(outplotx)-0.1, max(outplotx)+0.1)
         plt.grid(True)
@@ -368,3 +368,7 @@ def phased_plots(paths, filterCode, targets, period, phaseShift):
 
 
     return
+
+# Needed for windows to multiprocess appropriately
+if __name__ == "__main__":
+    pass
diff --git a/astrosource/utils.py b/astrosource/utils.py
index 33b1852..3cf5463 100644
--- a/astrosource/utils.py
+++ b/astrosource/utils.py
@@ -95,7 +95,7 @@ def get_targets(targetfile):
     # Remove any nan rows from targets
     targetRejecter=[]
 
-    if not (targets.shape[0] == 4 and targets.size == 4):
+    if not (targets.shape[0] == 5 and targets.size == 5):
         for z in range(targets.shape[0]):
           if isnan(targets[z][0]):
             targetRejecter.append(z)