diff --git a/allSecondDerivs.py b/allSecondDerivs.py
index 4d8c2e1..ae5ff69 100644
--- a/allSecondDerivs.py
+++ b/allSecondDerivs.py
@@ -39,14 +39,16 @@ for density in densities:
 
 ax=plt.axes()
 # ax.set_yscale("log")
-ax.plot(time, deviations[0])
-ax.plot(time, deviations[1])
-ax.plot(time, deviations[2])
-ax.plot(time, deviations[3])
-ax.plot(time, deviations[4])
-ax.plot(time, deviations[5])
-ax.plot(time, deviations[6])
-ax.plot(time, deviations[7])
-ax.plot(time, deviations[8])
+ax.plot(time, deviations[0],label="300")
+ax.plot(time, deviations[1],label="350")
+ax.plot(time, deviations[2],label="400")
+ax.plot(time, deviations[3],label="450")
+ax.plot(time, deviations[4],label="500")
+ax.plot(time, deviations[5],label="550")
+ax.plot(time, deviations[6],label="600")
+ax.plot(time, deviations[7],label="650")
+ax.plot(time, deviations[8],label="700")
+ax.legend()
 plt.show()
 
+
diff --git a/distanceVdensity.py b/distanceVdensity.py
index a57a752..5071ddd 100644
--- a/distanceVdensity.py
+++ b/distanceVdensity.py
@@ -5,19 +5,19 @@ import matplotlib.pyplot as plt
 
 columns = ["ID", "idx", "Mass", "Radius", "X", "Y", "Z", "vX", "vY", "vZ", "sX", "sY", "sZ", "Colour"]
 deviations = []
-distances = []
 densities = [300,350,400,450,500,550,600,650,700]
 breakupDistances = []
 theoreticalDensities=np.linspace(300,700,num=100)
 theoreticalDistances=[]
 for i in theoreticalDensities:
     distance=2.44*69911*(1330/i)**(1/3)
-    distance=distance/2
+    distance=distance
     theoreticalDistances.append(distance)
 
 for density in densities:
     densityDeviations = []
-    time = []
+    time = []    
+    distances = []
     for i in range (25,270):
         num = str(i).rjust(5, '0')
         file = "BTs/High-Res-"+str(density)+"/boom."+num+".bt"
@@ -46,19 +46,38 @@ for density in densities:
     breakupDistance=distances[maxIndex]
     breakupDistances.append(breakupDistance)
 
-# oneOverDensities = []
-# for i in densities:
-#     oneOver = 1/i
-#     oneOverDensities.append(oneOver)
+oneOverDensities = []
+for i in densities:
+    oneOver = 1/i
+    oneOverDensities.append(oneOver)
 
-# breakupDistances3 = []
-# for i in breakupDistances:
-#     cubed = i**3
-#     breakupDistances3.append(cubed)
+breakupDistances3 = []
+for i in breakupDistances:
+    cubed = i**3
+    breakupDistances3.append(cubed)
+    
+oneOverTheoreticalDensities = []
+for i in theoreticalDensities:
+    oneOver = 1/i
+    oneOverTheoreticalDensities.append(oneOver)
+
+theoreticalBreakupDistances3 = []
+for i in theoreticalDistances:
+    cubed = i**3
+    theoreticalBreakupDistances3.append(cubed)
+
+a,b = np.polyfit(oneOverDensities, breakupDistances3, 1)
+
+bestFit = []
+for i in oneOverDensities:
+    y = a*i+b
+    bestFit.append(y)
     
 ax=plt.axes()
 # ax.set_yscale("log")
-ax.plot(densities, breakupDistances)
-ax.plot(theoreticalDensities, theoreticalDistances)
+# ax.plot(densities, breakupDistances)
+plt.scatter(oneOverDensities, breakupDistances3)
+ax.plot(oneOverDensities, bestFit)
+ax.plot(oneOverTheoreticalDensities, theoreticalBreakupDistances3)
 plt.show()