## PhysicalParticles.jl

Physical vector and particle types for Julia
Author JuliaAstroSim
Popularity
2 Stars
Updated Last
2 Years Ago
Started In
November 2019

# PhysicalParticles.jl

Physical particle types for scientific simulation. Manipulate vectors as simple as numbers!

## Installation

`]add PhysicalParticles`

or

`using Pkg; Pkg.add("PhysicalParticles")`

or

`using Pkg; Pkg.add("https://github.com/JuliaAstroSim/PhysicalParticles.jl")`

To test the Package:

`]test PhysicalParticles`

## Documentation

• Devdocumentation of the in-development version.

## 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) * im
PVector2D{Complex{Float64}}(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") / 2
PVector(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)
1.5428387907456837e20 m

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}(-8.659560562354932e-17, -2.220446049250313e-16, 0.9999999999999998)

# Coordinate Transformations
julia> cylinderial(PVector(sqrt(2), sqrt(2), 1.0, u"m"))
(2.0 m, 0.7853981633974484, 1.0 m)

julia> cylinderial2xyz(2.0u"m", pi/4, 1.0u"m")
PVector(1.4142135623730951 m, 1.414213562373095 m, 1.0 m)

julia> spherical(PVector(sqrt(0.5), sqrt(0.5), 1.0, u"m"))
(1.4142135623730951 m, 0.7853981633974484, 0.7853981633974483)

julia> spherical2xyz(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> SPHGas()
SPHGas 0 GAS: 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, Entropy = 0.0, Density = 0.0, Hsml = 0.0, Left = 0.0, Right = 0.0, NumNgbFound = 0, RotVel = PVector{Float64}(0.0, 0.0, 0.0), DivVel = 0.0, CurlVel = 0.0, dHsmlRho = 0.0, Pressure = 0.0, DtEntropy = 0.0, MaxSignalVel = 0.0

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 M⊙ Gyr^-2, OldAcc = 0.0 kpc Gyr^-2

julia> b = SPHGas(uAstro)
SPHGas 0 GAS: 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 M⊙ Gyr^-2, OldAcc = 0.0 kpc Gyr^-2, Entropy = 0.0 kpc^2 M⊙ Gyr^-2 K^-1, Density = 0.0 M⊙ kpc^-3, Hsml = 0.0 kpc, Left = 0.0, Right = 0.0, NumNgbFound = 0, RotVel = PVector(0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1, 0.0 kpc Gyr^-1), DivVel = 0.0 Gyr^-1, CurlVel = 0.0 Gyr^-1, dHsmlRho = 0.0 kpc, Pressure = 0.0 M⊙ Gyr^-2 kpc^-1, DtEntropy = 0.0 kpc^2 M⊙ Gyr^-3 K^-1, MaxSignalVel = 0.0 kpc Gyr^-1

julia> distance(a,b)
0.0 m```

### StructArrays.jl support

`StructArray` provides a more efficient way to iterate on a field of particles:

```sArray = [Star() for i in 1:5]
sStruct = StructArray(sArray)

# Easier to set properties, and even faster!
sStruct.Mass[1] = 1000.0

assign_particles(sStruct, :Pos, randn_pvector(5))

mean(sStruct.Pos)```

### Random and Conversion

```julia> p = rand_pvector(3)
3-element Array{PVector{Float64},1}:
PVector{Float64}(0.899541890819791, 0.49609709458549345, 0.22817220536717397)
PVector{Float64}(0.21907343513386301, 0.39110699072427035, 0.3502946880565312)
PVector{Float64}(0.8107782153679699, 0.20218167820102884, 0.94236923352867)

julia> pu = rand_pvector(3, u"m")
3-element Array{PVector{Quantity{Float64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}}},1}:
PVector(0.5346672699901402 m, 0.6988269071898365 m, 0.8120077168096169 m)
PVector(0.46886820909936744 m, 0.9575982422487646 m, 0.10413358701332642 m)
PVector(0.0219005354136228 m, 0.327612194392396 m, 0.2837471711064179 m)

julia> p_Ball = [Ball(uSI) for i=1:3]
3-element Array{Ball{Int64},1}:
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_points(p_Ball, :Pos, pu)

julia> p_Ball
3-element Array{Ball{Int64},1}:
Ball 0: Pos = PVector(0.5346672699901402 m, 0.6988269071898365 m, 0.8120077168096169 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.46886820909936744 m, 0.9575982422487646 m, 0.10413358701332642 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.0219005354136228 m, 0.327612194392396 m, 0.2837471711064179 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 Array{PVector,1}:
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 Array{Ball{Int64},1}:
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 = 1000.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)```

There are differences among `center`, `pos_center`, `mass_center` and `median`:

• `center`: box center of particles
• `pos_center`: average position of particles
• `mass_center`: mass weighted average position of particles
• `median`: middle value of positions of particles

### Units

Units are supported by Unitful.jl and UnitfulAstro.jl

Set default units by

```const uSI = u"m,s,A,K,cd,kg,mol"
preferunits(uSI)```

or simply call `si()`. `astro()` and `cgs()` are implemented in the same way.

This would affect unit promotions in `Unitful` package and default outputs in related packages, by setting `Unitful.promotion` and `PhysicalParticles.uDefaults` respectively.

Interfaces to get basic units:

``````julia> getunits()
(m, s, A, K, cd, kg, mol)

julia> getunits(uAstro)
(kpc, Gyr, A, K, cd, M⊙, mol)

julia> getunits(nothing)
(nothing, nothing, nothing, nothing, nothing, nothing, nothing)

julia> getuLength()
m

julia> getuTime(uSI)
s

julia> getuCurrent(uCGS)
A

julia> getuTemperature(nothing)

julia> getuLuminosity()
cd

julia> getuMass()
kg

julia> getuAmount()
mol
``````

`axisunit` provides a convenient way to print units in the axis of plots:

```julia> axisunit(nothing)
""

julia> axisunit(u"m")
" [m]"

julia> axisunit("Time", u"Gyr")
"TIme [Gyr]"```

### Constants

Physical constants are imported from `CODATA2018` supported by PhysicalConstants.jl. However, constants in `PhysicalConstants` may cause type error if they are not converted to default units.

To prevent this problem, construct an immutable struct `Constant` corresponding to the provided `units`:

```julia> Constant()
Converted Constants:
G = 4.498502151469553e-6 kpc^3 M⊙^-1 Gyr^-2
m_e = 4.581240435253955e-61 M⊙
m_n = 8.423451938769546e-58 M⊙
m_p = 8.411856872862986e-58 M⊙
k_B = 7.2624677363918e-60 kpc^2 M⊙ K^-1 Gyr^-2
ACC0 = 3872.920970357523 kpc Gyr^-2

julia> Constant(uSI)
Converted Constants:
G = 6.6743e-11 m^3 kg^-1 s^-2
m_e = 9.1093837015e-31 kg
m_n = 1.67492749804e-27 kg
m_p = 1.67262192369e-27 kg
k_B = 1.380649e-23 kg m^2 K^-1 s^-2
ACC0 = 1.2e-10 m s^-2

julia> Constant(uCGS)
Converted Constants:
G = 6.674299999999999e-8 cm^3 g^-1 s^-2
m_e = 9.1093837015e-28 g
m_n = 1.67492749804e-24 g
m_p = 1.67262192369e-24 g
k_B = 1.380649e-16 g cm^2 K^-1 s^-2
ACC0 = 1.2e-8 cm s^-2

julia> using Unitful

julia> ustrip(Constant())
Converted Constants:
G = 4.498502151469553e-6
m_e = 4.581240435253955e-61
m_n = 8.423451938769546e-58
m_p = 8.411856872862986e-58
k_B = 7.2624677363918e-60
ACC0 = 3872.920970357523```

### Zerovalues

`ZeroValue` is useful for accumulated summation, array initialization, etc. Examples:

```ZeroValue(nothing)
ZeroValue()
ZeroValue(uSI)
ZeroValue(uCGS)```