From 6c5f6b8459b17aa2316902a7fa5c9e354ec854b3 Mon Sep 17 00:00:00 2001
From: Anson <anson@ansonbiggs.com>
Date: Sat, 20 Mar 2021 15:15:37 -0700
Subject: [PATCH] added monte carlo simulation

---
 thrust.jl | 54 ++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 34 insertions(+), 20 deletions(-)

diff --git a/thrust.jl b/thrust.jl
index 4284e77..c4841e9 100644
--- a/thrust.jl
+++ b/thrust.jl
@@ -1,37 +1,42 @@
 using Unitful
 using DataFrames
 using Plots
+# using Gadfly
 using UnitfulRecipes
 using Roots
 using CSV
+using MonteCarloMeasurements
 
 
+unsafe_comparisons(true)
+
 # Tank https://www.amazon.com/Empire-Paintball-BASICS-Pressure-Compressed/dp/B07B6M48SR/
-V = 90u"inch^3"
-P0 = 4500.0u"psi"
-Wtank = 2.3u"lb"
-Pmax = 100u"psi"
+V = (85 ± 5)u"inch^3"
+P0 = (4200.0 ± 300)u"psi"
+Wtank = (2.3 ± 0.2)u"lb"
+Pmax = (100 ± 5)u"psi"
 
 Wsolenoid = 1.5u"kg"
 
 # Params
-d_nozzle = (1 // 16) * u"inch"
+d_nozzle = ((1 // 16) ± 0.001)u"inch"
 a_nozzle = (pi / 4) * d_nozzle^2
-Poutmax = 800u"psi"
+Poutmax = (800 ± 50)u"psi"
 
 # Universal Stuff
-P_amb = 1u"atm"
-γ = 1.4
+P_amb = (1 ± 0.2)u"atm"
+γ = 1.4 ± 0.05
 R = 287.05u"J/(kg * K)"
-T = 300u"K"
+T = (300 ± 20)u"K"
 
 # Setting up while loop
-t = 1.0u"s"
+t = 0.0u"s"
 P = P0 |> u"Pa"
 M = V * (P / (R * T)) |> u"kg"
-df = DataFrame(Thrust = 0.0u"N", Pressure = P, Time = t, Mass = M)
-ts = 0.001u"s"
-while M > 0.01u"kg"
+# df = DataFrame(Thrust = (0 ± 0)u"N", Pressure = P, Time = t, Mass = M)
+df = DataFrame()
+ts = 10u"ms"
+while M > 0.1u"kg"
     # while t < 30u"s"
     # Calculate what is leaving tank
     P = minimum([P, Pmax])
@@ -44,16 +49,25 @@ while M > 0.01u"kg"
     # Calculate what is still in the tank
     M = M - ṁ * ts |> u"kg"
     P = (M * R * T) / V |> u"Pa"
-    t = t + ts
+    # df_step = DataFrame(Thrust = Thrust, Pressure = P, Time = t, Mass = M)
 
-    df_step = DataFrame(Thrust = Thrust, Pressure = P, Time = t, Mass = M)
-    append!(df, df_step)
+    if df == DataFrame()
+        df = DataFrame(Thrust = Thrust, Pressure = P, Time = t, Mass = M)
+    else
+        append!(df, DataFrame(Thrust = Thrust, Pressure = P, Time = t, Mass = M))
+    end
+    t = t + ts
 end
+impulse = sum(df.Thrust) * ts |> u"N*s"
 println("---------------------------\n\n\n\n")
-println("Total Impulse: ", sum(df.Thrust) * ts)
+println("Total Impulse: ", impulse)
 println("Burn Time: ", t)
+plot(impulse)
 println("Mass Total: ", Wtank + Wsolenoid + maximum(df.Mass))
 
-print(describe(df))
-plot(df.Time, df.Thrust, title = "Thrust Over Time")
-CSV.write("thrustdata.csv", df .|> ustrip)
+# plot(sum(df.Thrust) * ts |> u"N*s")
+# print(describe(df))
+# df.Time = round(u"s", df.Time)
+# plot(df.Time, df.Thrust, title = "Thrust Over Time")
+# CSV.write("thrustdata.csv", df .|> ustrip)
+# plot(df.Time .|> u"s" .|> ustrip, df.Thrust .|> u"N" .|> ustrip)> ustrip, df.Thrust .|> u"N" .|> ustrip)
\ No newline at end of file