finish first two problems

31 files changed, 209 insertions, 0 deletions

diff --git a/hw4/4-13.jl b/hw4/4-13.jl
new file mode 100644
index 0000000..2a5bd9f
--- /dev/null
+++ b/hw4/4-13.jl
@@ -0,0 +1,209 @@
+#!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
+
+# Simulate a solar system
+
+using Plots # for plotting trajectory
+using Measures # for adding margins to the plots (no cut-off labels)
+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
+
+m = .00001
+
+function SolarSystem()
+
+   bodies = Vector{body}()
+   push!(bodies, body("Sun", 1.0, zeros(3), zeros(3)))
+   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("Earth", m, [-2.84, 0.0, 0.0], [0.0, 0.0, 0.0])) 
+   # -2.82 = negative ejection unstable, -2.83 positive ejection unstable
+   # -2.85 = stable , with noticable deviation
+   # -2.9 = no ejeciton, stable, -2.85 = stable, with noticable deviation
+   # .90 yields normal motion, .88 yields chaotic motion (1 ties with earth), .888 yields chaotic motion (2 ties with earth)
+   #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 = 250.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
+body3 = 3 # planet
+# 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 of body 1 orbit
+xy_b1 = plot(sol[1, body1, :], sol[2, body1, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit Star 1", color="orange", label="Star 1")
+
+# Plot of body 2 orbit
+xy_b2 = plot(sol[1, body2, :], sol[2, body2, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit Star 2", color="green", label="Star 2")
+
+# Plot of body 3 orbit
+xy_planet = plot(sol[1, body3, :], sol[2, body3, :], xlabel = "x", ylabel = "y", legend = false, title = "Orbit Earth", color="blue", label="Planet")
+
+# Plot of the orbits overlapping
+xy_all = plot(sol[1, body1, :], sol[2, body1, :], xlabel = "x", ylabel = "y", title = "Orbits (earth, x_0=-10, m=$(m))", aspect_ratio=:equal, legend=:bottomright, color="orange", label="Star 1");
+plot!(xy_all, sol[1, body2, :], sol[2, body2, :], color="green", label="Star 2")
+plot!(xy_all, sol[1, body3, :], sol[2, body3, :], color="blue", label="Planet")
+
+
+plot(
+    xy_all, xy_b1, xy_b2, xy_planet,
+    aspect_ratio=:equal, 
+    layout=(1,4), size=(1200, 300), 
+    left_margin=7mm, bottom_margin=7mm
+    )
+
+savefig("hw4/4-13-m-$(m).png")
+
+
+# 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)
+
diff --git a/hw4/4-10-dθvdt.png b/hw4/images/4-10-dθvdt.png
index 9a197ac..9a197ac 100644
--- a/hw4/4-10-dθvdt.png
+++ b/hw4/images/4-10-dθvdt.png
diff --git a/hw4/4-10-slopes-unrounded.png b/hw4/images/4-10-slopes-unrounded.png
index a1e6601..a1e6601 100644
--- a/hw4/4-10-slopes-unrounded.png
+++ b/hw4/images/4-10-slopes-unrounded.png
diff --git a/hw4/4-10-slopes.png b/hw4/images/4-10-slopes.png
index 8df2664..8df2664 100644
--- a/hw4/4-10-slopes.png
+++ b/hw4/images/4-10-slopes.png
diff --git a/hw4/4-10.png b/hw4/images/4-10.png
index 8df2664..8df2664 100644
--- a/hw4/4-10.png
+++ b/hw4/images/4-10.png
diff --git a/hw4/images/4-13-2.84.png b/hw4/images/4-13-2.84.png
new file mode 100644
index 0000000..1a9dad7
--- /dev/null
+++ b/hw4/images/4-13-2.84.png
diff --git a/hw4/images/4-13-2.85.png b/hw4/images/4-13-2.85.png
new file mode 100644
index 0000000..990b0ae
--- /dev/null
+++ b/hw4/images/4-13-2.85.png
diff --git a/hw4/images/4-13-2.86.png b/hw4/images/4-13-2.86.png
new file mode 100644
index 0000000..84c1153
--- /dev/null
+++ b/hw4/images/4-13-2.86.png
diff --git a/hw4/images/4-13-282.png b/hw4/images/4-13-282.png
new file mode 100644
index 0000000..4af1eff
--- /dev/null
+++ b/hw4/images/4-13-282.png
diff --git a/hw4/images/4-13-283.png b/hw4/images/4-13-283.png
new file mode 100644
index 0000000..320c418
--- /dev/null
+++ b/hw4/images/4-13-283.png
diff --git a/hw4/images/4-13-m-0.01.png b/hw4/images/4-13-m-0.01.png
new file mode 100644
index 0000000..8c42707
--- /dev/null
+++ b/hw4/images/4-13-m-0.01.png
diff --git a/hw4/images/4-13-m-0.1.png b/hw4/images/4-13-m-0.1.png
new file mode 100644
index 0000000..eba9b3b
--- /dev/null
+++ b/hw4/images/4-13-m-0.1.png
diff --git a/hw4/images/4-13-m-0.4.png b/hw4/images/4-13-m-0.4.png
new file mode 100644
index 0000000..1ff5ddd
--- /dev/null
+++ b/hw4/images/4-13-m-0.4.png
diff --git a/hw4/images/4-13-m-0.5.png b/hw4/images/4-13-m-0.5.png
new file mode 100644
index 0000000..0f0780b
--- /dev/null
+++ b/hw4/images/4-13-m-0.5.png
diff --git a/hw4/images/4-13-m-0.6.png b/hw4/images/4-13-m-0.6.png
new file mode 100644
index 0000000..c74c384
--- /dev/null
+++ b/hw4/images/4-13-m-0.6.png
diff --git a/hw4/images/4-13-m-0.65.png b/hw4/images/4-13-m-0.65.png
new file mode 100644
index 0000000..83e040b
--- /dev/null
+++ b/hw4/images/4-13-m-0.65.png
diff --git a/hw4/images/4-13-m-0.7.png b/hw4/images/4-13-m-0.7.png
new file mode 100644
index 0000000..6c2a0d6
--- /dev/null
+++ b/hw4/images/4-13-m-0.7.png
diff --git a/hw4/images/4-13-m-0001.png b/hw4/images/4-13-m-0001.png
new file mode 100644
index 0000000..e72b159
--- /dev/null
+++ b/hw4/images/4-13-m-0001.png
diff --git a/hw4/images/4-13-m-001.png b/hw4/images/4-13-m-001.png
new file mode 100644
index 0000000..89704c7
--- /dev/null
+++ b/hw4/images/4-13-m-001.png
diff --git a/hw4/images/4-13-m-01.png b/hw4/images/4-13-m-01.png
new file mode 100644
index 0000000..915e74f
--- /dev/null
+++ b/hw4/images/4-13-m-01.png
diff --git a/hw4/images/4-13-m-1.0e-5.png b/hw4/images/4-13-m-1.0e-5.png
new file mode 100644
index 0000000..953a62c
--- /dev/null
+++ b/hw4/images/4-13-m-1.0e-5.png
diff --git a/hw4/images/4-13-m-1.png b/hw4/images/4-13-m-1.png
new file mode 100644
index 0000000..b4d09e3
--- /dev/null
+++ b/hw4/images/4-13-m-1.png
diff --git a/hw4/images/4-13-m-10.png b/hw4/images/4-13-m-10.png
new file mode 100644
index 0000000..c5731f4
--- /dev/null
+++ b/hw4/images/4-13-m-10.png
diff --git a/hw4/images/4-13-m-525.png b/hw4/images/4-13-m-525.png
new file mode 100644
index 0000000..61e0109
--- /dev/null
+++ b/hw4/images/4-13-m-525.png
diff --git a/hw4/images/4-13-m-527.png b/hw4/images/4-13-m-527.png
new file mode 100644
index 0000000..d7b7c54
--- /dev/null
+++ b/hw4/images/4-13-m-527.png
diff --git a/hw4/images/4-13-m-55.png b/hw4/images/4-13-m-55.png
new file mode 100644
index 0000000..84a5d49
--- /dev/null
+++ b/hw4/images/4-13-m-55.png
diff --git a/hw4/images/4-13-m-59.png b/hw4/images/4-13-m-59.png
new file mode 100644
index 0000000..7165eea
--- /dev/null
+++ b/hw4/images/4-13-m-59.png
diff --git a/hw4/images/4-13-m-6.png b/hw4/images/4-13-m-6.png
new file mode 100644
index 0000000..a10f582
--- /dev/null
+++ b/hw4/images/4-13-m-6.png
diff --git a/hw4/images/4-13-m-7.png b/hw4/images/4-13-m-7.png
new file mode 100644
index 0000000..113b80a
--- /dev/null
+++ b/hw4/images/4-13-m-7.png
diff --git a/hw4/images/4-13-m-8.png b/hw4/images/4-13-m-8.png
new file mode 100644
index 0000000..ae654b7
--- /dev/null
+++ b/hw4/images/4-13-m-8.png
diff --git a/hw4/images/4-13-test.png b/hw4/images/4-13-test.png
new file mode 100644
index 0000000..9c5991a
--- /dev/null
+++ b/hw4/images/4-13-test.png