function m671b_Vlasov-Poisson_PIC_hw3
% 
% PIC code for 1-D electron plasma two-stream instability
%
  clear; clf;

  q = -1;                 % electron charge
  Np = 200;            % number of particles
  Nm = 100;           % number of mesh cells
  dx = 1/Nm;          % mesh spacing
  xm = 0:dx:1;        % mesh points
  rho_0 = -q*2*Np; % background neutralizing charge density
  
  t_final = 1; dt = 0.001; Nstep = t_final/dt;   
  vmin = -3; vmax = 3;
  
  xp1 = zeros(1,Np); vp1 = zeros(1,Np); Fp1 = zeros(1,Np); 
  xp2 = zeros(1,Np); vp2 = zeros(1,Np); Fp2 = zeros(1,Np);
  phi = zeros(1,Nm+1); rho = zeros(1,Nm+1); E = zeros(1,Nm+1);
  
% initial data

  for i = 1:Np
    xp_beam = (i-0.5)/Np;
    xp1(i) = xp_beam + 0.02*sin(6*pi*xp_beam);
    vp1(i) = 1;
    xp2(i) = xp_beam + 0.02*sin(6*pi*xp_beam); 
    vp2(i) = -1;
  end
  
% map particles back to fundamental period

  for i = 1:Np
      if xp1(i) > 1; xp1(i) = xp1(i)-1; end
      if xp1(i) < 0; xp1(i) = xp1(i)+1; end
      if xp2(i) > 1; xp2(i) = xp2(i)-1; end
      if xp2(i) < 0; xp2(i) = xp2(i)+1; end
  end 
  
  xp1_old = xp1; vp1_old = vp1;
  xp2_old = xp2; vp2_old = vp2;
  
  str = sprintf('t = %0.5g',0);
  figure(1); plot(xp1,vp1,'bo',xp2,vp2,'ro'); 
  xlabel('x'); ylabel('v'); title(str); axis([ 0 1 vmin vmax])

% timestepping

  for step = 1:Nstep

% assign charge from particles to mesh points

  n = zeros(1,Nm+1);
  for j = 1:Nm+1
    n(j) = 0; 
    for i = 1:Np
      n(j) = n(j) + W(xp1(i)-xm(j),dx);
      n(j) = n(j) + W(xp2(i)-xm(j),dx);
    end
    rho(j) = q*n(j)/dx + rho_0;
  end
  rho(1) = q*(n(1)+n(Nm+1))/dx + rho_0;
  rho(Nm+1) = rho(1);
  
  str = sprintf('t = %0.5g',step*dt);
  figure(1); plot(xp1,vp1,'bo',xp2,vp2,'ro'); 
  xlabel('x'); ylabel('v'); title(str); axis([ 0 1 vmin vmax])

% solve Poisson equation for potential on the mesh
% input : rho(j), j = 1:Nm+1
% output : phi(j), j = 1:Nm+1

% determine electric field on the mesh

  E(1) = -(phi(2)-phi(Nm))/(2*dx);
  for j = 2:Nm
    E(j) = -(phi(j+1)-phi(j-1))/(2*dx);
  end
  E(Nm+1) = E(1);

% interpolate forces from mesh points to particles

  for i = 1:Np
    Fp1(i) = 0; Fp2(i) = 0;
    for j = 1:Nm+1
      Fp1(i) = Fp1(i) + W(xp1(i)-xm(j),dx)*E(j);
      Fp2(i) = Fp2(i) + W(xp2(i)-xm(j),dx)*E(j);
    end
    Fp1(i) = q * Fp1(i); Fp2(i) = q * Fp2(i);
  end

% update particle velocities and advect particles

  for i = 1:Np
    vp1(i) = vp1_old(i) + dt * Fp1(i);
    xp1(i) = xp1_old(i) + dt * vp1(i);
    vp2(i) = vp2_old(i) + dt * Fp2(i);
    xp2(i) = xp2_old(i) + dt * vp2(i);
  end

% map particles back to fundamental period

  for i = 1:Np
      if xp1(i) > 1; xp1(i) = xp1(i)-1; end
      if xp1(i) < 0; xp1(i) = xp1(i)+1; end
      if xp2(i) > 1; xp2(i) = xp2(i)-1; end
      if xp2(i) < 0; xp2(i) = xp2(i)+1; end
  end     
  xp1_old = xp1; vp1_old = vp1; xp2_old = xp2; vp2_old = vp2;

  end
  str = sprintf('t = %0.5g',Nstep*dt);
  figure(1); plot(xp1,vp1,'bo',xp2,vp2,'ro'); 
  xlabel('x'); ylabel('v'); title(str); axis([ 0 1 vmin vmax])
  
  function weight = W(x,dx)
  weight = 0;
% nearest mesh point
% if abs(x) < dx/2; weight = 1; end 
% cloud-in-cell
  if abs(x) < dx; weight = 1-abs(x)/dx; end
  
 
