# Metaprogramming ## The ParticleGroup example ```julia import Sobol using Plots, LinearAlgebra ``` ``` ParticleGroup{D,V}(n_particles, charge, mass) ``` * `D` : number of dimension in physical space * `V` : number of dimension in phase space * `n` : number of particles ```julia mutable struct ParticleGroup{D,V} n_particles :: Int64 data :: Array{Float64, 2} function ParticleGroup{D,V}(n) where {D, V} data = zeros( Float64, (D+V, n)) new( n, data) end end ``` --- Set position of ith particle of p to x ```julia @generated function set_x!( p :: ParticleGroup{D,V}, i, x :: Float64 ) where {D, V} :(p.data[1, i] = x) end ``` ``` set_x! (generic function with 1 method) ``` -- Set position of ith particle of p to x when x is a vector ```julia @generated function set_x!( p :: ParticleGroup{D,V}, i, x :: Vector{Float64} ) where {D, V} :(for j in 1:$D p.data[j, i] = x[j] end) end ``` ``` set_x! (generic function with 2 methods) ``` --- Set velocity of ith particle of p to v ```julia @generated function set_v!( p :: ParticleGroup{D,V}, i, v :: Float64 ) where {D, V} :(p.data[$D+1, i] = v) end ``` ``` set_v! (generic function with 1 method) ``` -- Set velocity of ith particle of p to v ```julia @generated function set_v!( p :: ParticleGroup{D,V}, i, v :: Vector{Float64} ) where {D, V} :(for j in 1:$V p.data[$D+j, i] = v[j] end) end ``` ``` set_v! (generic function with 2 methods) ``` --- Get position of ith particle of p ```julia @generated function get_x( p :: ParticleGroup{D,V}, i ) where {D, V} :(p.data[1:$D, i]) end ``` ``` get_x (generic function with 1 method) ``` Get velocity of ith particle of p ```julia @generated function get_v( p :: ParticleGroup{D,V}, i ) where {D, V} :(p.data[$D+1:$D+$V, i]) end ``` ``` get_v (generic function with 1 method) ``` --- Sampling from a probability distribution to initialize Landau damping ```julia function landau_sampling!( pg :: ParticleGroup{1,2}, alpha, kx ) function newton(r) x0, x1 = 0.0, 1.0 r *= 2π / kx while (abs(x1-x0) > 1e-12) p = x0 + alpha * sin( kx * x0) / kx f = 1 + alpha * cos( kx * x0) x0, x1 = x1, x0 - (p - r) / f end x1 end s = Sobol.SobolSeq(2) nbpart = pg.n_particles for i=1:nbpart v = sqrt(-2 * log( (i-0.5)/nbpart)) r1, r2 = Sobol.next!(s) θ = r1 * 2π set_x!(pg, i, newton(r2)) set_v!(pg, i, [v * cos(θ), v * sin(θ)]) end end ``` ``` landau_sampling! (generic function with 1 method) ``` --- ```julia n_particles = 10000 pg = ParticleGroup{1,2}( n_particles) alpha, kx = 0.1, 0.5 landau_sampling!(pg, alpha, kx) ``` -- ```julia xp = vcat([get_x(pg, i) for i in 1:pg.n_particles]...) vp = vcat([get_v(pg, i) for i in 1:pg.n_particles]'...) ``` ``` 10000×2 Array{Float64,2}: -4.4505 5.45029e-16 -7.70866e-16 -4.1964 2.4939e-16 4.07285 -2.82092 2.82092 2.77602 -2.77602 -2.73963 -2.73963 2.70897 2.70897 1.45172 3.50477 -1.43904 -3.47415 3.44671 -1.42768 ⋮ 0.0146796 0.0385387 -0.0137887 -0.0361999 0.0336995 -0.0128363 -0.0309982 0.0118074 -0.0122742 0.0273778 0.0108246 -0.0241444 -0.0204052 -0.0091482 0.0158054 0.00708599 0.00815159 0.00579259 ``` --- ```julia pp = plot(layout=(3,1)) histogram!(pp[1,1], xp, normalize=true, bins = 100, lab=:x) plot!(pp[1,1], x -> (1+alpha*cos(kx*x))/(2π/kx), 0., 2π/kx, lab="") histogram!(pp[2,1], vp[:,1], normalize=true, bins = 100, lab=:vx) plot!(pp[2,1], v -> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="") histogram!(pp[3,1], vp[:,2], normalize=true, bins = 100, lab=:vy) plot!(pp[3,1], v -> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="") ``` ![](particles.svg) --- ```julia histogram2d(vp[:,1], vp[:,2], normalize=true, bins=100) ``` ![](hist2d.svg) Next: [GEMPIC](/Numkin2019/05/build/index.html) *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*