finish hw3 code and start hw4. add stencil code via upload

1 files changed, 83 insertions, 0 deletions

diff --git a/hw4/SolarSystem1.jl b/hw4/SolarSystem1.jl
new file mode 100644
index 0000000..e62ee4f
--- /dev/null
+++ b/hw4/SolarSystem1.jl
@@ -0,0 +1,83 @@
+#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
+
+# 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] 
+
+   (G, δ) = param
+
+   r2 = dot(r, r)
+
+   a = -G * r * r2^(-1.5-δ) # acceleration with m = 1
+
+   drp[1:3] = p[1:3]
+   drp[4:6] = a[1:3]
+   
+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
+
+r0 = [1.0, 0.0, 0.0] # initial position in AU
+p0 = [0.0, 1.0*pi, 0.0] # initial velocity in AU / year
+rp0 = vcat(r0, p0) # 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 of position vs. time
+xt = plot(sample_times, sol[1, :], xlabel = "t", ylabel = "x(t)", legend = false, title = "x vs. t")
+
+# Plot of orbit
+xy = plot(sol[1, :], sol[2, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit", aspect_ratio=:equal)
+
+plot(xt, xy)