% dsmcne - Program to simulate a dilute gas using DSMC algorithm
% This version simulates planar Couette flow
clear all;  help dsmcne;   % Clear memory and print header

%* Initialize constants  (particle mass, diameter, etc.)
boltz = 1.3806e-23;    % Boltzmann's constant (J/K)
mass = 6.63e-26;       % Mass of argon atom (kg)
diam = 3.66e-10;       % Effective diameter of argon atom (m)
T = 273;               % Initial temperature (K)
density = 2.685e25;    % Number density of argon at STP (m^-3)
L = 1e-6;              % System size is one micron
Volume = L^3;          % Volume of the system (m^3)
npart = input('Enter number of simulation particles: ');
eff_num = density*Volume/npart;
fprintf('Each simulation particle represents %g atoms\n',eff_num);
mfp = Volume/(sqrt(2)*pi*diam^2*npart*eff_num);
fprintf('System width is %g mean free paths \n',L/mfp);
mpv = sqrt(2*boltz*T/mass);  % Most probable initial velocity 
vwall_m = input('Enter wall velocity as Mach number: ');
vwall = vwall_m * sqrt(5/3 * boltz*T/mass);
fprintf('Wall velocities are %g and %g m/s \n',-vwall,vwall);

%* Assign random positions and velocities to particles
rand('state',1);        % Initialize random number generators
randn('state',1);
x = L*rand(npart,1);    % Assign random positions
% Assign thermal velocities using Gaussian random numbers
v = sqrt(boltz*T/mass) * randn(npart,3);
% Add velocity gradient to the y-component
v(:,2) = v(:,2) + 2*vwall*(x(:)/L) - vwall;

%* Initialize variables used for evaluating collisions
ncell = 20;                   % Number of cells
tau = 0.2*(L/ncell)/mpv;      % Set timestep tau
vrmax = 3*mpv*ones(ncell,1);  % Estimated max rel. speed in a cell
selxtra = zeros(ncell,1);     % Used by collision routine "colider"
coeff = 0.5*eff_num*pi*diam^2*tau/(Volume/ncell);

%* Declare structure for lists used in sorting
sortData = struct('ncell', ncell,    ...
                  'npart', npart,    ...
                  'cell_n', zeros(ncell,1),  ...
                  'index', zeros(ncell,1),    ...
                  'Xref', zeros(npart,1));  

%* Initialize structure and variables used in statistical sampling
sampData = struct('ncell', ncell,    ...
                  'nsamp', 0,    ...
                  'ave_n', zeros(ncell,1), ...
                  'ave_u', zeros(ncell,3), ...
                  'ave_T', zeros(ncell,1));
tsamp = 0;                    % Total sampling time
dvtot = zeros(1,2);           % Total momentum change at a wall
dverr = zeros(1,2);           % Used to find error in dvtot

%* Loop for the desired number of time steps
colSum = 0;  strikeSum = [0 0];
nstep = input('Enter total number of timesteps: ');
for istep = 1:nstep
	
  %* Move all the particles
  [x, v, strikes, delv] = mover(x,v,npart,L,mpv,vwall,tau);
  strikeSum = strikeSum + strikes;

  %* Sort the particles into cells
  sortData = sorter(x,L,sortData);

  %* Evaluate collisions among the particles
  [v, vrmax, selxtra, col] = ...
          colider(v,vrmax,tau,selxtra,coeff,sortData);
  colSum = colSum + col;

  %* After initial transient, accumulate statistical samples
  if(istep > nstep/10) 
    sampData = sampler(x,v,npart,L,sampData);
    dvtot = dvtot + delv;
    dverr = dverr + delv.^2;
    tsamp = tsamp + tau;
  end

  %* Periodically display the current progress
  if( rem(istep,10) < 1 )
    fprintf('Finished %g of %g steps, Collisions = %g\n', ...
                                            istep,nstep,colSum);
    fprintf('Total wall strikes: %g (left)  %g (right)\n', ...
                                  strikeSum(1),strikeSum(2));
  end
end

%* Normalize the accumulated statistics
nsamp = sampData.nsamp;
ave_n = (eff_num/(Volume/ncell))*sampData.ave_n/nsamp;    
ave_u = sampData.ave_u/nsamp;
ave_T = mass/(3*boltz) * (sampData.ave_T/nsamp);
dverr = dverr/(nsamp-1) - (dvtot/nsamp).^2;
dverr = sqrt(dverr*nsamp);

%* Compute viscosity from drag force on the walls
force = (eff_num*mass*dvtot)/(tsamp*L^2);
ferr = (eff_num*mass*dverr)/(tsamp *L^2);
fprintf('Force per unit area is \n');
fprintf('Left wall:   %g +/- %g \n',force(1),ferr(1));  
fprintf('Right wall:  %g +/- %g \n',force(2),ferr(2));
vgrad = 2*vwall/L;  % Velocity gradient
visc = 1/2*(-force(1)+force(2))/vgrad;  % Average viscosity
viscerr = 1/2*(ferr(1)+ferr(2))/vgrad;  % Error
fprintf('Viscosity = %g +/- %g N s/m^2\n',visc,viscerr);
eta = 5*pi/32*mass*density*(2/sqrt(pi)*mpv)*mfp;
fprintf('Theoretical value of viscoisty is %g N s/m^2\n',eta);

%* Plot average density, velocity and temperature
figure(1); clf;
xcell = ((1:ncell)-0.5)/ncell * L;
plot(xcell,ave_n); xlabel('position');  ylabel('Number density');
figure(2); clf;
plot(xcell,ave_u); xlabel('position');  ylabel('Velocities');
legend('x-component','y-component','z-component');
figure(3); clf;
plot(xcell,ave_T); xlabel('position');  ylabel('Temperature');