Methods

function run(sim::Simulation)

This function does all the work for you:

1. prepare sim environment and preprocess data
2. evaluate force, initialize timesteps, output
3. mainloop

softlen(c::Collection, table::MVector)
softlen(p::AbstractParticle, table::MVector)
softlen(p::AbstractParticle, sim::Simulation)
softlen(c::Collection, sim::Simulation)

Retrive softening length from table

suggest_softlen(V::Number, N::Int64)
suggest_softlen(data::Union{Array, StructArray})
suggest_softlen(data::Union{Array, StructArray}, collection::Collection)
suggest_softlen(sim::Simulation)

return recommended softening length

See also suggest_softlen! and set_softlen!

suggest_softlen!(sim::Simulation)

Set all softening lengths to suggested value.

smooth_coef(x_coord, fitting_order, diff_order)

Fit the data with fitting_order polynomial and calculate the diff_order differential at x=0. Input the x coordinates of data. Return the coefficients for calculating the result from data.

Suppose the data is [u,v,w], we want to smooth v by fitting a line (1 order polynomial) and use it estimate a better v. Then set the x coordinates of [u,v,w] be [-1,0,1], fit it and calculate the result (0 order differential) at x=0, the result will be (u+v+w)/3. So smooth_coef([-1,0,1],1,0) will return [1/3,1/3/1/3].

Using Rational type or Sym type of SymPy can get the exact coefficients. For example, smooth_coef(Sym[0,1,2],2,1) get Sym[-3/2, 2, -1/2], which means the first order differential of [u,v,w] at u is -1.5u+2v-0.5w. Using this way can generate all data in https://en.wikipedia.org/wiki/Finite_difference_coefficient

Examples

• Linear extrapolation: smooth_coef([1,2],1,0) get [2, -1], so the left linear extrapolation of [u,v] is 2u-v.
• Quadratic extrapolation: smooth_coef([-3,-2,-1],2,0) get [1, -3, 3], so the right Quadratic extrapolation of [u,v,w] is u-3v+3w.
• First order central differential: smooth_coef(Rational[-1,0,1],2,1) get [-1//2, 0//1, 1//2], so the first order differential of [u,v,w] at v is (w-u)/2
• Second order central differential: smooth_coef([-1,0,1],2,2) get [1, -2, 1], so the second order differential of [u,v,w] at v is u+w-2v
• Five points quadratic smoothing: smooth_coef([-2,-1,0,1,2],2,0) get [-3/35, 12/35, 17/35, 12/35, -3/35], so [-3/35 12/35 17/35 12/35 -3/35]*[a,b,c,d,e] get the smoothed c.
• Four points quadratic interpolation: smooth_coef([-3,-1,1,3],2,0) get [-1/16, 9/16, 9/16, -1/16], so smooth_coef([-3,-1,1,3],2,0)'*[a,b,c,d] get the estimated value at the middle of b and c.

diff_mat(n, order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, lpoints=div(points,2), rpoints=points-lpoints-1, fitting_order=lpoints+rpoints, boundary=:Extrapolation, boundary_points=lpoints+rpoints+1, boundary_order=boundary_points-1, sparse=false)

Generate differential matrix.

Arguments

• n: Matrix size. If v is length n vector, diff_mat(n)*v calculate the differential of v.
• order: Differential order.
• T: Matrix element type. If set T=Rational or using SymPy and set T=Sym, diff_mat will return the exact value.
• dt: Numerical differential step size.
• points: Number of points for fitting polynomial to estimate differential. This argument is only for convenience. The real number of points is always lpoints+rpoints+1.
• lpoints: Number of points at left to the target point, which differential is calculated by the fitted polynomial.
• rpoints: Number of points at right to the target point. If lpoints==rpoints, then the differential is estimated as central finite difference. If lpoints==0, then it is normal forward finite difference. If rpoints==0, then it is backward finite difference.
• fitting_order: The order of the fitted polynomial for estimating differential.
• boundary: Boundary condition. Can be Dirichlet()(boundary value is zero), Periodic(assume data is periodic), :Extrapolation(boundary value is extrapolated according to boundary_points and boundary_order), :None(not deal with boundary, will return non-square matrix).
• boundary_points: Number of points for fitting polynomial to estimate differential at boundary. Normally it should not be much less than points, otherwise sometimes the current point may not be used to estimate the differential.
• boundary_order: The order of the fitted polynomial for points determined by boundary_points.
• sparse: If true, return sparse matrix instead of dense one.

Examples

k=5; x=rand(k);
diff_mat(k,1;points=3)*x #do 3 points 1st order central differential ((x[n+1]-x[n-1])/2).
diff_mat(k,2;points=3)*x #do 3 points 2nd order central differential (x[n+1]+x[n-1]-2x[n]).
diff_mat(k,1;points=2,lpoints=0)*x #do the normal 1st order forward differential (x[n+1]-x[n]).
diff_mat(k,1;lpoints=1,rpoints=0)*x #do the 1st order backward differential (x[n-1]-x[n]).

total_angular_momentum(sim::Simulation) -> Any

Compute total angular momentum vector of the system

total_potential(data::StructArray) -> Any

Sum potential energy of particles in data. Potentials need to be computed in advance.

Return nothing if empty.

total_momentum(sim::Simulation, axis::Symbol) -> Any

Compute momentum of the system in the direction of axis

total_momentum(sim::Simulation) -> Any

Compute total momentum vector of the system

total_kinetic(data::StructArray) -> Any

Sum kinetic energy: 0.5 * data.Mass[i] * data.Vel[i] * data.Vel[i]

total_kinetic(sim::Simulation) -> Any

Compute kinetic energy of particles on workers and return the sum

total_energy(sim::Simulation) -> Any

Total energy (kinetic + potential) of sim data

nu(y::Number, n::Number = 2)

Generalized MOND (MOdified Newtonian Dynamics) interpolation function. See also nu1, nu2

nu1(y::Number)

MOND (MOdified Newtonian Dynamics) interpolation function with index = 1. See also nu, nu2

nu1(y::Number)

MOND (MOdified Newtonian Dynamics) interpolation function with index = 2. See also nu, nu1

mond_Milgrom1983(sim::Simulation, data::StructArray)

Apply Milgrom 1983 formula of MOND (MOdified Newtonian Dynamics) to accelerations

QUMOND_PDM_density(m::MeshCartesianStatic, ACC0::Number)

Compute ρPDM on the RHS (right hand side) of QUMOND (QUasi-linear MOdified Newtonian Dynamics). Return ρPDM

QUMOND_phi(m::MeshCartesianStatic, ACC0::Number, G::Number)

First compute ρ_PDM, then solve QUMOND (QUasi-linear MOdified Newtonian Dynamics) equation on the mesh. Return modified potential

QUMOND_acc(m::MeshCartesianStatic, ACC0::Number, G::Number)

1. Compute ρ_PDM
2. Solve modified potential on the mesh
3. Compute acceleration by finite differencing the potential

Return acceleration

QUasi-linear MOdified Newtonian Dynamics

r_g(G, M, c)

r_g = 2 * G * M / c^2

where G is the gravity constant, M is the mass of the central object, c is light speed.

Schwarzchild radius.

r_g(G, M, c)

r_g = 2 * G * M / c^2

where G is the gravity constant, M is the mass of the central object, c is light speed.

Schwarzchild radius.

r_g(G, M, c)

r_g = 2 * G * M / c^2

where G is the gravity constant, M is the mass of the central object, c is light speed.

Schwarzchild radius.

pseudoNewtonianPotential(G, M, c, R)

pot = - G * M / (R - r_g(G, M, c))

Pseudo-Newtonian potential around a compat central object. Diverge at gravity radius r_g.

pseudoNewtonianPotential(G, M, rg)

pot = - G * M / (R - rg)

Pseudo-Newtonian potential around a compat central object. Diverge at gravity radius r_g.

pseudoNewtonianAcc(G, M, c, R, n::AbstractPoint)

acc = - G * M / (R - r_g(G, M, c))^2 * n

Pseudo-Newtonian acceleration around a compat central object. Diverge at gravity radius r_g.

pseudoNewtonianAcc(G, M, c, R, n::AbstractPoint)

acc = - G * M / (R - rg)^2 * n

Pseudo-Newtonian acceleration around a compat central object. Diverge at gravity radius r_g.

ellipticSemiMajor(G::Number, M::Number, r::Number, v::Number) -> Any

Length of semi-major axis of elliptic orbit

• G: gravitational constant
• M: mass of central object
• r: orbit radius at specific time
• v: velocity at radius r

ellipticPeriod(G::Number, M::Number, a::Number) -> Any

Orbital period of elliptic orbit

• G: gravitational constant
• M: mass of central object
• a: length of semi-major axis

eccentricity(G::Number, M::Number, v::Number, h::Number) -> Any

Eccentricity of elliptic orbit

• G: gravitational constant
• M: mass of central object
• v: velocity at radius r
• h: specific relative angular momentum at radius r

eccentricity(G::Number, M::Number, r::AbstractPoint, v::AbstractPoint) -> Any

Eccentricity vector of elliptic orbit

• G: gravitational constant
• M: mass of central object
• r: position at specific time
• v: velocity vector at position r

typicalvelocity(G, M, R)

v = sqrt(G*M/R)

Typical velocity of a spherical system.

meandensity(M, R)

ρ = 3*M/(4πR^3)

Mean density of a spherical system.

crosstime(R, v)

t_cross = R / v

The typical time needed to cross the system.

hubbletime(H0)

t_hubble = 1/H0

The age of the universe.

relaxtime(R, v, N)

trelax = N/(10 lnN) tcross

Relaxation time: The time over which the change in kinetic energy due to the long-range collisions has accumulated to a value that is comparable to the intrinsic kinetic energy of the particle.

interactiontime(R, v, N)

The typical time between two short-range interactions that cause a change in kinetic energy comparable to the intrinsic kinetic energy of the particle.

tinteraction = N * tcross

dynamicaltime(ρ, G)

t_dyn = sqrt(3 * π / (16 * G * ρ))

The time required to travel halfway across the system.

freefalltime(ρ, G)

t_ff = sqrt(3 * π / (32 * G * ρ))

The time it takes a sphere with zero pressure to collapse to a point.

orbitaltime(ρ, G)

t_orb = sqrt(3π/(Gρ))

The time it takes to complete a (circular) orbit.

struct LogInfo

• timers::Symbol

  Timer enum names to access timing (continuously starting from 1)

• timing::Vector{UInt64}

  Timings defined by timer enums

• analysers::Dict{String, Function}

  Analyse on the whole simulation

Examples #TODO

Use uppercase letters to avoid

Better printing of trait types

function writegadget2param(filename::String, ICfilename::String)

Params

• filename : name of param file
• ICfilename : name of initial conditions

alter_param(filename::String, param::String, value)

alter all lines containing param to param value

alter_line(filename::String, match::String, newline::String; all = false)

alter the first line containing match to newline

Keywords

• all::Bool = false : If true, alter all lines containing match to newline