No BSD License  

Highlights from
Numerical Methods for Physics

from Numerical Methods for Physics by Alejandro Garcia
Companion Software

mover(x,v,npart, ...
function [x,v,strikes,delv] = mover(x,v,npart, ...
                                xwall,mpv,vwall,tau);
% mover - Function to move particles by free flight
%         Also handles collisions with walls
% Inputs
%    x     - Positions of the particles
%    v     - Velocities of the particles
%    npart - Number of particles in the system
%    xwall - Positions of walls
%    mpv   - Most probable velocity off the wall
%    vwall - Wall velocities
%    tau   - Time step
% Output
%    x,v     - Updated positions and velocities
%    strikes - Number of particles striking each wall
%    delv    - Change of y-velocity at each wall
%     
strikes = [0 0];  delv = [0 0];
direction = [1 -1];   % Direction of particle leaving wall
x_old = x;            % Remember original position
% Update x position of particle
x(:) = x_old(:) + v(:,1)*tau;  
%x = rem(x+xwall(2),xwall(2)); % Periodic boundary conditions
%return;
for i=1:npart
  if( x(i) <= xwall(1) )
    flag=1; % Particle strikes left wall
  elseif( x(i) >= xwall(2) )
    flag=2; % Particle strikes right wall
  else
   flag=0;  % Particle strikes neither wall
  end
  if( flag > 0 )
    strikes(flag) = strikes(flag) + 1;
    delv(flag) = delv(flag) - v(i,2);
    % Reset velocities as biased Maxwellian
    v(i,1) = direction(flag)*sqrt(-log(1-rand(1))) * mpv;
    phi = 2*pi*rand(1);
    v_parll = sqrt(-log(1-rand(1))) * mpv;
    v(i,2) = v_parll*sin(phi) + vwall(flag);
    v(i,3) = v_parll*cos(phi);
    % Time of flight after leaving wall
    dtr = tau*(x(i)-xwall(flag))/(x(i)-x_old(i));   
    % Position after leaving wall
    x(i) = xwall(flag) + v(i,1)*dtr;  
    delv(flag) = delv(flag) + v(i,2);
  end
end

Contact us at files@mathworks.com