Title

# 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).*