Add perpendicular distance option

CIVfunctions.py

@@ -1,25 +1,30 @@
 import numpy as np
 import joblib
 
-def project(data, fit):
+def project(data, fit, perp_dist=False):
     #Use: Project scattered data (data) onto a line of best fit (fit)
-    #Returns: 2-D array of (x,y) locations on the line of your orthogonal projection
+    #Returns: 2-D array of (x,y) locations on the line of your orthogonal projection (and optionally the orthogonal distances from said line)
 
     #data: 2-D (N by 2) array of your data in some x-y space
     #fit: 2-D (N by 2) array of coordinates along your line of best fit
-
+    #perp_dist: Bool; if True, also return the data points' perpendicular distances from the best-fit curve    
+
     locs = [] #this will contain xy locations along the best-fit curve corresponding to each data point's orthogonal projection onto said curve - python list faster than numpy array
-
+    perp_distances = [] #to save perp distances from curve in    
+
     #Loop through each data point and compare with locations along fit
     for scat in data: 
         r2 = (scat[0]-fit[:,0])**2 + (scat[1]-fit[:,1])**2      #dist^2 of scat from each point along fit
         delta = [fit[np.argmin(r2), 0], fit[np.argmin(r2), 1]]  #save point along fit where dist^2 was minimum
-        locs.append(delta)  
+        locs.append(delta)
+        perp_distances.append(np.sqrt(r2.min()))  
 
+    if perp_dist:
+        return np.array(locs), np.array(perp_distances)
     return np.array(locs)
 
 
-def CIV_distance(data, fit, step=1, logEW=True, path="./"):
+def CIV_distance(data, fit, step=1, logEW=True, path="./", perp_dist=False):
     #fit: N-by-2 array containing coordinates of points along best fit line 
     #data: N-by-2 [[x,y]] array of data 
     #NOTE: This really just caters to this situation (assumes monotonically decreasing fit_original)
@@ -36,9 +41,13 @@ def CIV_distance(data, fit, step=1, logEW=True, path="./"):
 
     #2) data is now each point's orthogonal projection onto fit
     ind_projcut = -((fit.shape[0]-1)%step)
-    if ind_projcut!=0: data = project(data, fit[:ind_projcut])
-    else: data = project(data, fit)
-
+    if not perp_dist:
+        if ind_projcut!=0: data = project(data, fit[:ind_projcut])
+        else: data = project(data, fit)
+    else:
+        if ind_projcut!=0: data, perpendicular_dist = project(data, fit[:ind_projcut], perp_dist=True)
+        else: data, perpendicular_dist = project(data, fit, perp_dist=True)
+
     darr = [] #list to fill with distances along best-fit line for each point
 
     #3) Loop through each data point- start at tip of line and sum dist traveled until passing data point
@@ -52,5 +61,7 @@ def CIV_distance(data, fit, step=1, logEW=True, path="./"):
             if yp >= scat[1] >= y: #if we pass the projected y-coord, save the distance traveled
                 darr.append((d+dp)/2) #save the average between dprevious and d
                 break
-
+
+    if perp_dist:
+        return np.array(darr), perpendicular_dist
     return np.array(darr)