From 1b2da80b101f0490b38d6e32a1120642f8a1fad1 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Wed, 7 Feb 2024 00:57:19 -0500 Subject: do hw2 --- hw2/2-18.jl | 91 +++++++++++++++++++++++ hw2/2-18a.svg | 49 +++++++++++++ hw2/2-18b.png | Bin 0 -> 26182 bytes hw2/2-7.jl | 207 ++++++++++++++++++++++++++++++++++++++++++++++++++++ hw2/2-7a.svg | 51 +++++++++++++ hw2/2-7b.png | Bin 0 -> 45828 bytes hw2/2-7c.svg | 47 ++++++++++++ hw2/2-7d.png | Bin 0 -> 45294 bytes hw2/2-7e.svg | 49 +++++++++++++ hw2/2-7f.png | Bin 0 -> 45817 bytes hw2/2-9.jl | 136 ++++++++++++++++++++++++++++++++++ hw2/2-9a.svg | 55 ++++++++++++++ hw2/hw2-writeup.pdf | Bin 0 -> 763541 bytes 13 files changed, 685 insertions(+) create mode 100644 hw2/2-18.jl create mode 100644 hw2/2-18a.svg create mode 100644 hw2/2-18b.png create mode 100644 hw2/2-7.jl create mode 100644 hw2/2-7a.svg create mode 100644 hw2/2-7b.png create mode 100644 hw2/2-7c.svg create mode 100644 hw2/2-7d.png create mode 100644 hw2/2-7e.svg create mode 100644 hw2/2-7f.png create mode 100644 hw2/2-9.jl create mode 100644 hw2/2-9a.svg create mode 100644 hw2/hw2-writeup.pdf diff --git a/hw2/2-18.jl b/hw2/2-18.jl new file mode 100644 index 0000000..58bad47 --- /dev/null +++ b/hw2/2-18.jl @@ -0,0 +1,91 @@ +using Plots # for plotting trajectory + +# simulation parameters +Δt = 0.001 # time step +x_max = 60.0 # in feet +v_0 = 95.0 # mph +v_0 = v_0 * 5280 / 3600 # convert mph to feet/s +y_0 = 6.0 # in feet +ω = 0.0 # in rad/s +m = 149 # in grams +S_0_over_m = 4.1 * 10^(-4) # unitless +B_ref_over_m = 4.0 * 10^(-5) # in m-1, at 300K +g = 32.2 # in ft/s^2 + +function dynamics!( + x::Vector{Float64}, y::Vector{Float64}, + v_y::Vector{Float64}, v_x::Vector{Float64}, + t::Vector{Float64}) + while x[end] <= x_max + # decompose previous positions and velocities + x_i = x[end] + y_i = y[end] + v_x_i = v_x[end] + v_y_i = v_y[end] + + # calculate new positions + x_new = x_i + v_x_i * Δt + y_new = y_i + v_y_i * Δt + + # calculate drag force + v_i = sqrt(v_x_i^2 + v_y_i^2) + F_x = - B_ref_over_m * v_x_i * v_i + F_y = - g + S_0_over_m * v_x_i * ω + + # calculate new velocities + v_x_new = v_x_i + F_x * Δt + v_y_new = v_y_i + F_y * Δt + + # store new positions and velocities + push!(x, x_new) + push!(y, y_new) + push!(v_x, v_x_new) + push!(v_y, v_y_new) + push!(t, t[end] + Δt) + end +end + +# interpolate the last point that's past the x_max +function interpolate!(x::Vector{Float64}, y::Vector{Float64}) + if x[end] <= x_max + return + end + + x_i = x[end-1] + y_i = y[end-1] + + x_f = x[end] + y_f = y[end] + + m = (y_f - y_i) / (x_f - x_i) + b = y_i - m * x_i + + x_new = x_max + y_new = m * x_new + b + + x[end] = x_new + y[end] = y_new +end + +# run the simulation with ω = 0 +x = [0.0] +y = [y_0] +v_x = [v_0] # pitch only in x direction +v_y = [0.0] +t = [0.0] +dynamics!(x, y, v_y, v_x, t) +interpolate!(x, y) +plot(x, y, label="ω = 0", xlabel="x (feet)", ylabel="y (feet)", lw=2, title="Trajectory of a Fastball with Backspin") +println("y-value at x_max (60ft) for ω = 0:\t", y[end], " feet") + +# run the simulation with ω = 1000 rpm +ω = 1000 * 2 * π / 60 # convert rpm to rad/s +x = [0.0] +y = [y_0] +v_x = [v_0] # pitch only in x direction +v_y = [0.0] +t = [0.0] +dynamics!(x, y, v_y, v_x, t) +interpolate!(x, y) +println("y-value at x_max (60ft) for ω = 1000 rpm:\t", y[end], " feet") +plot!(x, y, label="ω = 1000 rpm", lw=2) diff --git a/hw2/2-18a.svg b/hw2/2-18a.svg new file mode 100644 index 0000000..26b2e7f --- /dev/null +++ b/hw2/2-18a.svg @@ -0,0 +1,49 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw2/2-18b.png b/hw2/2-18b.png new file mode 100644 index 0000000..2a0c0b8 Binary files /dev/null and b/hw2/2-18b.png differ diff --git a/hw2/2-7.jl b/hw2/2-7.jl new file mode 100644 index 0000000..d936cfb --- /dev/null +++ b/hw2/2-7.jl @@ -0,0 +1,207 @@ +using Plots # for plotting trajectory + +# simulation parameters +Δt = 0.01 # time step +y_min = 0.0 +Δθ = 0.01 # in degrees +θ_start = 0.0 # in degrees +θ_end = 90.0 # in degrees +v_0 = 700.0 # in m/s + +# constants +B_ref_over_m = 4.0 * 10^(-5) # in m-1, at 300K +T_ref = 300.0 # in kelvin +T_0 = 300.0 # in kelvin +g = 9.8 # in m/s^2 + +# isothermic parameters +y_o = 10 * 10^4 # k_BT/mg in meter + +# adiabatic parameters +α = 2.5 # for air +a = 6.5 * 10^(-3) # in kelvin/meter + +function isothermic!( + x::Vector{Float64}, y::Vector{Float64}, + v_y::Vector{Float64}, v_x::Vector{Float64}, + t::Vector{Float64}) + while y[end] >= y_min + # decompose previous positions and velocities + x_i = x[end] + y_i = y[end] + v_x_i = v_x[end] + v_y_i = v_y[end] + + # calculate new positions + x_new = x_i + v_x_i * Δt + y_new = y_i + v_y_i * Δt + + # calculate drag force + v_i = sqrt(v_x_i^2 + v_y_i^2) + F_drag = - B_ref_over_m * (T_0 / T_ref)^α * # temperature variation + ℯ^(-y_i/y_o) # density/altitude variation + F_drag_x = F_drag * v_x_i * v_i + F_drag_y = F_drag * v_y_i * v_i + + # calculate new velocities + v_x_new = v_x_i + F_drag_x * Δt + v_y_new = v_y_i + F_drag_y * Δt - g * Δt + + # store new positions and velocities + push!(x, x_new) + push!(y, y_new) + push!(v_x, v_x_new) + push!(v_y, v_y_new) + push!(t, t[end] + Δt) + end +end + +function adiabatic!( + x::Vector{Float64}, y::Vector{Float64}, + v_y::Vector{Float64}, v_x::Vector{Float64}, + t::Vector{Float64}) + while y[end] >= y_min + # decompose previous positions and velocities + x_i = x[end] + y_i = y[end] + v_x_i = v_x[end] + v_y_i = v_y[end] + + # calculate new positions + x_new = x_i + v_x_i * Δt + y_new = y_i + v_y_i * Δt + + # calculate drag force + v_i = sqrt(v_x_i^2 + v_y_i^2) + F_drag = - B_ref_over_m * (T_0 / T_ref)^α * # temperature variation + (1 - ((a * y_i) / T_0))^α # density/altitude variation + F_drag_x = F_drag * v_x_i * v_i + F_drag_y = F_drag * v_y_i * v_i + + # calculate new velocities + v_x_new = v_x_i + F_drag_x * Δt + v_y_new = v_y_i + F_drag_y * Δt - g * Δt + + # store new positions and velocities + push!(x, x_new) + push!(y, y_new) + push!(v_x, v_x_new) + push!(v_y, v_y_new) + push!(t, t[end] + Δt) + end +end + +function nodensity!( + x::Vector{Float64}, y::Vector{Float64}, + v_y::Vector{Float64}, v_x::Vector{Float64}, + t::Vector{Float64}) + while y[end] >= y_min + # decompose previous positions and velocities + x_i = x[end] + y_i = y[end] + v_x_i = v_x[end] + v_y_i = v_y[end] + + # calculate new positions + x_new = x_i + v_x_i * Δt + y_new = y_i + v_y_i * Δt + + # calculate drag force + v_i = sqrt(v_x_i^2 + v_y_i^2) + F_drag = - B_ref_over_m # coefficient of drag alone + F_drag_x = F_drag * v_x_i * v_i + F_drag_y = F_drag * v_y_i * v_i + + # calculate new velocities + v_x_new = v_x_i + F_drag_x * Δt + v_y_new = v_y_i + F_drag_y * Δt - g * Δt + + # store new positions and velocities + push!(x, x_new) + push!(y, y_new) + push!(v_x, v_x_new) + push!(v_y, v_y_new) + push!(t, t[end] + Δt) + end +end + +# interpolate the last point that's underground +function interpolate!(x::Vector{Float64}, y::Vector{Float64}) + if y[end] == 0 + return # no nothing if y is perfectly on 0 + end + + # calculate x_l, the interpolated x value at y=0 + r = -y[end-1] / y[end] + x_l = (x[end-1] + r * x[end]) / (1 + r) + + # set final values in the array to interpolated point on ground (y=0 ) + x[end] = x_l + y[end] = 0.0 +end + +# arrays for holidng range and angle data +Θ_steps = Int64((θ_end - θ_start) / Δθ) +ranges_adi = zeros(Θ_steps + 1) +angles_adi = zeros(Θ_steps + 1) +ranges_iso = zeros(Θ_steps + 1) +angles_iso = zeros(Θ_steps + 1) +ranges_nod = zeros(Θ_steps + 1) +angles_nod = zeros(Θ_steps + 1) + +for i in 1:Θ_steps+1 + # arrays to store steps + θ = θ_start + (i-1) * Δθ + + local x_adi = [0.0] + local y_adi = [0.0] + local v_x_adi = [v_0 * cosd(θ)] + local v_y_adi = [v_0 * sind(θ)] + local t_adi = [0.0] + + local x_iso = [0.0] + local y_iso = [0.0] + local v_x_iso = [v_0 * cosd(θ)] + local v_y_iso = [v_0 * sind(θ)] + local t_iso = [0.0] + + local x_nod = [0.0] + local y_nod = [0.0] + local v_x_nod = [v_0 * cosd(θ)] + local v_y_nod = [v_0 * sind(θ)] + local t_nod = [0.0] + + # simulate the trajectory + adiabatic!(x_adi, y_adi, v_y_adi, v_x_adi, t_adi) + isothermic!(x_iso, y_iso, v_y_iso, v_x_iso, t_iso) + nodensity!(x_nod, y_nod, v_y_nod, v_x_nod, t_nod) + + # interpolate the last point that's underground + interpolate!(x_adi, y_adi) + interpolate!(x_iso, y_iso) + interpolate!(x_nod, y_nod) + + # store the range and angle + ranges_adi[i] = x_adi[end] + angles_adi[i] = θ + ranges_iso[i] = x_iso[end] + angles_iso[i] = θ + ranges_nod[i] = x_nod[end] + angles_nod[i] = θ +end + +max_i_adi = argmax(ranges_adi) +max_i_iso = argmax(ranges_iso) +max_i_nod = argmax(ranges_nod) + +# print the max range and angle +println("Max range (adiabatic): ", ranges_adi[max_i_adi], " at angle ", angles_adi[max_i_adi], " (T=", T_0, "K)") +println("Max range (isothermic): ", ranges_iso[max_i_iso], " at angle ", angles_iso[max_i_iso], " (T=", T_0, "K)") +println("Max range (no density model): ", ranges_nod[max_i_nod], " at angle ", angles_nod[max_i_nod], " (T=", T_0, "K)") + +# plot the range over angles +plot_title = "Range over angle (T=$T_0 K)" +plot(angles_adi, ranges_adi, xlabel="angle (degrees)", ylabel="range (m)", title=plot_title, label="adiabatic", lw=2, color=:blue) +plot!(angles_iso, ranges_iso, label="isothermic", lw=2, color=:red) +plot!(angles_nod, ranges_nod, label="no density model", lw=2, color=:green) + diff --git a/hw2/2-7a.svg b/hw2/2-7a.svg new file mode 100644 index 0000000..7217e7b --- /dev/null +++ b/hw2/2-7a.svg @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw2/2-7b.png b/hw2/2-7b.png new file mode 100644 index 0000000..c864577 Binary files /dev/null and b/hw2/2-7b.png differ diff --git a/hw2/2-7c.svg b/hw2/2-7c.svg new file mode 100644 index 0000000..12c218e --- /dev/null +++ b/hw2/2-7c.svg @@ -0,0 +1,47 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw2/2-7d.png b/hw2/2-7d.png new file mode 100644 index 0000000..fcf9f17 Binary files /dev/null and b/hw2/2-7d.png differ diff --git a/hw2/2-7e.svg b/hw2/2-7e.svg new file mode 100644 index 0000000..6d6ffcf --- /dev/null +++ b/hw2/2-7e.svg @@ -0,0 +1,49 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw2/2-7f.png b/hw2/2-7f.png new file mode 100644 index 0000000..39139fc Binary files /dev/null and b/hw2/2-7f.png differ diff --git a/hw2/2-9.jl b/hw2/2-9.jl new file mode 100644 index 0000000..1190c4a --- /dev/null +++ b/hw2/2-9.jl @@ -0,0 +1,136 @@ +using Plots # for plotting trajectory + +# simulation parameters +Δt = 0.01 # time step +y_min = 0.0 +θ_to_use = [45, 35] # in degrees +v_0 = 700.0 # in m/s + +# constants +B_ref_over_m = 4.0 * 10^(-5) # in m-1, at 300K +T_ref = 300.0 # in kelvin +T_0 = 300.0 # in kelvin +g = 9.8 # in m/s^2 + +# isothermic parameters +y_o = 10 * 10^4 # k_BT/mg in meter + +# adiabatic parameters +α = 2.5 # for air +a = 6.5 * 10^(-3) # in kelvin/meter + +function adiabatic!( + x::Vector{Float64}, y::Vector{Float64}, + v_y::Vector{Float64}, v_x::Vector{Float64}, + t::Vector{Float64}) + while y[end] >= y_min + # decompose previous positions and velocities + x_i = x[end] + y_i = y[end] + v_x_i = v_x[end] + v_y_i = v_y[end] + + # calculate new positions + x_new = x_i + v_x_i * Δt + y_new = y_i + v_y_i * Δt + + # calculate drag force + v_i = sqrt(v_x_i^2 + v_y_i^2) + F_drag = - B_ref_over_m * (T_0 / T_ref)^α * # temperature variation + (1 - ((a * y_i) / T_0))^α # density/altitude variation + F_drag_x = F_drag * v_x_i * v_i + F_drag_y = F_drag * v_y_i * v_i + + # calculate new velocities + v_x_new = v_x_i + F_drag_x * Δt + v_y_new = v_y_i + F_drag_y * Δt - g * Δt + + # store new positions and velocities + push!(x, x_new) + push!(y, y_new) + push!(v_x, v_x_new) + push!(v_y, v_y_new) + push!(t, t[end] + Δt) + end +end + +function nodensity!( + x::Vector{Float64}, y::Vector{Float64}, + v_y::Vector{Float64}, v_x::Vector{Float64}, + t::Vector{Float64}) + while y[end] >= y_min + # decompose previous positions and velocities + x_i = x[end] + y_i = y[end] + v_x_i = v_x[end] + v_y_i = v_y[end] + + # calculate new positions + x_new = x_i + v_x_i * Δt + y_new = y_i + v_y_i * Δt + + # calculate drag force + v_i = sqrt(v_x_i^2 + v_y_i^2) + F_drag = - B_ref_over_m # coefficient of drag alone + F_drag_x = F_drag * v_x_i * v_i + F_drag_y = F_drag * v_y_i * v_i + + # calculate new velocities + v_x_new = v_x_i + F_drag_x * Δt + v_y_new = v_y_i + F_drag_y * Δt - g * Δt + + # store new positions and velocities + push!(x, x_new) + push!(y, y_new) + push!(v_x, v_x_new) + push!(v_y, v_y_new) + push!(t, t[end] + Δt) + end +end + +# interpolate the last point that's underground +function interpolate!(x::Vector{Float64}, y::Vector{Float64}) + if y[end] == 0 + return # no nothing if y is perfectly on 0 + end + + # calculate x_l, the interpolated x value at y=0 + r = -y[end-1] / y[end] + x_l = (x[end-1] + r * x[end]) / (1 + r) + + # set final values in the array to interpolated point on ground (y=0 ) + x[end] = x_l + y[end] = 0.0 +end + +# setup empty plot to add to +p = plot(xlabel="x (m)", ylabel="y (m)", title="Cannon Shell Trajectory", xlim=(0, 30000), xticks=0:5000:30000, legend=:topright, lw=2) +for θ in θ_to_use + # arrays to store the trajectory + x = [0.0] + y = [0.0] + v_x = [v_0 * cosd(θ)] + v_y = [v_0 * sind(θ)] + t = [0.0] + + # run the simulation + adiabatic!(x, y, v_y, v_x, t) + interpolate!(x, y) + plot_label = "adiabatic, θ = $θ" + plot!(x, y, label=plot_label, lw=2) + + # reset arrays + x = [0.0] + y = [0.0] + v_x = [v_0 * cosd(θ)] + v_y = [v_0 * sind(θ)] + t = [0.0] + + nodensity!(x, y, v_y, v_x, t) + interpolate!(x, y) + plot_label = "nodensity, θ = $θ" + plot!(x, y, label=plot_label, lw=2, linestyle=:dash) +end + +# display the plot +display(p) diff --git a/hw2/2-9a.svg b/hw2/2-9a.svg new file mode 100644 index 0000000..7593221 --- /dev/null +++ b/hw2/2-9a.svg @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/hw2/hw2-writeup.pdf b/hw2/hw2-writeup.pdf new file mode 100644 index 0000000..7f8bf5c Binary files /dev/null and b/hw2/hw2-writeup.pdf differ -- cgit v1.2.3-70-g09d2