83 lines
2.3 KiB
Python
83 lines
2.3 KiB
Python
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)
|
|
|
|
|
|
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)
|
|
|
|
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)
|
|
plt.scatter(oneOverDensities, breakupDistances3)
|
|
ax.plot(oneOverDensities, bestFit)
|
|
ax.plot(oneOverTheoreticalDensities, theoreticalBreakupDistances3)
|
|
plt.show()
|
|
|