-- Gkyl ------------------------------------------------------------------------
-- 1x1v simulation of collisionless damping of an electron Langmuir wave,
-- using kinetic ions and electrons. We use normalize units.
--------------------------------------------------------------------------------
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()

permitt  = 1.0   -- Permittivity of free space.
permeab  = 1.0   -- Permeability of free space.
eV       = 1.0   -- Elementary charge, or Joule-eV conversion factor.
elcMass  = 1.0   -- Electron mass.
ionMass  = 1.0   -- Ion mass.

nElc = 1.0    -- Electron number density.
nIon = nElc   -- Ion number density.
Te   = 1.0    -- Electron temperature.
Ti   = Te     -- Ion temperature.

vtElc   = math.sqrt(eV*Te/elcMass)                   -- Electron thermal speed.
vtIon   = math.sqrt(eV*Ti/ionMass)                   -- Ion thermal speed.
wpe     = math.sqrt((eV^2)*nElc/(permitt*elcMass))   -- Plasma frequency.
lambdaD = vtElc/wpe                                  -- Debye length.

-- Amplitude and wavenumber of sinusoidal perturbation.
pertA = 1.0e-3
pertK = .750/lambdaD

-- Maxwellian in (x,vx)-space, given the density (denU), bulk flow
-- velocity (flowU), mass and temperature (temp).
local function maxwellian1D(x, vx, den, flowU, mass, temp)
   local v2   = (vx - flowU)^2
   local vtSq = temp/mass
   return (den/math.sqrt(2*math.pi*vtSq))*math.exp(-v2/(2*vtSq))
end

plasmaApp = Plasma.App {
   tEnd         = 20.0/wpe,           -- End time.
   nFrame       = 20,                 -- Number of output frames.
   lower        = {-math.pi/pertK},   -- Lower boundary of configuration space.
   upper        = { math.pi/pertK},   -- Upper boundary of configuration space.
   cells        = {64},               -- Configuration space cells.
   polyOrder    = 1,                  -- Polynomial order.
   periodicDirs = {1},                -- Periodic directions.
   
   elc = Plasma.Species {
      charge = -eV, mass = elcMass,
      lower = {-6.0*vtElc},      -- Velocity space lower boundary.
      upper = { 6.0*vtElc},      -- Velocity space upper boundary.
      cells = {64},              -- Number of cells in velocity space.
      init = function (t, xn)    -- Initial conditions.
	 local x, v = xn[1], xn[2]
	 return (1+pertA*math.cos(pertK*x))*maxwellian1D(x, v, nElc, 0.0, elcMass, Te) 
      end,
      evolve = true, -- Evolve species?
   },

   ion = Plasma.Species {
      charge = eV, mass = ionMass,
      lower = {-6.0*vtIon},      -- Velocity space lower boundary.
      upper = { 6.0*vtIon},      -- Velocity space upper boundary.
      cells = {64},              -- Number of cells in velocity space.
      init  = function (t, xn)   -- Initial conditions.
	 local x, v = xn[1], xn[2]
	 return maxwellian1D(x, v, nIon, 0.0, ionMass, Ti) 
      end,
      evolve = true, -- Evolve species?
   },

   field = Plasma.Field {
      epsilon0 = permitt, mu0 = permeab,
      init = function (t, xn)   -- Initial conditions.
         local Ex, Ey, Ez = -pertA*math.sin(pertK*xn[1])/pertK, 0.0, 0.0
         local Bx, By, Bz = 0.0, 0.0, 0.0
         return Ex, Ey, Ez, Bx, By, Bz
      end,
      evolve = true, -- Evolve field?
   },
}
-- Run application.
plasmaApp:run()