1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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)
|
