diff options
Diffstat (limited to 'hw4/4-10.jl')
| -rw-r--r-- | hw4/4-10.jl | 52 |
1 files changed, 35 insertions, 17 deletions
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 |