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

12 files changed, 370 insertions, 95 deletions

diff --git a/.DS_Store b/.DS_Store
index f47adcd..ada1a6c 100644
--- a/.DS_Store
+++ b/.DS_Store
diff --git a/hw3/DrivenPendulum.jl b/hw3/3-12.jl
index 53a1af9..c83e4c3 100644
--- a/hw3/DrivenPendulum.jl
+++ b/hw3/3-12.jl
@@ -1,3 +1,6 @@
+# FOR PROBLEM 3.12
+# author: sotech117
+
 #!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
 
 # Simulate driven pendulum to find chaotic regime
diff --git a/hw3/3-13.jl b/hw3/3-13.jl
index aacfa44..dc8dec6 100644
--- a/hw3/3-13.jl
+++ b/hw3/3-13.jl
@@ -1,3 +1,6 @@
+# FOR PROBLEM 3.13
+# author: sotech117
+
 #!/Applications/Julia-1.8.app/Contents/Resources/julia/bin/julia
 
 # Simulate driven pendulum to find chaotic regime
diff --git a/hw3/Hamiltonian.jl b/hw3/Hamiltonian.jl
index bc8f31a..7cefe55 100644
--- a/hw3/Hamiltonian.jl
+++ b/hw3/Hamiltonian.jl
@@ -1,11 +1,12 @@
+# FOR PHYS2600 HW3
+# author: sotech117
+
 using DifferentialEquations
 using Plots
-using Polyhedra # for finding the volume of the phase-space ensemble
-import CDDLib
-lib = CDDLib.Library()
 
 g = 9.8
-sim_time = 100.0
+t_span = (0.0, 1000.0)
+p = g
 
 function tendency!(du, u, p, t)
     q_1, q_2, p_1, p_2 = u
@@ -17,10 +18,10 @@ function tendency!(du, u, p, t)
     dq_2 = (-p_1*cos(q_1 - q_2) + 2*p_2) / denominator
 
     k_1 = p_1*p_2*sin(q_1 - q_2) / denominator
-    k_2 = sin(2*(q_1 - q_2)) * (p_1*p_1 + 2*p_2*p_2 - 2*p_1*p_2*cos(q_1 - q_2)) / (2*denominator*denominator)
+    k_2 = (p_1*p_1 + 2*p_2*p_2 - 2*p_1*p_2*cos(q_1 - q_2)) / (2*denominator*denominator)
 
-    dp_1 = -2 * sin(q_1) - k_1 + k_2
-    dp_2 = - sin(q_2) + k_1 - k_2 
+    dp_1 = -2 * sin(q_1) - k_1 + k_2 * sin(2*(q_1 - q_2))
+    dp_2 = - sin(q_2) + k_1 - k_2 * sin(2*(q_1 - q_2))
     
     du[1] = dq_1
     du[2] = dq_2
@@ -38,21 +39,16 @@ function energy(u, p)
 end
 
 function build_ensemble(q_start, q_end, p_start, p_end, n)
-    q_1s = []
-    q_2s = []
-    p_1s = []
-    p_2s = []
+    q_1s = zeros(n)
+    q_2s = zeros(n)
+    p_1s = zeros(n)
+    p_2s = zeros(n)
     for i in 1:n
-        q_1 = q_start + (q_end - q_start) * rand()
-        q_2 = q_start + (q_end - q_start) * rand()
-        p_1 = p_start + (p_end - p_start) * rand()
-        p_2 = p_start + (p_end - p_start) * rand()
-        push!(q_1s, q_1)
-        push!(q_2s, q_2)
-        push!(p_1s, p_1)
-        push!(p_2s, p_2)
+        q_1s[i] = q_start + (q_end - q_start) * rand()
+        q_2s[i] = q_start + (q_end - q_start) * rand()
+        p_1s[i] = p_start + (p_end - p_start) * rand()
+        p_2s[i] = p_start + (p_end - p_start) * rand()
     end
-
     return [q_1s, q_2s, p_1s, p_2s]
 end
 
@@ -79,97 +75,110 @@ function run_simulation(u_0, tspan, p)
     return sol
 end
 
-function return_convex_hull(ensemble, library)
-    # plot the first phase space
-    vertices = zeros(length(ensemble[1]), 4)
-    for i in 1:length(ensemble[1])
-        vertices[i, 1] = ensemble[1][i]
-        vertices[i, 2] = ensemble[2][i]
-        vertices[i, 3] = ensemble[3][i]
-        vertices[i, 4] = ensemble[4][i]
-    end
-    v_1 = vrep(vertices)
-    p_1 = polyhedron(v_1, library)
-    vol_1 = Polyhedra.volume(p_1)
-
-    return p_1, vol_1
-end
-
 # returns the new ensemble after t_span time
 function simulate_ensembles(ensemble, tspan, p)
-    q_1s = []
-    q_2s = []
-    p_1s = []
-    p_2s = []
+    q_1s = zeros(length(ensemble[1]))
+    q_2s = zeros(length(ensemble[2]))
+    p_1s = zeros(length(ensemble[3]))
+    p_2s = zeros(length(ensemble[4]))
+    println("length ensemble = ", length(ensemble[1]))
     for i in 1:length(ensemble[1])
         u_0 = [ensemble[1][i], ensemble[2][i], ensemble[3][i], ensemble[4][i]]
-        println("init cond" ,u_0)
         sol = run_simulation(u_0, tspan, p)
-        q_1s = sol[1, :]
-        q_2s = sol[2, :]
-        p_1s = sol[3, :]
-        p_2s = sol[4, :]
-
-        push!(q_1s, q_1s[end])
-        push!(q_2s, q_2s[end])
-        push!(p_1s, p_1s[end])
-        push!(p_2s, p_2s[end])
-    end
 
+        q_1s[i] = sol[1, end]
+        q_2s[i] = sol[2, end]
+        p_1s[i] = sol[3, end]
+        p_2s[i] = sol[4, end]
+    end
+    println("q1_s", length(q_1s))
     return [q_1s, q_2s, p_1s, p_2s]
 end
 
-function calculate_ensemble_volume(ensemble)
-    dΓ = 1;
+function sum_all_positions_and_momenta(ensemble)
+    q_1 = 0
+    q_2 = 0
+    p_1 = 0
+    p_2 = 0
     for i in 1:length(ensemble[1])
-        dΓ *= (ensemble[1][i]) * (ensemble[2][i]) * (ensemble[3][i]) * ensemble[4][i]
+        q_1 += ensemble[1][i]
+        q_2 += ensemble[2][i]
+        p_1 += ensemble[3][i]
+        p_2 += ensemble[4][i]
     end
-    return dΓ
+    return [q_1, q_2, p_1, p_2]
 end
 
+function estimate_phase_space_volume(ensemble, n)
+    # make 100000 random points in the 4d phase space from [-pi, pi] x [-pi, pi] x [-3, 3] x [-3, 3]
+    points = build_ensemble(-pi, pi, -3, 3, 100000)
+
+    # iterate over the points and find the number of points that hit the ensemble
+    count = 0
+    for i in 1:length(points[1])
+        for j in 1:length(ensemble[1])
+            if abs(points[1][i] - ensemble[1][j]) < n && abs(points[2][i] - ensemble[2][j]) < n && abs(points[3][i] - ensemble[3][j]) < n && abs(points[4][i] - ensemble[4][j]) < n
+                count += 1
+                break
+            end
+        end
+    end
+
+    # return the volume
+    return (count / 100000) # note units don't matter, just comparing
+end
+
+
 # PROBLEM 1 -> show that two systems with the same intial energy can have different chaotic results
-# u_1_0 = [0, 0, 6.26, 0.0]
-# u_2_0 = [0.0, 0.0, 0.0, 6.26]
+u_1_0 = [0, 0, 6.26, 0.0]
+u_2_0 = [0.0, 0.0, 0.0, 6.26]
 
-# # print the initial energies
-# println("energy 1 = ", energy(u_1_0, p))
-# println("energy 2 = ", energy(u_2_0, p))
+# print the initial energies
+println("energy 1 = ", energy(u_1_0, p))
+println("energy 2 = ", energy(u_2_0, p))
 
-# sol_1 = run_simulation(u_1_0, tspan, p)
-# sol_2 = run_simulation(u_2_0, tspan, p)
+sol_1 = run_simulation(u_1_0, t_span, p)
+sol_2 = run_simulation(u_2_0, t_span, p)
 
-# # plot the phase space
-# p1 = scatter(sol_1[1, :], sol_1[3, :], label="q_1", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space i=1 for Pendulum 1", legend=false, color=:blue, msw=0, ms=.75)
-# p2 = scatter(sol_2[1, :], sol_2[3, :], label="q_2", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space i=1 for Pendulum 2", legend=false, color=:red, msw=0, ms=.75)
+# plot the phase space
+p1 = scatter(sol_1[1, :], sol_1[3, :], label="q_1", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space i=1 for Pendulum 1", legend=false, color=:blue, msw=0, ms=.75)
+p2 = scatter(sol_2[1, :], sol_2[3, :], label="q_2", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space i=1 for Pendulum 2", legend=false, color=:red, msw=0, ms=.75)
 
-# # plot the other phase space
-# p3 = scatter(sol_1[2, :], sol_1[4, :], label="q_1", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space i=2 for Pendulum 1", legend=false, color=:blue, msw=0, ms=.75)
-# p4 = scatter(sol_2[2, :], sol_2[4, :], label="q_2", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space i=2 for Pendulum 2", legend=false, color=:red, msw=0, ms=.75)
+# plot the other phase space
+p3 = scatter(sol_1[2, :], sol_1[4, :], label="q_1", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space i=2 for Pendulum 1", legend=false, color=:blue, msw=0, ms=.75)
+p4 = scatter(sol_2[2, :], sol_2[4, :], label="q_2", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space i=2 for Pendulum 2", legend=false, color=:red, msw=0, ms=.75)
 
-# # plot the trajectories
-# p5 = scatter(sol_1[1, :], sol_1[2, :], label="q_1 (radians)", xlabel="q_1 (radians)", ylabel="q_2 (radians)", title="Relative Trajectory for Pend1", legend=false, color=:blue, msw=0, ms=.75)
-# p6 = scatter(sol_2[1, :], sol_2[2, :], label="q_2 (radians)", xlabel="q_1 (radians)", ylabel="q_2 (radians)", title="Relative Trajectory for Pend2", legend=false, color=:red, msw=0, ms=.75)
+# plot the trajectories
+p5 = scatter(sol_1[1, :], sol_1[2, :], label="q_1 (radians)", xlabel="q_1 (radians)", ylabel="q_2 (radians)", title="Relative Trajectory for Pend1", legend=false, color=:blue, msw=0, ms=.75)
+p6 = scatter(sol_2[1, :], sol_2[2, :], label="q_2 (radians)", xlabel="q_1 (radians)", ylabel="q_2 (radians)", title="Relative Trajectory for Pend2", legend=false, color=:red, msw=0, ms=.75)
 
-# # plot the trajectories over time
-# p7 = plot(sol_1.t, sol_1[1, :], label="pendulum 1", xlabel="time (s)", ylabel="q_1 (radians)", title="q_1 vs time", legend=:topright)
-# plot!(p7, sol_2.t, sol_2[1, :], label="pendulum 2", xlabel="time (s)", ylabel="q_1 (radians)", title="q_1 vs time")
+# plot the trajectories over time
+p7 = plot(sol_1.t, sol_1[1, :], label="pendulum 1", xlabel="time (s)", ylabel="q_1 (radians)", title="q_1 vs time", legend=:topright)
+plot!(p7, sol_2.t, sol_2[1, :], label="pendulum 2", xlabel="time (s)", ylabel="q_1 (radians)", title="q_1 vs time")
 
-# p8 = plot(sol_1.t, sol_1[2, :], label="pendulum 1", xlabel="time (s)", ylabel="q_2 (radians)", title="q_2 vs time", legend=:topright)
-# plot!(p8, sol_2.t, sol_2[2, :], label="pendulum 2", xlabel="time (s)", ylabel="q_2 (radians)", title="q_2 vs time")
+p8 = plot(sol_1.t, sol_1[2, :], label="pendulum 1", xlabel="time (s)", ylabel="q_2 (radians)", title="q_2 vs time", legend=:topright)
+plot!(p8, sol_2.t, sol_2[2, :], label="pendulum 2", xlabel="time (s)", ylabel="q_2 (radians)", title="q_2 vs time")
 
-# plot(p1, p2, p3, p4, p5, p6, p7, p8, layout=(4, 2), size=(1000, 1000))
+plt = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout=(4, 2), size=(1000, 1000))
 
-# savefig("hw3/double_pendulum.png")
+savefig(plt, "hw3/double_pendulum.png")
 
 # PROBLEM 2 -> show that the volume of the phase space is conserved
 
-# plot the initial ensemble in each phase space
-ensemble = build_ensemble(0.0, pi, -1.0, 2.0, 10)
-println("4d volume before = ", return_convex_hull(ensemble, lib)[2])
-
+ensemble = build_ensemble(-1, 1, -1, 1, 1000)
+# plot these points
+s1 = scatter(ensemble[1], ensemble[3], label="t=0", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space (i=1)", color=:blue, msw=.5)
+s2 = scatter(ensemble[2], ensemble[4], label="t=0", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space (i=2)", color=:red, msw=.5)
 # simulate the ensemble
-new_ensemble = simulate_ensembles(ensemble, tspan, p)
-println("4d volume after = ", return_convex_hull(new_ensemble, lib)[2])
+new_ensemble = simulate_ensembles(ensemble, t_span, p)
+# plot these points
+scatter!(s1, new_ensemble[1], new_ensemble[3], label="t=1000", xlabel="q_1 (radians)", ylabel="p_1 (momentum)", title="Phase Space (i=1)", color=:green, msw=.5)
+scatter!(s2, new_ensemble[2], new_ensemble[4], label="t=1000", xlabel="q_2 (radians)", ylabel="p_2 (momentum)", title="Phase Space (i=2)", color=:green, msw=.5)
+
+# print the volumes
+println("\n\nOriginal Volume:\t", estimate_phase_space_volume(ensemble, .1))
+println("Final Volume:\t", estimate_phase_space_volume(new_ensemble, .1))
+
+plt_2 = plot(s1, s2)
+savefig(plt_2, "hw3/ensemble.png")
 
-# println("energy before = ", energy(ensemble, p))
-# println("volume before = ", calculate_ensemble_volume(ensemble))
diff --git a/hw3/Lorenz.jl b/hw3/Lorenz.jl
index 4a7fa0e..b1dd8f0 100644
--- a/hw3/Lorenz.jl
+++ b/hw3/Lorenz.jl
@@ -1,3 +1,6 @@
+# FOR PROBLEM 3.25
+# author: sotech117
+
 using DifferentialEquations
 using Plots
 
@@ -25,8 +28,7 @@ r_maxes = []
 z_maxes = []
 r_mins = []
 z_mins = []
-let 
-    r_diverge_1 = -1
+let
     for i in 1:r_steps
         r = r_values[i]
         p = [10.0, r, 8/3] # parameters of the Lorentz 63 system
@@ -49,12 +51,6 @@ let
                 push!(r_mins, r_values[i])
                 push!(z_mins, z_values[j])
             end
-
-            # calcualte dz
-            # if r_diverge_1 < 0 && length(z_mins) > 1 && length(z_maxes) > 1 && abs(z_mins[end] - z_maxes[end]) > 20.0
-            #     r_diverge_1 = r
-            #     println("Divergence starts at r = $r")
-            # end
         end
 
         if i % 500 == 0
@@ -70,6 +66,5 @@ let
     vline!([122], label="Critical r value ≈ 122", lc=:orange, lw=1.5, ls=:dash)
     vline!([146], label="Critical r value ≈ 146", lc=:purple, lw=1.5, ls=:dash)
 
-
     savefig("hw3/test11.png")
 end
\ No newline at end of file
diff --git a/hw3/bigdata.png b/hw3/bigdata.png
new file mode 100644
index 0000000..589f43a
--- /dev/null
+++ b/hw3/bigdata.png
diff --git a/hw3/double_pendulum.png b/hw3/double_pendulum.png
index e95959a..1024c62 100644
--- a/hw3/double_pendulum.png
+++ b/hw3/double_pendulum.png
diff --git a/hw3/ensemble.png b/hw3/ensemble.png
new file mode 100644
index 0000000..dd9fb70
--- /dev/null
+++ b/hw3/ensemble.png
diff --git a/hw3/l2.jl b/hw3/l2.jl
index 64a9cc0..2e1519c 100644
--- a/hw3/l2.jl
+++ b/hw3/l2.jl
@@ -1,3 +1,6 @@
+# FOR PROBLEM 3.14
+# author: sotech117
+
 using DifferentialEquations
 using Plots
 
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)
+