GitHub

Package Guide

PhysicalParticles.jl is useful for particle based scientific simulations

Installation

From the Julia REPL, type ] to enter the Pkg REPL mode and run

pkg> add PhysicalParticles

or add from git repository

pkg> add https://github.com/JuliaAstroSim/PhysicalParticles.jl

Test the package by

pkg> test PhysicalParticles

Basic Usage

Vectors

julia> using PhysicalParticles, Unitful, UnitfulAstro
julia> a = PVector()PVector{Float64}(0.0, 0.0, 0.0)
julia> b = PVector(1.0u"m", 2.0u"m", 3.0u"m")PVector(1.0 m, 2.0 m, 3.0 m)
julia> c = PVector2D(u"m/s")PVector2D(0.0 m s^-1, 0.0 m s^-1)
julia> uconvert(u"m", PVector(1.0, 1.0, 1.0, u"km"))PVector(1000.0 m, 1000.0 m, 1000.0 m)
julia> PVector(BigFloat)PVector{BigFloat}(0.0, 0.0, 0.0)
julia> PVector2D(BigInt, u"m")PVector2D(0 m, 0 m)
julia> PVector(1.0, 1.0) * imPVector2D{ComplexF64}(0.0 + 1.0im, 0.0 + 1.0im)
julia> b * 2.0u"s"PVector(2.0 m s, 4.0 m s, 6.0 m s)
julia> b + PVector(2.0, 2.0, 2.0, u"m") / 2PVector(2.0 m, 3.0 m, 4.0 m)
julia> norm(PVector2D(3.0f0,4.0f0))5.0f0
julia> normalize(PVector(3.0, 4.0))PVector2D{Float64}(0.6, 0.8)
julia> d = PVector(3u"kpc", 4u"kpc")PVector2D(3 kpc, 4 kpc)
julia> norm(d)5.0 kpc
julia> distance(PVector2D(0.0, 0.0), PVector2D(3.0, 4.0))5.0
julia> rotate(PVector(1.0, 0.0), 0.5pi)PVector2D{Float64}(6.123233995736766e-17, 1.0)
julia> rotate(PVector(1.0, 0.0, 0.0), 0.0, 0.0, 0.5pi)PVector{Float64}(6.123233995736766e-17, 1.0, 0.0)
julia> rotate_z(PVector(1.0, 0.0, 0.0), 90.0u"°")PVector{Float64}(0.0, 1.0, 0.0)
julia> rotate(PVector(1.0, 0.0, 0.0), 0.0, 0.0, 90.0u"°", PVector(-1.0, 0.0, 0.0))PVector{Float64}(-1.0, 2.0, 0.0)
julia> rotate(PVector(0.0, 1.0, 0.0), PVector(0.0, 1.0, 1.0), pi)PVector{Float64}(0.0, -2.220446049250313e-16, 0.9999999999999998)
julia> cartesian2cylinderial(PVector(sqrt(2), sqrt(2), 1.0, u"m"))(2.0 m, 0.7853981633974484, 1.0 m)
julia> cylinderial2cartesian(2.0u"m", pi/4, 1.0u"m")PVector(1.4142135623730951 m, 1.414213562373095 m, 1.0 m)
julia> cartesian2spherical(PVector(sqrt(0.5), sqrt(0.5), 1.0, u"m"))(1.4142135623730951 m, 0.7853981633974484, 0.7853981633974483)
julia> spherical2cartesian(sqrt(2)u"m", pi/4, pi/4)PVector(0.7071067811865476 m, 0.7071067811865475 m, 1.0000000000000002 m)
julia> zero(PVector{Float64})PVector{Float64}(0.0, 0.0, 0.0)
julia> iszero(PVector(u"m"))true
julia> isnan(PVector(NaN, NaN))true
julia> PVector2D(1.0, 1.0) ≈ PVector2D(1.0 + 1.0e-8, 1.0 + 1.0e-8)true
julia> ustrip(PVector(1.0, 1.0, 1.0, u"km"))PVector{Float64}(1.0, 1.0, 1.0)
julia> ustrip(u"m", PVector(1.0, 1.0, 1.0, u"km"))PVector{Float64}(1000.0, 1000.0, 1000.0)

Particles

We provide 2D version for each type below, for example, the 2D version of Ball is Ball2D:

julia> Massless()Massless 0: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0)
julia> Massless(PVector(0.0, 0.0, 0.0), PVector(), 1)Massless 1: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0)
julia> Massless2D(uCGS)Massless 0: Pos = PVector2D(0.0 cm, 0.0 cm), Vel = PVector2D(0.0 cm s^-1, 0.0 cm s^-1)
julia> Ball()Ball 0: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0), Acc = PVector{Float64}(0.0, 0.0, 0.0), Mass = 0.0
julia> Ball(PVector(0.0u"m", 0.0u"m", 0.0u"m"), PVector(u"m/s"), PVector(u"m/s^2"), 0.0u"kg", 1)Ball 1: Pos = PVector(0.0 m, 0.0 m, 0.0 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
julia> Star()Star 0 STAR: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0), Acc = PVector{Float64}(0.0, 0.0, 0.0), Mass = 0.0, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0, OldAcc = 0.0
julia> Star2D(uAstro)Star 0 STAR: Pos = PVector2D(0.0 kpc, 0.0 kpc), Vel = PVector2D(0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1), Acc = PVector2D(0.0 kpc Gyr^-2, 0.0 kpc Gyr^-2), Mass = 0.0 M⊙, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0 kpc^2 M⊙ Gyr^-2, OldAcc = 0.0 kpc Gyr^-2
julia> a = Star(uAstro)Star 0 STAR: Pos = PVector(0.0 kpc, 0.0 kpc, 0.0 kpc), Vel = PVector(0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1), Acc = PVector(0.0 kpc Gyr^-2, 0.0 kpc Gyr^-2, 0.0 kpc Gyr^-2), Mass = 0.0 M⊙, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0 kpc^2 Gyr^-2, OldAcc = 0.0 kpc Gyr^-2
julia> b = Star(uAstro)Star 0 STAR: Pos = PVector(0.0 kpc, 0.0 kpc, 0.0 kpc), Vel = PVector(0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1), Acc = PVector(0.0 kpc Gyr^-2, 0.0 kpc Gyr^-2, 0.0 kpc Gyr^-2), Mass = 0.0 M⊙, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0 kpc^2 Gyr^-2, OldAcc = 0.0 kpc Gyr^-2
julia> distance(a,b)0.0 kpc

Mutate array of immutable particles

To simulate particle dynamics on GPU (bit types), and for better performance, particle data are stored in immutable sturcts. One of the shortcomings is that you have to allocate a new object when changing member data, however, memory operation is optimized when doing this in array:

julia> using BangBang
julia> p = Star()Star 0 STAR: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0), Acc = PVector{Float64}(0.0, 0.0, 0.0), Mass = 0.0, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0, OldAcc = 0.0
julia> p = setproperties!!(p, Mass = 321.0)Star 0 STAR: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0), Acc = PVector{Float64}(0.0, 0.0, 0.0), Mass = 321.0, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0, OldAcc = 0.0
julia> a = [Star()]1-element Vector{Star{Float64, Float64, Float64, Float64, Float64, Int64}}:
 Star 0 STAR: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0), Acc = PVector{Float64}(0.0, 0.0, 0.0), Mass = 0.0, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0, OldAcc = 0.0
julia> a[1] = setproperties!!(a[1], Mass = 321.0) # In general, `Mass` is mutated right at its original memory address.Star 0 STAR: Pos = PVector{Float64}(0.0, 0.0, 0.0), Vel = PVector{Float64}(0.0, 0.0, 0.0), Acc = PVector{Float64}(0.0, 0.0, 0.0), Mass = 321.0, Ti_endstep = 0, Ti_begstep = 0, Potential = 0.0, OldAcc = 0.0

Random and Conversion

julia> p = rand(PVector{Float64}, 3)3-element Vector{PVector{Float64}}:
 PVector{Float64}(0.4675296210769132, 0.14071960569245578, 0.36080797806507536)
 PVector{Float64}(0.8600082438947183, 0.11430952909041481, 0.08145226029334951)
 PVector{Float64}(0.23058725796424084, 0.20410350748520034, 0.521572995526429)
julia> pu = rand(PVector{Float64}, 3) * u"m"3-element Vector{PVector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}:
 PVector(0.985808685834883 m, 0.17405280557579272 m, 0.5779277729115215 m)
 PVector(0.281736232974783 m, 0.09004641264897861 m, 0.8514818414483671 m)
 PVector(0.20279789714935403 m, 0.590023609483166 m, 0.011192218037889057 m)
julia> p_Ball = [Ball(uSI) for i=1:3]3-element Vector{Ball{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}, Unitful.Quantity{Float64, 𝐋 𝐓^-2, Unitful.FreeUnits{(m, s^-2), 𝐋 𝐓^-2, nothing}}, Unitful.Quantity{Float64, 𝐌, Unitful.FreeUnits{(kg,), 𝐌, nothing}}, Int64}}:
 Ball 0: Pos = PVector(0.0 m, 0.0 m, 0.0 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
 Ball 0: Pos = PVector(0.0 m, 0.0 m, 0.0 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
 Ball 0: Pos = PVector(0.0 m, 0.0 m, 0.0 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
julia> assign_particles(p_Ball, :Pos, pu)
julia> p_Ball3-element Vector{Ball{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}, Unitful.Quantity{Float64, 𝐋 𝐓^-2, Unitful.FreeUnits{(m, s^-2), 𝐋 𝐓^-2, nothing}}, Unitful.Quantity{Float64, 𝐌, Unitful.FreeUnits{(kg,), 𝐌, nothing}}, Int64}}:
 Ball 0: Pos = PVector(0.985808685834883 m, 0.17405280557579272 m, 0.5779277729115215 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
 Ball 0: Pos = PVector(0.281736232974783 m, 0.09004641264897861 m, 0.8514818414483671 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
 Ball 0: Pos = PVector(0.20279789714935403 m, 0.590023609483166 m, 0.011192218037889057 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 0.0 kg
julia> pconvert([1.0, 2.0, 3.0])PVector{Float64}(1.0, 2.0, 3.0)
julia> pconvert([1.0u"m" 4.0u"m";
                 2.0u"m" 5.0u"m";
                 3.0u"m" 6.0u"m"])2-element StructArray(::Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, ::Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, ::Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}) with eltype PVector:
 PVector(1.0 m, 2.0 m, 3.0 m)
 PVector(4.0 m, 5.0 m, 6.0 m)

Extent

julia> p = [Ball(PVector(-1.0u"m", 1.0u"m", 1.0u"m"), PVector(u"m/s"), PVector(u"m/s^2"), 1.0u"kg", 1),
            Ball(PVector(3.0u"m", -3.0u"m", -3.0u"m"), PVector(u"m/s"), PVector(u"m/s^2"), 3000.0u"g", 2)]2-element Vector{Ball{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}}, Unitful.Quantity{Float64, 𝐋 𝐓^-2, Unitful.FreeUnits{(m, s^-2), 𝐋 𝐓^-2, nothing}}, M, Int64} where M}:
 Ball 1: Pos = PVector(-1.0 m, 1.0 m, 1.0 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 1.0 kg
 Ball 2: Pos = PVector(3.0 m, -3.0 m, -3.0 m), Vel = PVector(0.0 m s^-1, 0.0 m s^-1, 0.0 m s^-1), Acc = PVector(0.0 m s^-2, 0.0 m s^-2, 0.0 m s^-2), Mass = 3000.0 g
julia> minimum_x(p)-1.0 m
julia> maximum_x(p)3.0 m
julia> center(p)PVector(1.0 m, -1.0 m, -1.0 m)
julia> pos_center(p)PVector(1.0 m, -1.0 m, -1.0 m)
julia> mass_center(p)PVector(2.0 m, -2.0 m, -2.0 m)
julia> median(p, :Pos)PVector(1.0 m, -1.0 m, -1.0 m)
julia> extent(p)Extent: xMin = -1.0 m, xMax = 3.0 m, yMin = -3.0 m, yMax = 1.0 m, zMin = -3.0 m, zMax = 1.0 m, SideLength = 4.0 m, Center = PVector(1.0 m, -1.0 m, -1.0 m)