From ec622e808d84956a75c2b9a9826479cb5737e04a Mon Sep 17 00:00:00 2001 From: sotech117 Date: Wed, 28 Feb 2024 01:17:16 -0500 Subject: finish first two problems --- "hw4/4-10-d\316\270vdt.png" | Bin 29668 -> 0 bytes hw4/4-10-slopes-unrounded.png | Bin 203611 -> 0 bytes hw4/4-10-slopes.png | Bin 195826 -> 0 bytes hw4/4-10.png | Bin 195826 -> 0 bytes hw4/4-13.jl | 209 +++++++++++++++++++++++++++++++++++ "hw4/images/4-10-d\316\270vdt.png" | Bin 0 -> 29668 bytes hw4/images/4-10-slopes-unrounded.png | Bin 0 -> 203611 bytes hw4/images/4-10-slopes.png | Bin 0 -> 195826 bytes hw4/images/4-10.png | Bin 0 -> 195826 bytes hw4/images/4-13-2.84.png | Bin 0 -> 60254 bytes hw4/images/4-13-2.85.png | Bin 0 -> 78298 bytes hw4/images/4-13-2.86.png | Bin 0 -> 65022 bytes hw4/images/4-13-282.png | Bin 0 -> 54948 bytes hw4/images/4-13-283.png | Bin 0 -> 52838 bytes hw4/images/4-13-m-0.01.png | Bin 0 -> 55617 bytes hw4/images/4-13-m-0.1.png | Bin 0 -> 69452 bytes hw4/images/4-13-m-0.4.png | Bin 0 -> 46827 bytes hw4/images/4-13-m-0.5.png | Bin 0 -> 44175 bytes hw4/images/4-13-m-0.6.png | Bin 0 -> 66389 bytes hw4/images/4-13-m-0.65.png | Bin 0 -> 50937 bytes hw4/images/4-13-m-0.7.png | Bin 0 -> 54637 bytes hw4/images/4-13-m-0001.png | Bin 0 -> 54541 bytes hw4/images/4-13-m-001.png | Bin 0 -> 57222 bytes hw4/images/4-13-m-01.png | Bin 0 -> 56169 bytes hw4/images/4-13-m-1.0e-5.png | Bin 0 -> 68171 bytes hw4/images/4-13-m-1.png | Bin 0 -> 44533 bytes hw4/images/4-13-m-10.png | Bin 0 -> 57975 bytes hw4/images/4-13-m-525.png | Bin 0 -> 45469 bytes hw4/images/4-13-m-527.png | Bin 0 -> 43497 bytes hw4/images/4-13-m-55.png | Bin 0 -> 44530 bytes hw4/images/4-13-m-59.png | Bin 0 -> 47383 bytes hw4/images/4-13-m-6.png | Bin 0 -> 47447 bytes hw4/images/4-13-m-7.png | Bin 0 -> 54522 bytes hw4/images/4-13-m-8.png | Bin 0 -> 53484 bytes hw4/images/4-13-test.png | Bin 0 -> 41774 bytes 35 files changed, 209 insertions(+) delete mode 100644 "hw4/4-10-d\316\270vdt.png" delete mode 100644 hw4/4-10-slopes-unrounded.png delete mode 100644 hw4/4-10-slopes.png delete mode 100644 hw4/4-10.png create mode 100644 hw4/4-13.jl create mode 100644 "hw4/images/4-10-d\316\270vdt.png" create mode 100644 hw4/images/4-10-slopes-unrounded.png create mode 100644 hw4/images/4-10-slopes.png create mode 100644 hw4/images/4-10.png create mode 100644 hw4/images/4-13-2.84.png create mode 100644 hw4/images/4-13-2.85.png create mode 100644 hw4/images/4-13-2.86.png create mode 100644 hw4/images/4-13-282.png create mode 100644 hw4/images/4-13-283.png create mode 100644 hw4/images/4-13-m-0.01.png create mode 100644 hw4/images/4-13-m-0.1.png create mode 100644 hw4/images/4-13-m-0.4.png create mode 100644 hw4/images/4-13-m-0.5.png create mode 100644 hw4/images/4-13-m-0.6.png create mode 100644 hw4/images/4-13-m-0.65.png create mode 100644 hw4/images/4-13-m-0.7.png create mode 100644 hw4/images/4-13-m-0001.png create mode 100644 hw4/images/4-13-m-001.png create mode 100644 hw4/images/4-13-m-01.png create mode 100644 hw4/images/4-13-m-1.0e-5.png create mode 100644 hw4/images/4-13-m-1.png create mode 100644 hw4/images/4-13-m-10.png create mode 100644 hw4/images/4-13-m-525.png create mode 100644 hw4/images/4-13-m-527.png create mode 100644 hw4/images/4-13-m-55.png create mode 100644 hw4/images/4-13-m-59.png create mode 100644 hw4/images/4-13-m-6.png create mode 100644 hw4/images/4-13-m-7.png create mode 100644 hw4/images/4-13-m-8.png create mode 100644 hw4/images/4-13-test.png diff --git "a/hw4/4-10-d\316\270vdt.png" "b/hw4/4-10-d\316\270vdt.png" deleted file mode 100644 index 9a197ac..0000000 Binary files "a/hw4/4-10-d\316\270vdt.png" and /dev/null differ diff --git a/hw4/4-10-slopes-unrounded.png b/hw4/4-10-slopes-unrounded.png deleted file mode 100644 index a1e6601..0000000 Binary files a/hw4/4-10-slopes-unrounded.png and /dev/null differ diff --git a/hw4/4-10-slopes.png b/hw4/4-10-slopes.png deleted file mode 100644 index 8df2664..0000000 Binary files a/hw4/4-10-slopes.png and /dev/null differ diff --git a/hw4/4-10.png b/hw4/4-10.png deleted file mode 100644 index 8df2664..0000000 Binary files a/hw4/4-10.png and /dev/null differ diff --git a/hw4/4-13.jl b/hw4/4-13.jl new file mode 100644 index 0000000..2a5bd9f --- /dev/null +++ b/hw4/4-13.jl @@ -0,0 +1,209 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +# Simulate a solar system + +using Plots # for plotting trajectory +using Measures # for adding margins to the plots (no cut-off labels) +using DifferentialEquations # for solving ODEs +using LinearAlgebra # for dot and cross products + +G = 4.0*pi^2 # time scale = year and length scale = AU + +mutable struct body + name::String # name of star or planet + m::Float64 # mass + r::Vector{Float64} # position vector + p::Vector{Float64} # momentum vector +end + +function angularMomentum(b::body) + r = b.r + p = b.p + return cross(r, p) +end + +function kineticEnergy(b::body) + p = b.p + m = b.m + return dot(p, p) / (2.0 * m) +end + +function rp(b::body) + return [b.r; b.p] +end + +function force(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * body1.m * body2.m * r * rSquared^(-1.5) +end + +function potentialEnergy(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * body1.m * body2.m * rSquared^(-0.5) +end + +mutable struct SolarSystem + bodies::Vector{body} + numberOfBodies::Int64 + phaseSpace::Matrix{Float64} # 6N dimensional phase space +end + +m = .00001 + +function SolarSystem() + + bodies = Vector{body}() + push!(bodies, body("Sun", 1.0, zeros(3), zeros(3))) + push!(bodies, body("Star", 1.0, [1.0, 0.0, 0.0], [0.0, sqrt(1.95) * 2.0 * pi, 0.0])) + push!(bodies, body("Earth", m, [-2.84, 0.0, 0.0], [0.0, 0.0, 0.0])) + # -2.82 = negative ejection unstable, -2.83 positive ejection unstable + # -2.85 = stable , with noticable deviation + # -2.9 = no ejeciton, stable, -2.85 = stable, with noticable deviation + # .90 yields normal motion, .88 yields chaotic motion (1 ties with earth), .888 yields chaotic motion (2 ties with earth) + #push!(bodies, body("Jupiter", 1.0, [3.0, 0.0, 0.0], [0.0, 0.25 * pi, 0.0])) + numberOfBodies = size(bodies)[1] + + phaseSpace = zeros(6, 0) + for b in bodies + phaseSpace = [phaseSpace rp(b)] + end + + return SolarSystem(bodies, numberOfBodies, phaseSpace) + +end + +function TotalAngularMomentum(s::SolarSystem) + L = zeros(3) + for b in s.bodies + L += angularMomentum(b) + end + + return L +end + +function TotalEnergy(s::SolarSystem) + ke = 0.0 + pe = 0.0 + + for body1 in s.bodies + ke += kineticEnergy(body1) + for body2 in s.bodies + if (body1 != body2) + pe += 0.5 * potentialEnergy(body1, body2) + end + end + end + + return pe + ke +end + +function ZeroOutLinearMomentum!(s::SolarSystem) + + totalLinearMomentum = zeros(3) + totalMass = 0.0 + for body in s.bodies + totalLinearMomentum += body.p + totalMass += body.m + end + + s.phaseSpace = zeros(6, 0) + for body in s.bodies + body.p -= body.m * totalLinearMomentum / totalMass + s.phaseSpace = [s.phaseSpace rp(body)] + end + + return nothing +end + +function tendency!(dps, ps, s::SolarSystem, t) + + i = 1 # update phase space of individual bodies + for b in s.bodies + b.r = ps[1:3, i] + b.p = ps[4:6, i] + i += 1 + end + + # find velocities of bodies and forces on them. O(N^2) computational cost + N = s.numberOfBodies + for i in 1:N + b1 = s.bodies[i] + dps[1:3, i] = b1.p / b1.m #dr/dt + dps[4:6, i] = zeros(3) + for j in 1:i-1 + b2 = s.bodies[j] + f = force(b1, b2) # call only once + dps[4:6, i] += f + dps[4:6, j] -= f # Newton's 3rd law + end + end + + return nothing +end + +s = SolarSystem() +ZeroOutLinearMomentum!(s) +println(typeof(s)) +println("Initial total energy = ", TotalEnergy(s)) +println("Initial total angular momentum = ", TotalAngularMomentum(s)) +println("Number of bodies = ", s.numberOfBodies) +for b in s.bodies + println("body name = ", b.name) +end + +t_final = 250.0 # final time of simulation +tspan = (0.0, t_final) # span of time to simulate +prob = ODEProblem(tendency!, s.phaseSpace, tspan, s) # specify ODE +sol = solve(prob, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy + +sample_times = sol.t +println("\n\t Results") +println("Final time = ", sample_times[end]) +println("Final total energy = ", TotalEnergy(s)) +println("Final total angular momentum = ", TotalAngularMomentum(s)) + +body1 = 1 +body2 = 2 +body3 = 3 # planet +# Plot of position vs. time +xt = plot(sample_times, sol[1, body1, :], xlabel = "t", ylabel = "x(t)", legend = false, title = "x vs. t") + +# Plot of orbit +xy = plot([(sol[1, body1, :], sol[2, body1, :]), (sol[1, body2, :], sol[2, body2, :])], colors = ("yellow", "green"), xlabel = "x", ylabel = "y", legend = false, title = "Orbit") + +# Plot of body 1 orbit +xy_b1 = plot(sol[1, body1, :], sol[2, body1, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit Star 1", color="orange", label="Star 1") + +# Plot of body 2 orbit +xy_b2 = plot(sol[1, body2, :], sol[2, body2, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit Star 2", color="green", label="Star 2") + +# Plot of body 3 orbit +xy_planet = plot(sol[1, body3, :], sol[2, body3, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit Earth", color="blue", label="Planet") + +# Plot of the orbits overlapping +xy_all = plot(sol[1, body1, :], sol[2, body1, :], xlabel = "x", ylabel = "y", title = "Orbits (earth, x_0=-10, m=$(m))", aspect_ratio=:equal, legend=:bottomright, color="orange", label="Star 1"); +plot!(xy_all, sol[1, body2, :], sol[2, body2, :], color="green", label="Star 2") +plot!(xy_all, sol[1, body3, :], sol[2, body3, :], color="blue", label="Planet") + + +plot( + xy_all, xy_b1, xy_b2, xy_planet, + aspect_ratio=:equal, + layout=(1,4), size=(1200, 300), + left_margin=7mm, bottom_margin=7mm + ) + +savefig("hw4/4-13-m-$(m).png") + + +# samples = 1000 +# interval = floor(Int,size(sol.t)[1] / samples) +# animation = @animate for i=1:samples-1 +# plot([(sol[1, body1, i*interval], sol[2, body1, i*interval]), (sol[1, body1, (i+1)*interval], sol[2, body1, (i+1)*interval])], +# aspect_ratio=:equal, colors = ("yellow", "green"), xlabel = "x", ylabel = "y", legend = false, title = "Orbit", xlims=(-1, 1), ylims = (-1, 1)) +# end + +# gif(animation, fps=15) + diff --git "a/hw4/images/4-10-d\316\270vdt.png" "b/hw4/images/4-10-d\316\270vdt.png" new file mode 100644 index 0000000..9a197ac Binary files /dev/null and "b/hw4/images/4-10-d\316\270vdt.png" differ diff --git a/hw4/images/4-10-slopes-unrounded.png b/hw4/images/4-10-slopes-unrounded.png new file mode 100644 index 0000000..a1e6601 Binary files /dev/null and b/hw4/images/4-10-slopes-unrounded.png differ diff --git a/hw4/images/4-10-slopes.png b/hw4/images/4-10-slopes.png new file mode 100644 index 0000000..8df2664 Binary files /dev/null and b/hw4/images/4-10-slopes.png differ diff --git a/hw4/images/4-10.png b/hw4/images/4-10.png new file mode 100644 index 0000000..8df2664 Binary files /dev/null and b/hw4/images/4-10.png differ diff --git a/hw4/images/4-13-2.84.png b/hw4/images/4-13-2.84.png new file mode 100644 index 0000000..1a9dad7 Binary files /dev/null and b/hw4/images/4-13-2.84.png differ diff --git a/hw4/images/4-13-2.85.png b/hw4/images/4-13-2.85.png new file mode 100644 index 0000000..990b0ae Binary files /dev/null and b/hw4/images/4-13-2.85.png differ diff --git a/hw4/images/4-13-2.86.png b/hw4/images/4-13-2.86.png new file mode 100644 index 0000000..84c1153 Binary files /dev/null and b/hw4/images/4-13-2.86.png differ diff --git a/hw4/images/4-13-282.png b/hw4/images/4-13-282.png new file mode 100644 index 0000000..4af1eff Binary files /dev/null and b/hw4/images/4-13-282.png differ diff --git a/hw4/images/4-13-283.png b/hw4/images/4-13-283.png new file mode 100644 index 0000000..320c418 Binary files /dev/null and b/hw4/images/4-13-283.png differ diff --git a/hw4/images/4-13-m-0.01.png b/hw4/images/4-13-m-0.01.png new file mode 100644 index 0000000..8c42707 Binary files /dev/null and b/hw4/images/4-13-m-0.01.png differ diff --git a/hw4/images/4-13-m-0.1.png b/hw4/images/4-13-m-0.1.png new file mode 100644 index 0000000..eba9b3b Binary files /dev/null and b/hw4/images/4-13-m-0.1.png differ diff --git a/hw4/images/4-13-m-0.4.png b/hw4/images/4-13-m-0.4.png new file mode 100644 index 0000000..1ff5ddd Binary files /dev/null and b/hw4/images/4-13-m-0.4.png differ diff --git a/hw4/images/4-13-m-0.5.png b/hw4/images/4-13-m-0.5.png new file mode 100644 index 0000000..0f0780b Binary files /dev/null and b/hw4/images/4-13-m-0.5.png differ diff --git a/hw4/images/4-13-m-0.6.png b/hw4/images/4-13-m-0.6.png new file mode 100644 index 0000000..c74c384 Binary files /dev/null and b/hw4/images/4-13-m-0.6.png differ diff --git a/hw4/images/4-13-m-0.65.png b/hw4/images/4-13-m-0.65.png new file mode 100644 index 0000000..83e040b Binary files /dev/null and b/hw4/images/4-13-m-0.65.png differ diff --git a/hw4/images/4-13-m-0.7.png b/hw4/images/4-13-m-0.7.png new file mode 100644 index 0000000..6c2a0d6 Binary files /dev/null and b/hw4/images/4-13-m-0.7.png differ diff --git a/hw4/images/4-13-m-0001.png b/hw4/images/4-13-m-0001.png new file mode 100644 index 0000000..e72b159 Binary files /dev/null and b/hw4/images/4-13-m-0001.png differ diff --git a/hw4/images/4-13-m-001.png b/hw4/images/4-13-m-001.png new file mode 100644 index 0000000..89704c7 Binary files /dev/null and b/hw4/images/4-13-m-001.png differ diff --git a/hw4/images/4-13-m-01.png b/hw4/images/4-13-m-01.png new file mode 100644 index 0000000..915e74f Binary files /dev/null and b/hw4/images/4-13-m-01.png differ diff --git a/hw4/images/4-13-m-1.0e-5.png b/hw4/images/4-13-m-1.0e-5.png new file mode 100644 index 0000000..953a62c Binary files /dev/null and b/hw4/images/4-13-m-1.0e-5.png differ diff --git a/hw4/images/4-13-m-1.png b/hw4/images/4-13-m-1.png new file mode 100644 index 0000000..b4d09e3 Binary files /dev/null and b/hw4/images/4-13-m-1.png differ diff --git a/hw4/images/4-13-m-10.png b/hw4/images/4-13-m-10.png new file mode 100644 index 0000000..c5731f4 Binary files /dev/null and b/hw4/images/4-13-m-10.png differ diff --git a/hw4/images/4-13-m-525.png b/hw4/images/4-13-m-525.png new file mode 100644 index 0000000..61e0109 Binary files /dev/null and b/hw4/images/4-13-m-525.png differ diff --git a/hw4/images/4-13-m-527.png b/hw4/images/4-13-m-527.png new file mode 100644 index 0000000..d7b7c54 Binary files /dev/null and b/hw4/images/4-13-m-527.png differ diff --git a/hw4/images/4-13-m-55.png b/hw4/images/4-13-m-55.png new file mode 100644 index 0000000..84a5d49 Binary files /dev/null and b/hw4/images/4-13-m-55.png differ diff --git a/hw4/images/4-13-m-59.png b/hw4/images/4-13-m-59.png new file mode 100644 index 0000000..7165eea Binary files /dev/null and b/hw4/images/4-13-m-59.png differ diff --git a/hw4/images/4-13-m-6.png b/hw4/images/4-13-m-6.png new file mode 100644 index 0000000..a10f582 Binary files /dev/null and b/hw4/images/4-13-m-6.png differ diff --git a/hw4/images/4-13-m-7.png b/hw4/images/4-13-m-7.png new file mode 100644 index 0000000..113b80a Binary files /dev/null and b/hw4/images/4-13-m-7.png differ diff --git a/hw4/images/4-13-m-8.png b/hw4/images/4-13-m-8.png new file mode 100644 index 0000000..ae654b7 Binary files /dev/null and b/hw4/images/4-13-m-8.png differ diff --git a/hw4/images/4-13-test.png b/hw4/images/4-13-test.png new file mode 100644 index 0000000..9c5991a Binary files /dev/null and b/hw4/images/4-13-test.png differ -- cgit v1.2.3-70-g09d2