finish hw4

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