using Unitful
using DataFrames
using Plots
theme(:ggplot2);
using CSV
using Measurements
using Measurements: value, uncertainty
using Printf

# Tank https://www.amazon.com/Empire-Paintball-BASICS-Pressure-Compressed/dp/B07B6M48SR/
V = (85 ± 5)u"inch^3"
P0 = (4200.0 ± 300)u"psi"
Wtank = (2.3 ± 0.2)u"lb"
Pmax = (250 ± 50)u"psi" # Max Pressure that can come out the nozzle

Wsolenoid = 1.5u"kg"

# Params
d_nozzle = ((1 // 18) ± 0.001)u"inch"
a_nozzle = (pi / 4) * d_nozzle^2

# Universal Stuff
P_amb = (1 ± 0.2)u"atm"
γ = 1.4 ± 0.05
R = 287.05u"J/(kg * K)"
T = (300 ± 20)u"K"

t = 0.0u"s"
P = P0 |> u"Pa"
M = V * (P / (R * T)) |> u"kg"
ts = 1u"ms"
df = DataFrame(Thrust = (0 ± 0)u"N", Pressure = P, Time = t, Mass = M)

while M > 0.005u"kg"
    # Calculate what is leaving tank
    P = minimum([P, Pmax])
    ve = sqrt((2 * γ / (γ - 1)) * R * T * (1 - P_amb / P)^((γ - 1) / γ)) |> u"m/s"
    ρ = P / (R * T) |> u"kg/m^3"
    ṁ = ρ * a_nozzle * ve |> u"kg/s"

    Thrust = ṁ * ve + a_nozzle * (P - P_amb) |> u"N"

    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)
    append!(df, df_step)
end
final_time = t |> u"s"
impulse = sum(df.Thrust) * ts |> u"N*s"
println("---------------------------\n\n\n\n")
println("Burn Time: ", t)

println("Mass Total: ", Wtank + Wsolenoid + maximum(df.Mass))

print(describe(df))

thrust_values = df.Thrust .|> ustrip .|> value
thrust_uncertainties = df.Thrust .|> ustrip .|> uncertainty

plot(
    df.Time .|> ustrip,
    thrust_values,
    title = "Thrust Over Time",
    ribbon = (thrust_uncertainties, thrust_uncertainties),
    fillalpha = 0.2,
    label = "Thrust",
    xlabel = "Time (s)",
    ylabel = "Thrust (N)",
)

out = DataFrame(
    Thrust = thrust_values,
    Uncertainty = thrust_uncertainties,
    Time = df.Time .|> u"s" .|> ustrip,
)
CSV.write("AirProp.csv", out)

### Save data to readme.md

savefig("ThrustCurve.png")

include("prop_compare.jl")

b = IOBuffer();
t = TextDisplay(b);
display(t, "text/html", describe(df));
table = String(take!(b)); # https://stackoverflow.com/a/60443621/8774114

t = TextDisplay(b);
display(t, "text/html", describe(df));
table = String(take!(b)); # https://stackoverflow.com/a/60443621/8774114

readme = Printf.@sprintf """
# Air Prop Thrust 

`readme auto generated, run air_prop_sim.jl to recreate.`

![Thrust Over Time](ThrustCurve.png)

Total Impulse: %s

Burn TIme: %s

More Data:

%s

## Comparison to Rocket Motors

![Propulsion Comparison](prop_compare.png)

https://www.thrustcurve.org/

""" impulse final_time table

write("readme.md", readme);