In [1]:

    

%matplotlib inline
%load_ext autoreload
%autoreload 2


    

In [2]:

    

import numpy as np
from scipy.stats import norm, rayleigh, binom
from picsheath import *
import matplotlib.pyplot as plt
from matplotlib import animation
from collections import namedtuple


    

In [3]:

    

# Choose parameters based on debye length

# Particle temperatures
KT  = 0.1
# Number of particles
N_e = 40000
# Number of cells spanned by a debye length
ld_cell = 100.0
nx = 10000
# Number of particle pairs injected per time step
n_pairs = 2
nt = 20000
L  = 2*np.pi
L_source = L/2.0

dx = L/(nx-1)
dt = 0.01
B0 = 0
me = 1.
mi = 1836.

n0  = N_e/L
ld  = dx*ld_cell
q   = np.sqrt(KT/(2*n0*ld**2))
N_i = N_e


    

In [4]:

    

n0 = N_e/L
ld = np.sqrt(KT/(n0*q**2)/2.0)
print "Debye Length: ", round(ld, 3)
print "dx:           ", round(dx, 4)
print "N cells/ld:   ", round(ld/dx, 3)
print "N_e/ld:       ", round(N_e/(L/ld))
print "q:            ", round(q, 4)
print "N:            ", N_e+N_i
print "N inject:     ", nt*n_pairs*2


    

    




Debye Length:  0.063
dx:            0.0006
N cells/ld:    100.0
N_e/ld:        400.0
q:             0.0446
N:             80000
N inject:      80000

In [5]:

    

n0e = np.linspace(0., L, N_e+1)[:-1]
n0i = np.linspace(0., L, N_i+1)[:-1]

sample_e = lambda size: norm.rvs(size=size, scale=np.sqrt(KT/me)) 
sample_i = lambda size: norm.rvs(size=size, scale=np.sqrt(KT/mi))
vx0i = sample_i(N_i)
vx0e = sample_e(N_e)

def emmert_source(size, scale):
    d = binom.rvs(1, .5, size=size)
    d[d==0] = -1
    d[d==1] = 1
    sample  = d*rayleigh.rvs(scale=scale, size=size)
    return sample

source_e = lambda size: emmert_source(size, np.sqrt(KT/me)) 
source_i = lambda size: emmert_source(size, np.sqrt(KT/mi))

electron = Species(-q, me, N_e, n0e,
                   vx0e, np.zeros(N_e),
                   sample_e, source_e)
ion      = Species(q, mi, N_i, n0i,
                   vx0i, np.zeros(N_i),
                   sample_i, source_i)

results = pic(electron, ion, nx, dx, nt, dt, L, B0,
              ['sig', 'N', 'phi_wall'],
              n_pairs=n_pairs, L_source=L_source)

#el  = results['el']
phi      = results['phi']
phi_wall = results['phi_wall']
sig      = results['sig_all']
Na       = results['N_all']
nphi_wall = q*phi_wall/KT

x_vals = np.arange(nx)*dx
t_vals = np.arange(nt+1)*dt


    

In [6]:

    

plt.plot(t_vals, sig)


    

    
Out[6]:





[<matplotlib.lines.Line2D at 0x7fe6d68f5ad0>]






    

In [7]:

    

plt.plot(x_vals, q*phi/KT)
plt.axvline(L_source, color='g', linestyle='--')
plt.axvline(L-ld, color='r', linestyle='--')
plt.xlabel("x")
plt.ylabel("$e\phi/KT$")
plt.savefig("sheath-profile.pdf")


    

In [8]:

    

i = 1000
avg_nphiw = np.mean(nphi_wall[-i:])
plt.plot(nphi_wall[-i:])
plt.ylabel("$e\phi/KT$")
plt.axhline(avg_nphiw, color='g', linestyle='--')
plt.savefig("emmert-source-phiw.pdf")
avg_nphiw


    

    
Out[8]:





-2.8500838671537156






    

In [9]:

    

plt.plot(Na)


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0x7fe6d6856f90>]






    

In [9]: