diff --git a/distanceVDensity.py b/allSecondDerivs.py
similarity index 70%
rename from distanceVDensity.py
rename to allSecondDerivs.py
index 918ee85..4d8c2e1 100644
--- a/distanceVDensity.py
+++ b/allSecondDerivs.py
@@ -6,12 +6,13 @@ import matplotlib.pyplot as plt
 columns = ["ID", "idx", "Mass", "Radius", "X", "Y", "Z", "vX", "vY", "vZ", "sX", "sY", "sZ", "Colour"]
 deviations = []
 distances = []
-time = []
 densities = ["300","350","400","450","500","550","600","650","700"]
 breakUpDistances = []
 
 for density in densities:
-    for i in range (1,320):
+    densityDeviations = []
+    time = []
+    for i in range (25,270):
         num = str(i).rjust(5, '0')
         file = "BTs/High-Res-"+density+"/boom."+num+".bt"
         data = pd.read_csv(file, sep=' ', names=columns)
@@ -26,24 +27,26 @@ for density in densities:
         deviation = data["distance"].std()
         distance = data["distance"].mean()
         time.append(i)
-        deviations.append(deviation)
+        densityDeviations.append(deviation)
         distances.append(distance)
 
 
     del time[-1]
     del time[-1]
-    del distances[-1]
-    del distances[-1]
-    firstDeriv=np.diff(deviations)
+    firstDeriv=np.diff(densityDeviations)
     secondDeriv=np.diff(firstDeriv)
-
-    max_index=np.argmax(secondDeriv)
-
-    print(distances[max_index])
-    breakUpDistances.append(distances[max_index])
+    deviations.append(secondDeriv)
 
 ax=plt.axes()
 # ax.set_yscale("log")
-ax.plot(densities, breakUpDistances)
+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])
 plt.show()
 
diff --git a/distanceVdensity.py b/distanceVdensity.py
new file mode 100644
index 0000000..a57a752
--- /dev/null
+++ b/distanceVdensity.py
@@ -0,0 +1,64 @@
+import math
+import pandas as pd
+import numpy as np
+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
+    theoreticalDistances.append(distance)
+
+for density in densities:
+    densityDeviations = []
+    time = []
+    for i in range (25,270):
+        num = str(i).rjust(5, '0')
+        file = "BTs/High-Res-"+str(density)+"/boom."+num+".bt"
+        data = pd.read_csv(file, sep=' ', names=columns)
+        data["X"] = data["X"]*1.5e8
+        data["Y"] = data["Y"]*1.5e8
+        data["Z"] = data["Z"]*1.5e8
+        data["X2"] = data["X"]**2
+        data["Y2"] = data["Y"]**2
+        data["Z2"] = data["Z"]**2
+        data["distance2"] = data["X2"]+data["Y2"]+data["Z2"]
+        data["distance"] = data["distance2"]**0.5
+        deviation = data["distance"].std()
+        distance = data["distance"].mean()
+        time.append(i)
+        densityDeviations.append(deviation)
+        distances.append(distance)
+
+
+    del time[-1]
+    del time[-1]
+    firstDeriv=np.diff(densityDeviations)
+    secondDeriv=np.diff(firstDeriv)
+
+    maxIndex=np.argmax(secondDeriv)
+    breakupDistance=distances[maxIndex]
+    breakupDistances.append(breakupDistance)
+
+# oneOverDensities = []
+# for i in densities:
+#     oneOver = 1/i
+#     oneOverDensities.append(oneOver)
+
+# breakupDistances3 = []
+# for i in breakupDistances:
+#     cubed = i**3
+#     breakupDistances3.append(cubed)
+    
+ax=plt.axes()
+# ax.set_yscale("log")
+ax.plot(densities, breakupDistances)
+ax.plot(theoreticalDensities, theoreticalDistances)
+plt.show()
+
diff --git a/energy.py b/energy.py
index 9dfbf71..d9538bf 100644
--- a/energy.py
+++ b/energy.py
@@ -12,6 +12,7 @@ args=parser.parse_args()
 columns = ["ID", "idx", "Mass", "Radius", "X", "Y", "Z", "vX", "vY", "vZ", "sX", "sY", "sZ", "Colour"]
 Kenergies = []
 Genergies = []
+energies = []
 time = []
 
 for i in range (1,401):
@@ -37,8 +38,9 @@ for i in range (1,401):
     data["vel"] = data["distance2"]**0.5
     data["KE"] = 0.5*data["Mass"]*data["vel2"]
     data["GPE"] = (6.67e-11*1.89e27*data["Mass"])/(data["distance"])
-    # data["Energy"] = data["KE"]+data["GPE"]
-    # energy = data["Energy"].sum()
+    data["Energy"] = data["KE"]+data["GPE"]
+    energy = data["Energy"].sum()
+    energies.append(energy)
     time.append(i)
     KE = data["KE"].sum()
     Kenergies.append(KE)
@@ -49,5 +51,6 @@ ax=plt.axes()
 # ax.set_yscale("log")
 ax.plot(time, Kenergies, label="Kinetic")
 ax.plot(time, Genergies, label="Gravitational")
+ax.plot(time, energies, label="Total")
 ax.legend()
 plt.show()