From 02756d17bca6f2b3bafa3f7b9fb6e5af438e94a0 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Thu, 29 Feb 2024 01:52:17 -0500 Subject: finish hw4 --- hw4/.DS_Store | Bin 0 -> 8196 bytes hw4/4-10.jl | 52 ++++-- hw4/4-13.jl | 2 + hw4/4-17.jl | 287 +++++++++++++++++++++++++++++++++ hw4/4-19.jl | 282 ++++++++++++++++++++++++++++++++ hw4/SolarSystem3.jl | 274 +++++++++++++++++++++++++++++++ hw4/comphysHW4-final.pdf | Bin 0 -> 1292987 bytes hw4/images/.DS_Store | Bin 0 -> 6148 bytes "hw4/images/4-10-d\316\270vdt-ldl.png" | Bin 0 -> 32065 bytes hw4/images/4-10.png | Bin 195826 -> 112631 bytes hw4/images/4-17.png | Bin 0 -> 44694 bytes hw4/images/4-19-.224.png | Bin 0 -> 159051 bytes hw4/images/4-19-0.png | Bin 0 -> 120320 bytes hw4/images/4-19-final.png | Bin 0 -> 46468 bytes hw4/images/4-19-log.png | Bin 0 -> 41271 bytes hw4/images/4-19.png | Bin 0 -> 71163 bytes hw4/images/A_alphachange.png | Bin 0 -> 29380 bytes hw4/images/l_exponent.svg | 144 +++++++++++++++++ 18 files changed, 1024 insertions(+), 17 deletions(-) create mode 100644 hw4/.DS_Store create mode 100644 hw4/4-17.jl create mode 100644 hw4/4-19.jl create mode 100644 hw4/SolarSystem3.jl create mode 100644 hw4/comphysHW4-final.pdf create mode 100644 hw4/images/.DS_Store create mode 100644 "hw4/images/4-10-d\316\270vdt-ldl.png" create mode 100644 hw4/images/4-17.png create mode 100644 hw4/images/4-19-.224.png create mode 100644 hw4/images/4-19-0.png create mode 100644 hw4/images/4-19-final.png create mode 100644 hw4/images/4-19-log.png create mode 100644 hw4/images/4-19.png create mode 100644 hw4/images/A_alphachange.png create mode 100644 hw4/images/l_exponent.svg (limited to 'hw4') diff --git a/hw4/.DS_Store b/hw4/.DS_Store new file mode 100644 index 0000000..a5beba9 Binary files /dev/null and b/hw4/.DS_Store differ diff --git a/hw4/4-10.jl b/hw4/4-10.jl index fa63926..4335c17 100644 --- a/hw4/4-10.jl +++ b/hw4/4-10.jl @@ -1,5 +1,6 @@ #!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia # FOR PROBLEM 4.10 +# author: sotech117 # Simulate planet around stationary star @@ -107,6 +108,11 @@ function linear_regression(x, y) return (a, b) end +function ldlVec(p, r, k) + return cross(p, angularMomentum(vcat(r, p))) - k*(normalize(r)) +end + + function do_α_simulations(α_values, a, e, GM_S, t_final, δ) plots = [] dθvdt = [] @@ -119,6 +125,8 @@ function do_α_simulations(α_values, a, e, GM_S, t_final, δ) p0 = [0.0, vy1_0, 0.0] # initial velocity in AU / year rp0 = vcat(r0, p0) # initial condition in phase space + ldl_0 = normalize(ldlVec(p0, r0, GM_S)) + tspan = (0.0, t_final) # span of time to simulate @@ -168,21 +176,35 @@ function do_α_simulations(α_values, a, e, GM_S, t_final, δ) color=:green, right_margin=5mm ) - # Add the linear regression of the change in angles - (a, b) = linear_regression(sol.t[i_changes], all_θ[i_changes]) - a = round(a, digits=2) - b = round(b, digits=2) - plot!(dθ, sol.t[i_changes], a * sol.t[i_changes] .+ b, label = "θ ≈ $a * t + $b", left_margin=7mm) - # Push to the plots array - push!(plots, dθ) - - # Store the slope into the return array - push!(dθvdt, a) + # # Add the linear regression of the change in angles + # (a, b) = linear_regression(sol.t[i_changes], all_θ[i_changes]) + # a = round(a, digits=2) + # b = round(b, digits=2) + # plot!(dθ, sol.t[i_changes], a * sol.t[i_changes] .+ b, label = "θ ≈ $a * t + $b", left_margin=7mm) + # # Push to the plots array + # push!(plots, dθ) + + # # Store the slope into the return array + # push!(dθvdt, a) + + # find the final ldl vector + rf = sol[1:3, end] + pf = sol[4:6, end] + ldl_f = normalize(ldlVec(pf, rf, GM_S)) + t_final = sol.t[end] + println("A_t=0 = ", ldl_0) + println("A_t=$(t_final) = ", ldl_f) + println("θ between A_t=0 and A_t = ", acos(dot(ldl_0, ldl_f))) + dθ = acos(dot(ldl_0, ldl_f)) * rad_to_deg + println("Δθ/Δt for α=$(α) = ", (dθ / t_final), "\n") + + push!(dθvdt, (dθ / t_final)) end return plots, dθvdt end + # Run the simulations (plots, dθvdt) = do_α_simulations(α_values, a, e, GM_S, t_final, δ) @@ -194,19 +216,15 @@ p_orbits = plot(plots..., layout=(4, 2), size=(800, 1000)) savefig(p_orbits, "hw4/4-10.png") # Plot the change in dθ/dt over α -p_dθvdt = scatter(α_values, dθvdt, xlabel="α", ylabel="dθ/dt (degrees/year)", title="Percession Rate vs. α", lw=2, label="dθ/dt (from regression slopes)") +p_dθvdt = scatter(α_values, dθvdt, xlabel="α", ylabel="Δθ/Δt (degrees/year)", title="Percession Rate vs. α", lw=2, label="Δθ/Δt (from Laplace-Runge-Lenz vector)") # Perform a linear regression on the data (a, b) = linear_regression(α_values, dθvdt) a = round(a, digits=2) b = round(b, digits=2) plot!(p_dθvdt, α_values, a * α_values .+ b, label = "dθ/dt ≈ $a * α + $b", left_margin=7mm) -savefig(p_dθvdt, "hw4/4-10-dθvdt.png") - +savefig(p_dθvdt, "hw4/4-10-dθvdt-ldl.png") percession_rate = a * 1.1e-8 # degrees per year to arcsec per year # convert from degrees per year to arcsec per century percession_rate = percession_rate * 3600 * 100 # arcsec per century -println("Percession rate = ", percession_rate, " arcsec/century") - - - +println("Percession rate = ", percession_rate, " arcsec/century") \ No newline at end of file diff --git a/hw4/4-13.jl b/hw4/4-13.jl index 2a5bd9f..b942ef7 100644 --- a/hw4/4-13.jl +++ b/hw4/4-13.jl @@ -1,4 +1,6 @@ #!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia +# FOR PROBLEM 4.13 +# author: sotech117 # Simulate a solar system diff --git a/hw4/4-17.jl b/hw4/4-17.jl new file mode 100644 index 0000000..eb2b3b3 --- /dev/null +++ b/hw4/4-17.jl @@ -0,0 +1,287 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia +# FOR PROBLEM 4.17 +# author: sotech117 + +# Simulate a solar system + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs +using LinearAlgebra # for dot and cross products + +println("Number of threads = ", Threads.nthreads()) + +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 + v::Vector{Float64} # velocity vector +end + +function star(name="Sun", m = 1.0, r = zeros(3), v = zeros(3)) + return body(name, m, r, v) +end + +function planet(name="Earth", m=3.0e-6, a=1.0, ϵ=0.017, i=0.0) + perihelion = (1.0 - ϵ) * a + aphelion = (1.0 + ϵ) * a + speed = sqrt(G * (1.0 + ϵ)^2 / (a * (1.0 - ϵ^2))) + phi = 0.0 + r = [perihelion * cos(i*pi/180.0) * cos(phi), perihelion * cos(i*pi/180.0) * sin(phi), perihelion * sin(i*pi/180.0)] + v = [-speed * sin(phi), speed * cos(phi), 0.0] + return body(name, m, r, v) +end + +function momentum(b::body) + return b.m * b.v +end + +function angularMomentum(b::body) + return b.m * cross(b.r, b.v) +end + +function kineticEnergy(b::body) + v = b.v + m = b.m + return 0.5 * m * dot(v, v) +end + +function potentialEnergy(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * body1.m * body2.m / sqrt(rSquared) +end + +function rv(b::body) + return [b.r; b.v] +end + +function force(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * r * rSquared^(-1.5) + #return -G * r / (rSquared * sqrt(rSquared)) +end + +mutable struct SolarSystem + bodies::Vector{body} + numberOfBodies::Int64 + phaseSpace::Matrix{Float64} # 6N dimensional phase space +end + +function SolarSystem() + + bodies = Vector{body}() + push!(bodies, star()) # default = Sun + # push!(bodies, planet("Venus", 2.44e-6, 0.72, 0.0068, 3.39)) + # push!(bodies, planet("Earth", 3.0e-6, 1.0, 0.017, 0.0)) + push!(bodies, planet("Jupiter", 0.00095, 5.2, 0.049, 1.3)) + # push!(bodies, planet("Saturn", 0.000285, 9.58, 0.055, 2.48)) + push!(bodies, body("Astroid 1", 1e-13, [3.0, 0.0, 0.0], [0.0, 3.628, 0.0])) + push!(bodies, body("Astroid 2 ", 1e-13, [3.276, 0.0, 0.0], [0.0, 3.471, 0.0])) + # push!(bodies, body("Astroid 3 (no gap)", 1e-13,[3.7, 0.0, 0.0], [0.0, 3.267, 0.0])) + numberOfBodies = size(bodies)[1] + + phaseSpace = zeros(6, 0) + for b in bodies + phaseSpace = [phaseSpace rv(b)] + end + + return SolarSystem(bodies, numberOfBodies, phaseSpace) + +end + +function TotalMass(s::SolarSystem) + M = 0.0 + for b in s.bodies + M += b.m + end + return M +end + +function TotalLinearMomentum(s::SolarSystem) + P = zeros(3) + for b in s.bodies + P += momentum(b) + end + return P +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 CenterOfMassFrame!(s::SolarSystem) + + M = TotalMass(s) + P = TotalLinearMomentum(s) + V = P / M + + s.phaseSpace = zeros(6, 0) + for body in s.bodies + body.v -= V #boost to COM frame + s.phaseSpace = [s.phaseSpace rv(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.v = 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.v #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] += b2.m * f + dps[4:6, j] -= b1.m * f # Newton's 3rd law + end + end + + return nothing +end + +function parallel_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.v = ps[4:6, i] + i += 1 + end + + # find velocities of bodies and forces on them. O(N^2) computational cost + N = s.numberOfBodies + Threads.@threads for i in 1:N + b1 = s.bodies[i] + dps[1:3, i] = b1.v #dr/dt + dps[4:6, i] = zeros(3) + for j in 1:N + if (j != i) + b2 = s.bodies[j] + f = force(b1, b2) + dps[4:6, i] += b2.m * f + end + end + end + + return nothing +end + +s = SolarSystem() +CenterOfMassFrame!(s) +println(typeof(s)) +println("Initial total energy = ", TotalEnergy(s)) +println("Initial total linear momentum = ", TotalLinearMomentum(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 = 200.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, maxiters=1e8, reltol=1e-10, abstol=1e-10) # solve using Tsit5 algorithm to specified accuracy + +println("\n\t Results") +println("Final time = ", sol.t[end]) +println("Final total energy = ", TotalEnergy(s)) +println("Final total linear momentum = ", TotalLinearMomentum(s)) +println("Final total angular momentum = ", TotalAngularMomentum(s)) + +function eccentricity(body::Int64, solution) + n = size(solution[1, 1, :])[1] + samples = 1000 + interval = floor(Int, n/samples) + ϵ = zeros(samples-1) + time = zeros(samples-1) + for i in 1:samples-1 + r2 = sol[1, body, i*interval:(i+1)*interval].^2 + sol[2, body, i*interval:(i+1)*interval].^2 + sol[3, body, i*interval:(i+1)*interval].^2 + aphelion = sqrt(maximum(r2)) + perihelion = sqrt(minimum(r2)) + ϵ[i] = (aphelion - perihelion) / (aphelion + perihelion) + time[i] = solution.t[i*interval] + end + return (time, ϵ) +end + +function obliquity(body::Int64, solution) + n = size(solution[1, 1, :])[1] + samples = 1000 + interval = floor(Int, n/samples) + ob = zeros(samples) + time = zeros(samples) + for i in 1:samples + r = solution[1:3, body, i*interval] + v = solution[4:6, body, i*interval] + ell = cross(r, v) + norm = sqrt(dot(ell, ell)) + ob[i] = (180.0 / pi) * acos(ell[3]/norm) + time[i] = solution.t[i*interval] + end + return (time, ob) +end + + +body1 = 2 # Jupiter +body2 = 3 # Astroid 1 +body3 = 4 # Astroid 2 +# body4 = 5 # Astroid 3 + +# Plot of orbit +xy = scatter([(sol[1, body1, :], sol[2, body1, :]), (sol[1, body2, :], sol[2, body2, :]), + (sol[1, body3, :], sol[2, body3, :]), + # (sol[1, body4, :], sol[2, body4, :]) + ], + xlabel = "x (AU)", ylabel = "y (AU)", title = "Orbit", + label = ["Jupiter" "Asteroid not in gap (3.0AU)" "Asteroid 2/1 gap (3.27AU)" "Asteroid not in gap (3.7AU)"], + aspect_ratio=:equal, markersize=.5, markerstrokewidth = 0, legend = :bottomright, + colors=[:red :blue :green :black] + ) + +# Plot of obliquity +tilt = obliquity(body2, sol) +obliquityPlot = plot(tilt, xlabel = "t", ylabel = "tilt", legend = false, title = "obliquity") + +# Plot of eccentricity +eccentric = eccentricity(body2, sol) +eccentricityPlot = plot(eccentric, xlabel = "t", ylabel = "ϵ", legend = false, title = "eccentricity") + +#plot(obliquityPlot, eccentricityPlot) +plot(xy, aspect_ratio=:equal) + +savefig("4-17.png") diff --git a/hw4/4-19.jl b/hw4/4-19.jl new file mode 100644 index 0000000..3208b07 --- /dev/null +++ b/hw4/4-19.jl @@ -0,0 +1,282 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia +# FOR PROBLEM 4.19 +# author: sotech117 + +# Simulate planet around stationary star + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs +using LinearAlgebra +#using Unitful +#using UnitfulAstro # astrophysical units + +G = 4.0*pi^2 # time scale = year and length scale = AU +δ = 0.0 # deviation from inverse-square law + +param = (G, δ) # parameters of force law + +function tendency!(drp, rp, param, t) + + # 6d phase space + r = rp[1:3] + p = rp[4:6] + θ = rp[7] + + + ω = rp[8] + + (G, δ) = param + + r2 = dot(r, r) + + a = -G * r * r2^(-1.5-δ) # acceleration with m = 1 + + x_c = r[1] + y_c = r[2] + + τ = -3.0 * G * r2^(-2.5 - δ) * (x_c * sin(θ) - y_c * cos(θ)) * (x_c * cos(θ) + y_c * sin(θ)) + + drp[1:3] = p[1:3] + drp[4:6] = a[1:3] + drp[7] = ω + drp[8] = τ +end + +function energy(rp, param) + + r = rp[1:3] + p = rp[4:6] + + (G, δ) = param + + r2 = dot(r, r) + + pe = -G * r2^(-0.5 - δ)/(1.0 + 2.0 * δ) + ke = 0.5 * dot(p, p) + + return pe + ke + +end + +function angularMomentum(rp) + + r = rp[1:3] + p = rp[4:6] + + return cross(r, p) + +end + +# take a list and reduce theta to the interval [-π, π] +function clean_θ(θ::Vector{Float64}) + rθ = [] + for i in 1:length(θ) + tmp = θ[i] % (2 * π) + if tmp > π + tmp = tmp - 2 * π + elseif tmp < -π + tmp = tmp + 2 * π + end + push!(rθ, tmp) + end + return rθ +end + +function simulate(r0x, p0y, θ01) + r0 = [r0x, 0.0, 0.0] # initial position in HU + p0 = [0.0, p0y, 0.0] # initial velocity in HU / year + θ0 = [θ01] + pθ0 = [0.0] + rp0 = vcat(r0, p0, θ0, pθ0) # initial condition in phase space + t_final = 10.0 # final time of simulation + + tspan = (0.0, t_final) # span of time to simulate + + + prob = ODEProblem(tendency!, rp0, tspan, param) # 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("Initial energy = ", energy(sol[:,1], param)) + println("Final energy = ", energy(sol[:, end], param)) + println("Initial angular momentum = ", angularMomentum(sol[:,1])) + println("Final angular momentum = ", angularMomentum(sol[:, end])) + + # Plot θ over time + sol[7, :] = clean_θ(sol[7, :]) + + return sol +end + +function calculate_vy1_0(a, e, GM_S) + return sqrt(GM_S * (1.0 - e)/(a * (1.0 + e))) +end + +function calculate_Δθ(sol1, sol2, ticks=1000) + # generate a equal set of times from the smaller of the two simulations + time_0 = min(sol1.t[1], sol2.t[1]) + time_f = min(sol1.t[end], sol2.t[end]) + # make an even step fucntion over the length of these times + sample_times = range(time_0, time_f, length=ticks) + + # interpolate the two simulations to the same time + diff = [] + for i in 1:length(sample_times)-1 + t = sample_times[i] + + # find the index that is before the time + idx1 = -1 + for j in 1:length(sol1.t) + if sol1.t[j] > t + idx1 = j - 1 + break + end + end + idx2 = -1 + for j in 1:length(sol2.t) + if sol2.t[j] > t + idx2 = j - 1 + break + end + end + + # interpolate the values between the two timestamps + θ1 = sol1[7, idx1] + (sol1[7, idx1+1] - sol1[7, idx1]) * (t - sol1.t[idx1])/(sol1.t[idx1+1] - sol1.t[idx1]) + θ2 = sol2[7, idx2] + (sol2[7, idx2+1] - sol2[7, idx2]) * (t - sol2.t[idx2])/(sol2.t[idx2+1] - sol2.t[idx2]) + + push!(diff, sqrt((θ1 - θ2)^2)) + end + + return (sample_times[1:end-1], diff) +end + + +function linear_regression(x, y) + n = length(x) + x̄ = sum(x) / n + ȳ = sum(y) / n + a = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄).^2) + b = ȳ - a * x̄ + return (a, b) +end + + +function find_l_exponent(eccentricity, GM_S) + a = 1.0 + vy1_0 = calculate_vy1_0(a, eccentricity, GM_S) + sol1 = simulate(1.0, vy1_0, 0.0) + sol2 = simulate(1.0, vy1_0, 0.01) # small change in θ_0 + (ts, diff_θ) = calculate_Δθ(sol1, sol2) + + plt = plot(ts, diff_θ, yscale=:ln, xlabel="time (s)", ylabel="Δθ (radians)", title="Δθ vs time", legend=:bottomright, lw=2, label="\$ Δ θ ≡ √{(θ_1 - θ_2)^2} \$") + + # linear regression of the log of the difference + (a, b) = linear_regression(ts, log.(diff_θ)) + a_rounded = round(a, digits=4) + b = round(b, digits=4) + plt = plot!(ts, exp.(a*ts .+ b), label="exp($a_rounded*t + $b)", lw=1.2, color=:red, linestyle=:dash) + + return a, plt +end + +function process_data!(les, es, plt) + # find the first index where the l exponent is greater than .2 + idx = -1 + for i in 1:length(les) + if les[i] > .2 + idx = i + break + end + end + + println("e=$(es[idx]) & l=$(les[idx])") + + # add a vertical line at this value of e + e_rounded = round(es[idx], digits=2) + plot!(plt, [es[idx], es[idx]], [-0, .75], label="e : L > 0 ≈ $(e_rounded)", color=:red, lw=.75, linestyle=:dash) + + # perform a log regression only on values where the l exponent is greater than .2 + # (a, b) = linear_regression(es[idx:end], log.(les[idx:end])) + # a_rounded = round(a, digits=2) + # b = round(b, digits=2) + # plot!(plt, es[idx:end], exp.(a*es[idx:end] .+ b), label="L ≈ exp($a_rounded*e + $b)", lw=2, color=:orange, linestyle=:dash) + + return +end + +# simualte two elipical +# eccentricity = .224 +# (a, plt) = find_l_exponent(eccentricity, G) + +# generate a span of eccentricities + +# eccentricities = range(0.0, .985, length=750) # todo: add more len for final render +# l_exponents = [] +# for e in eccentricities +# (a, _) = find_l_exponent(e, G) +# push!(l_exponents, a) +# end + +# # do a normal plot +# # p1 = plot( +# # eccentricities, l_exponents, xlabel="eccentricity (e)", ylabel="Lyapunov Exponents (L)", title="Lyapunov Exponent vs Eccentricity", +# # lw=.8, label="Lyapunov Exponents (connected)") +# p1 = scatter(eccentricities, l_exponents, +# title="Lyapunov Exponent vs Eccentricity", xlabel="Eccentricity (e)", ylabel="Lyapunov Exponent (L)", +# color=:blue, markerstrokealpha=0.0, markersize=2, msw=0.0, +# label="Lyapunov Exponents") + +# # process the data +# process_data!(l_exponents, eccentricities, p1) + +# savefig(p1, "hw4/4-19-final.png") + + +eccentricity = .224 +vy_elipitical = calculate_vy1_0(1.0, eccentricity, G) +println("vy1_0 = ", vy_elipitical) +sol1 = simulate(1.0, vy_elipitical, 0.0) +sol2 = simulate(1.0, vy_elipitical, 0.01) # small change in θ_0 +(ts, diff_θ) = calculate_Δθ(sol1, sol2) +# plot the θ over time +plt_1 = plot(sol1.t, sol1[7, :], xlabel="time (Hyperion years)", ylabel="θ (radians)", title="θ vs time \$(θ_0 = 0, e = .224\$)", legend=false, lw=2, label="θ_0 = 0.0") +# plot the ω over time +plt_2 = plot(sol2.t, sol2[7, :], xlabel="time (Hyperion years)", ylabel="θ (radians)", title="θ vs time \$(θ_0 = .01, e = .224)\$", legend=false, lw=2, label="ω_0 = 0.0") +# plot the difference in θ over time +plot_3 = plot(ts, diff_θ, yscale=:ln, xlabel="time (Hyperion years) ", ylabel="Δθ (radians)", title="Δθ vs time", legend=:bottomright, lw=2, label="\$ Δ θ ≡ √{(θ_1 - θ_2)^2} \$") + +# apply a linear regression to the log of the difference +(a, b) = linear_regression(ts, log.(diff_θ)) +a_rounded = round(a, digits=2) +b = round(b, digits=2) +plot_3 = plot!(ts, exp.(a*ts .+ b), label="exp($a_rounded*t + $b)", lw=1.2, color=:red, linestyle=:dash) + +plt_combined = plot(plt_1, plt_2, plot_3, layout=(3, 1), size=(800, 1000)) +savefig(plt_combined, "hw4/4-19-.224.png") + +# run another simuation, but only a circle +e = 0 +vy_circular = calculate_vy1_0(1.0, e, G) +println("vy1_0 = ", vy_circular) +sol1 = simulate(1.0, vy_circular, 0.0) +sol2 = simulate(1.0, vy_circular, 0.01) # small change in θ_0 +(ts, diff_θ) = calculate_Δθ(sol1, sol2) +# plot the θ over time +plt_1 = plot(sol1.t, sol1[7, :], xlabel="time (Hyperion years)", ylabel="θ (radians)", title="θ vs time (\$θ_0 = 0, e = 0\$)", lw=2, legend=false) +# plot the ω over time +plt_2 = plot(sol2.t, sol2[7, :], xlabel="time (Hyperion years)", ylabel="θ (radians)", title="θ vs time (\$θ_0 = .01, e = 0)\$", lw=2, legend = false) +# plot the difference in θ over time +plt = plot(ts, diff_θ, yscale=:ln, xlabel="time (Hyperion years)", ylabel="Δθ (radians)", title="Δθ vs time", legend=:bottomright, lw=2, label="\$ Δ θ ≡ √{(θ_1 - θ_2)^2} \$") + +# apply a linear regression to the log of the difference +(a, b) = linear_regression(ts, log.(diff_θ)) +a_rounded = round(a, digits=2) +b = round(b, digits=2) +plt = plot!(ts, exp.(a*ts .+ b), label="exp($a_rounded*t + $b)", lw=1.2, color=:red, linestyle=:dash) + +plt_combined = plot(plt_1, plt_2, plt, layout=(3, 1), size=(800, 1000)) + +savefig(plt_combined, "hw4/4-19-0.png") + diff --git a/hw4/SolarSystem3.jl b/hw4/SolarSystem3.jl new file mode 100644 index 0000000..d423669 --- /dev/null +++ b/hw4/SolarSystem3.jl @@ -0,0 +1,274 @@ +#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia + +# Simulate a solar system + +using Plots # for plotting trajectory +using DifferentialEquations # for solving ODEs +using LinearAlgebra # for dot and cross products + +println("Number of threads = ", Threads.nthreads()) + +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 + v::Vector{Float64} # velocity vector +end + +function star(name="Sun", m = 1.0, r = zeros(3), v = zeros(3)) + return body(name, m, r, v) +end + +function planet(name="Earth", m=3.0e-6, a=1.0, ϵ=0.017, i=0.0) + perihelion = (1.0 - ϵ) * a + aphelion = (1.0 + ϵ) * a + speed = sqrt(G * (1.0 + ϵ)^2 / (a * (1.0 - ϵ^2))) + phi = 0.0 + r = [perihelion * cos(i*pi/180.0) * cos(phi), perihelion * cos(i*pi/180.0) * sin(phi), perihelion * sin(i*pi/180.0)] + v = [-speed * sin(phi), speed * cos(phi), 0.0] + return body(name, m, r, v) +end + +function momentum(b::body) + return b.m * b.v +end + +function angularMomentum(b::body) + return b.m * cross(b.r, b.v) +end + +function kineticEnergy(b::body) + v = b.v + m = b.m + return 0.5 * m * dot(v, v) +end + +function potentialEnergy(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * body1.m * body2.m / sqrt(rSquared) +end + +function rv(b::body) + return [b.r; b.v] +end + +function force(body1::body, body2::body) + r = body1.r - body2.r + rSquared = dot(r, r) + return -G * r * rSquared^(-1.5) + #return -G * r / (rSquared * sqrt(rSquared)) +end + +mutable struct SolarSystem + bodies::Vector{body} + numberOfBodies::Int64 + phaseSpace::Matrix{Float64} # 6N dimensional phase space +end + +function SolarSystem() + + bodies = Vector{body}() + push!(bodies, star()) # default = Sun + push!(bodies, planet("Venus", 2.44e-6, 0.72, 0.0068, 3.39)) + push!(bodies, planet("Earth", 3.0e-6, 1.0, 0.017, 0.0)) + push!(bodies, planet("Jupiter", 0.00095, 5.2, 0.049, 1.3)) + push!(bodies, planet("Saturn", 0.000285, 9.58, 0.055, 2.48)) + numberOfBodies = size(bodies)[1] + + phaseSpace = zeros(6, 0) + for b in bodies + phaseSpace = [phaseSpace rv(b)] + end + + return SolarSystem(bodies, numberOfBodies, phaseSpace) + +end + +function TotalMass(s::SolarSystem) + M = 0.0 + for b in s.bodies + M += b.m + end + return M +end + +function TotalLinearMomentum(s::SolarSystem) + P = zeros(3) + for b in s.bodies + P += momentum(b) + end + return P +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 CenterOfMassFrame!(s::SolarSystem) + + M = TotalMass(s) + P = TotalLinearMomentum(s) + V = P / M + + s.phaseSpace = zeros(6, 0) + for body in s.bodies + body.v -= V #boost to COM frame + s.phaseSpace = [s.phaseSpace rv(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.v = 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.v #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] += b2.m * f + dps[4:6, j] -= b1.m * f # Newton's 3rd law + end + end + + return nothing +end + +function parallel_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.v = ps[4:6, i] + i += 1 + end + + # find velocities of bodies and forces on them. O(N^2) computational cost + N = s.numberOfBodies + Threads.@threads for i in 1:N + b1 = s.bodies[i] + dps[1:3, i] = b1.v #dr/dt + dps[4:6, i] = zeros(3) + for j in 1:N + if (j != i) + b2 = s.bodies[j] + f = force(b1, b2) + dps[4:6, i] += b2.m * f + end + end + end + + return nothing +end + +s = SolarSystem() +CenterOfMassFrame!(s) +println(typeof(s)) +println("Initial total energy = ", TotalEnergy(s)) +println("Initial total linear momentum = ", TotalLinearMomentum(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 = 1000.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, maxiters=1e8, reltol=1e-10, abstol=1e-10) # solve using Tsit5 algorithm to specified accuracy + +println("\n\t Results") +println("Final time = ", sol.t[end]) +println("Final total energy = ", TotalEnergy(s)) +println("Final total linear momentum = ", TotalLinearMomentum(s)) +println("Final total angular momentum = ", TotalAngularMomentum(s)) + +function eccentricity(body::Int64, solution) + n = size(solution[1, 1, :])[1] + samples = 1000 + interval = floor(Int, n/samples) + ϵ = zeros(samples-1) + time = zeros(samples-1) + for i in 1:samples-1 + r2 = sol[1, body, i*interval:(i+1)*interval].^2 + sol[2, body, i*interval:(i+1)*interval].^2 + sol[3, body, i*interval:(i+1)*interval].^2 + aphelion = sqrt(maximum(r2)) + perihelion = sqrt(minimum(r2)) + ϵ[i] = (aphelion - perihelion) / (aphelion + perihelion) + time[i] = solution.t[i*interval] + end + return (time, ϵ) +end + +function obliquity(body::Int64, solution) + n = size(solution[1, 1, :])[1] + samples = 1000 + interval = floor(Int, n/samples) + ob = zeros(samples) + time = zeros(samples) + for i in 1:samples + r = solution[1:3, body, i*interval] + v = solution[4:6, body, i*interval] + ell = cross(r, v) + norm = sqrt(dot(ell, ell)) + ob[i] = (180.0 / pi) * acos(ell[3]/norm) + time[i] = solution.t[i*interval] + end + return (time, ob) +end + + +body1 = 2 +body2 = 3 # Earth +body3 = 4 +body4 = 5 + +# Plot of orbit +xy = plot([(sol[1, body1, :], sol[2, body1, :]), (sol[1, body2, :], sol[2, body2, :]), + (sol[1, body3, :], sol[2, body3, :]), (sol[1, body4, :], sol[2, body4, :])], colors = ("yellow", "green"), xlabel = "x", ylabel = "y", legend = false, title = "Orbit") + +# Plot of obliquity +tilt = obliquity(body2, sol) +obliquityPlot = plot(tilt, xlabel = "t", ylabel = "tilt", legend = false, title = "obliquity") + +# Plot of eccentricity +eccentric = eccentricity(body2, sol) +eccentricityPlot = plot(eccentric, xlabel = "t", ylabel = "ϵ", legend = false, title = "eccentricity") + +#plot(obliquityPlot, eccentricityPlot) +plot(xy, aspect_ratio=:equal) + diff --git a/hw4/comphysHW4-final.pdf b/hw4/comphysHW4-final.pdf new file mode 100644 index 0000000..b37a6ec Binary files /dev/null and b/hw4/comphysHW4-final.pdf differ diff --git a/hw4/images/.DS_Store b/hw4/images/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/hw4/images/.DS_Store differ diff --git "a/hw4/images/4-10-d\316\270vdt-ldl.png" "b/hw4/images/4-10-d\316\270vdt-ldl.png" new file mode 100644 index 0000000..3566481 Binary files /dev/null and "b/hw4/images/4-10-d\316\270vdt-ldl.png" differ diff --git a/hw4/images/4-10.png b/hw4/images/4-10.png index 8df2664..f4803ef 100644 Binary files a/hw4/images/4-10.png and b/hw4/images/4-10.png differ diff --git a/hw4/images/4-17.png b/hw4/images/4-17.png new file mode 100644 index 0000000..4daca05 Binary files /dev/null and b/hw4/images/4-17.png differ diff --git a/hw4/images/4-19-.224.png b/hw4/images/4-19-.224.png new file mode 100644 index 0000000..4904474 Binary files /dev/null and b/hw4/images/4-19-.224.png differ diff --git a/hw4/images/4-19-0.png b/hw4/images/4-19-0.png new file mode 100644 index 0000000..a868a5b Binary files /dev/null and b/hw4/images/4-19-0.png differ diff --git a/hw4/images/4-19-final.png b/hw4/images/4-19-final.png new file mode 100644 index 0000000..4fab26d Binary files /dev/null and b/hw4/images/4-19-final.png differ diff --git a/hw4/images/4-19-log.png b/hw4/images/4-19-log.png new file mode 100644 index 0000000..b06e7df Binary files /dev/null and b/hw4/images/4-19-log.png differ diff --git a/hw4/images/4-19.png b/hw4/images/4-19.png new file mode 100644 index 0000000..ce12d9c Binary files /dev/null and b/hw4/images/4-19.png differ diff --git a/hw4/images/A_alphachange.png b/hw4/images/A_alphachange.png new file mode 100644 index 0000000..10e01a4 Binary files /dev/null and b/hw4/images/A_alphachange.png differ diff --git a/hw4/images/l_exponent.svg b/hw4/images/l_exponent.svg new file mode 100644 index 0000000..ad76bd5 --- /dev/null +++ b/hw4/images/l_exponent.svg @@ -0,0 +1,144 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file -- cgit v1.2.3-70-g09d2