diff --git a/julia/leapfrog.jl b/julia/leapfrog.jl index 3feaa4d..7d7be41 100644 --- a/julia/leapfrog.jl +++ b/julia/leapfrog.jl @@ -1,19 +1,14 @@ using StaticArrays +using LinearAlgebra: norm -const dims = 3 -const PosVector = SVector{dims, Float64} -const zerov = zero(SVector{dims}) - -@inline function gravitate!(x::Vector{PosVector}, a::Vector{PosVector}, - m::Vector) +function gravitate!(x::Vector{T}, a::Vector{T},m::AbstractVector) where T G = 6.6742367e-11 - fill!(a, zerov) - @inbounds for i in 1:length(m) - for j in 1:length(m) + fill!(a, zero(T)) + for i in eachindex(x) + for j in eachindex(x) if i == j continue end - d = x[i] - x[j] prefactor = -G/norm(d)^3 * m[j] a[i] += prefactor*d @@ -21,37 +16,35 @@ const zerov = zero(SVector{dims}) end end -function main(;dt::Real=0.08, tmax::Real=3.6525e8) - - # Arrays (initialized to zero by default) - x = [zerov, - @SVector([0.0162, 6.57192058353e-15, 5.74968548652e-16])] - - v = [zerov, - @SVector([-1.48427302304e-14, 0.0399408809121, 0.00349437429104])] - - a = [zerov, zerov] - - m = [0.08, 3.0e-6] # M_Sun - +function main( + x::Vector{T},v::Vector{T},m::AbstractVector; + dt::Real=0.08, tmax::Real=3.6525e8 +) where T + a = similar(x) hdt = dt / 2 - - @inbounds for k in 1:round(Int, tmax/dt) - for i in 1:length(m) - x[i] += hdt*v[i] - end + for k in 1:round(Int,tmax/dt) + @. x += hdt * v gravitate!(x, a, m) - for i in 1:length(m) - v[i] += a[i]*dt - x[i] += v[i]*hdt - end + @. v += a*dt + @. x += v*hdt end - println("Positions:", x) end -println("Running up to tmax=10 to compile...") -@time main(tmax=1e2) println("Running full benchmark...") -@time main() + +# Dimension of particles +T = SVector{3,Float64} + +# Initial positions and velocities +x = [zero(T),T(0.0162, 6.57192058353e-15, 5.74968548652e-16)] +v = [zero(T),T(-1.48427302304e-14, 0.0399408809121, 0.00349437429104)] + +# Masses +m = @SVector [0.08, 3.0e-6] + +# Run simulation +main(x,v,m) + + diff --git a/julia/target/optimized/leapfrog b/julia/target/optimized/leapfrog index d242e9e..d99b0b3 100755 --- a/julia/target/optimized/leapfrog +++ b/julia/target/optimized/leapfrog @@ -1,3 +1,3 @@ #!/bin/bash # OSX: /Applications/Julia-0.5.app/Contents/Resources/julia/bin/julia -julia -O3 leapfrog.jl +julia -O3 --check-bounds=no leapfrog.jl \ No newline at end of file