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 = []
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
    theoreticalDistances.append(distance)

for density in densities:
    densityDeviations = []
    time = []    
    distances = []
    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)

    for i in range(len(densityDeviations)):
        if densityDeviations[i] > 0.75:
            index = i
            break

    breakupDistances.append(distances[index])


    
ax=plt.axes()
plt.scatter(densities, breakupDistances)
ax.plot(theoreticalDensities, theoreticalDistances)
plt.show()