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

3 files changed, 262 insertions, 0 deletions

diff --git a/hw4/4-10.jl b/hw4/4-10.jl
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/hw4/4-10.jl
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)
diff --git a/hw4/SolarSystem2.jl b/hw4/SolarSystem2.jl
new file mode 100644
index 0000000..b2e02b7
--- /dev/null
+++ b/hw4/SolarSystem2.jl
@@ -0,0 +1,179 @@
+#!/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
+
+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
+
+function SolarSystem()
+
+   bodies = Vector{body}()
+   push!(bodies, body("Sun", 1.0, zeros(3), zeros(3)))
+   #push!(bodies, body("Earth", 1.0, [-1.0, 0.0, 0.0], [0.0, 1.0 * pi, 0.0]))
+   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("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 = 10.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
+# 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(xy, aspect_ratio=:equal)
+=#
+
+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)
+