Computer Science COMP 199 (Winter term)
Excursions in Computing Science

---

%MATLABpak09cII:
%
%	2D kinetic gas simulator and support programs in order of appearance
%	(Anything ending Sim.m is a modification of collSim.m: compare!)
%
%	molCollWallCountBot.m			moleculeCollide.m
%	molCollWallWhen.m			
%
%	condif.m				implies.m
%	nonnegrows.m				minrows.m
%
%	moleculeCollWhen.m			moleculeInit.m
%
%	genDirection.m				genNormal.m
%
%	collSim.m
%
%	molDatHisto.m				statSim.m
%
%	momenCumHisto.m				molDataDistr.m
%
%	cumuSim.m
%
%	FTfluct.m				FTfluctAgg.m
%	momAvgHisto.m				FTfluctVar.m
%
%	manyInits.m
%
%	energyHisDis.m				entrSim.m
%
%	molCollWallWhen.m (reprise)		molCollWallCountBot.m (reprise)
%
%	tmptureSim.m
%
%	smooth2enerIgn.m
%
%	presSim.m				press2Sim.m
%
%	chmpotSim.m
%

% function p = molCollWallCountBot(pb,1,1,abcd,r,c,0)		THM	101025
% given vectors momentum pb, box boundaries abcd, position c; scalar radius r
% find new momentum p
function p = molCollWallCountBot(pb,tmult,bmult,abcd,r,c,cntB)
% tmult = 1; bmult = 1; cntB = 0;		% for future refinement
abr = abcd - [-r,r,-r,r];	% extend wall to _centre_ of molecule
p = pb;				% initialize: momentum out = momentum in
if c(1) < abr(1) | c(1) > abr(2) | c(2) < abr(3) | c(2) > abr(4)
  display(['molCollWallCountBot: ' num2str(c') ' out of box ' num2str(abr)])
elseif (c(1)==abr(1) | c(1)==abr(2)) & (c(2)==abr(3) | c(2)==abr(4)), p = -pb;
					% 4 corners: reverse momentum
elseif c(1)==abr(1) | c(1)==abr(2), p(1) = -pb(1);%side wall: reverse x-momentum
elseif c(2) == abr(4), p(2) = -pb(2);	%top: reverse y-momentum
elseif c(2) == abr(3), p(2) = -pb(2);	%bot: reverse y-momentum
else display('molCollWallCountBot: no wall collision')
end %if

% function [p1,p2] = moleculeCollide(m1,m2,p1b,p2b,0,C12)  THM		101025
% input masses (units don't matter: only ratios used)
%	initial momenta (typically opposite signs): column vectors
%	C12 is column vector connecting centres at collision time = c2 - c1
% output momentum column vectors
function [p1,p2] = moleculeCollide(m1,m2,p1b,p2b,mhf,C12)
mhf = 0;		% for future refinement
pc = p1b + p2b;		% momentum of CoM
a = m1/(m1+m2);		% dimensionless ratio for conversion to CoM
p1bc = p1b - a*pc;	% CoM momentum of first molecule
p2bc = p2b - (1-a)*pc;	% CoM momentum of second molecule	% p2bc = - p1bc
c12 = C12/sqrt(C12'*C12);	% normalize centre-to-centre vector
dim=size(c12);
reflect = 2*c12*c12' - eye(dim(1));	% reflection matrix: Week 7c Note 8
p1c = -reflect*p1bc;	% CoM momentum after collision
p2c = -p1c;		% CoM momentum after collision
p1 = p1c + a*pc;	% convert back from CoM coordinates
p2 = p2c + (1-a)*pc;	% ditto

% function [t,c] = molCollWallWhen(cb,p,M,r,abcd)	THM		101025
% find
%   time of collision, t
%   centre after collision:
%     [a-r;y] for right, [r-a;y] for left, [x;b-r] for top, [x,r-b] for bottom
%     - could use sum of appropriate for corner, but let's default to wall
%       if a tie, then schedule an immediate hit on top or bottom
% given
%   centre of molecule at time 0, column vector cb
%   momentum, column vector p
%   mass (amu), M
%   vector abcd describing box (abcd(1),abcd(2))*(abcd(3),abcd(4))
% uses	v = condif(cond,then,els);
%	z = implies(x,y);
%	found = nonnegrows(matrix);
%	found = minrows(matrix);
function [t,c] = molCollWallWhen(cb,p,M,r,abcd)
 m = M*10.36;		% convert from amu to milli-eV (ps/nm)^2
 v = p/m;		% velocity
% distances from location  cb  to wall or middle horizontal diaphragm
 dap = abcd(2) - r - cb(1);	% positive distance to vertical wall
 dav = abcd(1) + r - cb(1);	% negative distance to vertical wall
 dbp = abcd(4) - r - cb(2);	% positive distance to horizontal wall
 dbv = abcd(3) + r - cb(2);	% negative distance to horizontal wall
% if distance==0 then depend on sign of  v ; time is dist/v; new molecule pos, c
 record(1,:) = condif(implies(dap==0,v(1)>0),[dap/v(1),abcd(2)-r,cb(2)+v(2)*dap/v(1)],[inf,0,0]);
 record(2,:) = condif(implies(dav==0,v(1)<0),[dav/v(1),abcd(1)+r,cb(2)+v(2)*dav/v(1)],[inf,0,0]);
 record(3,:) = condif(implies(dbp==0,v(2)>0),[dbp/v(2),cb(1)+v(1)*dbp/v(2),abcd(4)-r],[inf,0,0]);
 record(4,:) = condif(implies(dbv==0,v(2)<0),[dbv/v(2),cb(1)+v(1)*dbv/v(2),abcd(3)+r],[inf,0,0]);
 next = minrows(nonnegrows(record));	% 2 identical if corner, else 1
 sizeNext = size(next);			% display test only
 t = next(1,1);			% time is smallest of the 5 times
 c = [next(1,2);next(1,3)];	% new position of molecule 4 corresponding event

% function v = condif(cond,then,els)			THM		101025
% conditional expression: getting around a MATLAB omission
function v = condif(cond,then,els)
if cond, v = then; else v = els; end % if

% function z = implies(x,y)				THM		101025
% the Boolean expression x => y ("if x then y") is formally equivalent to
%			~x|y, which is the same as y|~x ("y or not x")
function z = implies(x,y)
z = ~x | y;

% function found = nonnegrows(matrix)			THM		101025
% find rows of matrix whose first element is nonnegative
function found = nonnegrows(matrix)
sizeMatrix = size(matrix);
found = [];
for j = 1:sizeMatrix(1)
  if matrix(j,1)>=0, found = [found; matrix(j,:)]; end % if
end % for j

% function found = minrows(matrix)			THM		101025
% return all rows of a rectangular matrix with smallest first element
function found = minrows(matrix)
sizeMatrix = size(matrix);
found = matrix(1,:);
for j = 2:sizeMatrix(1)
  if matrix(j,1) < found(1,1)
    found = matrix(j,:);
  elseif matrix(j,1) == found(1,1)
    found = [found; matrix(j,:)];
  end % if compare
end % for j

% function [t,c1,c2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2) THM 101025
% find
%   time of collision, t (if imag(t)~=0 then no collision)
%   centres of each molecule at collision, column vectors c1, c2
% given
%   centres of each molecule at time 0, column vectors c1b, c2b
%   masses and momenta of each molecule, M1, M2, column vectors p1b, p2b
%   radii of each molecule, r1, r2
% uses x = quadratic(a,b,c)
function [t,c1,c2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2)
 m1 = M1*10.36;	% convert from amu to merts (meV(ps/nm)^2)
 m2 = M2*10.36;	% convert from amu to merts (meV(ps/nm)^2)
v1b = p1b/m1;		% find velocity of first molecule (nm/ps)
v2b = p2b/m2;		% find velocity of second molecule (nm/ps)
v = v1b - v2b;		% make second molecule stationary
c = c1b - c2b;		% centre second molecule at origin
r = r1 + r2;		% expand 2nd molecule to radius r, shrink 1st to point
if v'*v==0, t = inf; c1 = c1b; c2 = c2b; return, end % if	% parallel vels
t = min(quadratic(v'*v,2*c'*v,c'*c - r^2));		% display test only
c1 = c + v1b*t + c2b;	% move 1st molecule, then transform 2nd away from origin
c2 = v2b*t + c2b;	% transform 2nd away from origin
if t<-eps*10^3 | abs(imag(t))>eps	% flag coincident or no collision
  t = inf;
end % if t<0 |

% function molData = moleculeInit(N,T,startRange,r,M)		THM	101026
% given number of molecules N, temperature T, molecular radius r (nm) and
%  mass M (AMU);
%     NB 40 degC == 313.15 degK == 27 meV (26.9935 using k = .0862)
%  startRange is array (posLo, posHi, stdev, momLo, momHi) for the
%  initial positions of the molecules (angles in degrees, stdev in nm)
% find N*7 array molData, columns cx,cy,px,py,r,M,E
% uses genDirection(range)	(uniformly distributed angular direction)
%      genNormal(stdev)	(Normally distributed position on line, given stdev)
function molData = moleculeInit(N,T,startRange,r,M)
startDirRange = [startRange(1),startRange(2)]*pi/180;
stdev = startRange(3);
startLoc = [startRange(4),startRange(5)];
momDirRange = [startRange(6),startRange(7)]*pi/180;
%for k = 1:N
k = 1;
% duplicates = [];			% track duplicates: test only
while true		% loop over N molecules, but pass over any duplicates
  dirStart = genDirection(startDirRange);
  c = genNormal(stdev)*[cos(dirStart),sin(dirStart)] + startLoc;
  dirMom = genDirection(momDirRange);		     % changed 100317
  % uniform momenta (changed 100602,100606):
  E = (1 + (2*rand(1)-1)*sqrt(2)/10)*0.0862*T;	% avg +/- (-1 -- 1)rt2/10
  p = sqrt(2*10.36*M*E)*[cos(dirMom),sin(dirMom)];
  moldata = [c(1),c(2),p(1),p(2),r,M,E];
  duplicate = false;
  for j = 1:k-1			% a very slow way of checking for duplicates
    if min(moldata==molData(j,:))		% min = and for boolean array
%      duplicates = [duplicates;[k,j,moldata]];% track duplicates: test only
      duplicate = true;
      break					% full duplicate found
	% NB it doesn't matter if partial duplicate, e.g., positions coincide
    end % if min()
  end % for j
  if ~duplicate | k==1
      molData(k,:) = moldata;			% record unique initial molecule
      k = k + 1;
  end % if ~duplicate
  if k>N, break, end %if
end % while true (aka  for k)

% function a = genDirection(range)		THM			100301
% generate random direction in range(1)--range(2)
function a = genDirection(range)
lowangle = range(1);
anglerange = range(2) - lowangle;
a = lowangle + anglerange*rand(1);

% function d = genNormal(stdev)			THM			100301
% given standerd deviation stdev
% find a distance d from centre
% uses erfNormal(stdev,decpl)  to give cumulative sum of Normal
function d = genNormal(stdev)
% Following \cite[p.263]{Ross:prob} to generate random variable X from uniform u:
%       generate uniform random u on [0,1];
%	p = erfNormal^{-1}(u)
if stdev == 0, d = 0; else
  u = rand(1);			% uniform random
  for d = 0:0.03*stdev:3*stdev
    if erfNormal(d,2) > u, break, end %if
  end % for d
  s = rand(1);
  if s<=0.5, signd = -1; else signd = 1; end
  d = signd*d;
end % if stdev

% function molData = collSim(N,T,r,M,ab)	THM			101026
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return nothing (collision simulation only, no stats); molData is sop;see below
% uses	molData = moleculeInit(N,T,startRange,r,M)
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
function molData = collSim(N,T,r,M,ab)
tic					% elapsed time

% initialize
startRange = [0,0,5,0,0,40,50];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
molData = moleculeInit(N,T,startRange,r,M)	% display test only
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if mod(happIndx,stepChunk)==0
      elapsed = toc;				% elapsed time
      chunkNo = chunkNo + 1;
      if firstChunk
	firstChunk = false;
      end % if firstChunk
      beep
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/n '],'s');
      if isempty(reply)
          reply = 'y';
      end
      if reply~='y'	% not 'y'
        break
      else tic					% elapsed time
      end % if reply
    end % if mod()
  end % if collision
end % while(true)

% function molDatAverage = molDataHisto(molData,molDatAverage) THM 101028
% plot histogram of statistics from statSim()
% given
%	molData, N*7 array of [cx,cy,px,py,4,M,E] to which we'll add p
%	molDatAverage, initially [], average of all these
% find
%	molDatAverage = [molDatAverage;molDatAvg]
function molDatAverage = molDataHisto(molData,molDatAverage)
molData(:,8) = sqrt(molData(:,3).^2 + molData(:,4).^2);	% momentum magnitude
sizeMolData = size(molData);
% firstInvocation = size(molDatAverage)==[0,0];% initially [0,0] (or invokNo==1)
sizeMolDataAvg = size(molDatAverage);
invokNo = sizeMolDataAvg(1) + 1;
firstInvocation = invokNo==1;			% initially
molDatAvg = sum(molData)/sizeMolData(1)
molDatAverage = [molDatAverage;molDatAvg];
%molDataStdev = sqrt(sum(molData.^2-molDatAverage.^2)/sizeMolData(1))
molDatBucket = zeros(10,5);
bw = 10;		% bucket width
for k = 1:4				% momentum (1,2), energy (3), |p| (4)
  k1 = condif(k<3,k+2,k+4);		% 1->3 (px), 2->4 (py), 3->7 (E), 4->8
  bs = condif(k<3,-5*bw,0);		% bucket start
  for j = 1: sizeMolData(1)
    if molData(j,k1) < bs + bw, molDatBucket(1,k) = molDatBucket(1,k) + 1;
    elseif molData(j,k1) < bs + 2*bw, molDatBucket(2,k) = molDatBucket(2,k) + 1;
    elseif molData(j,k1) < bs + 3*bw, molDatBucket(3,k) = molDatBucket(3,k) + 1;
    elseif molData(j,k1) < bs + 4*bw, molDatBucket(4,k) = molDatBucket(4,k) + 1;
    elseif molData(j,k1) < bs + 5*bw, molDatBucket(5,k) = molDatBucket(5,k) + 1;
    elseif molData(j,k1) < bs + 6*bw, molDatBucket(6,k) = molDatBucket(6,k) + 1;
    elseif molData(j,k1) < bs + 7*bw, molDatBucket(7,k) = molDatBucket(7,k) + 1;
    elseif molData(j,k1) < bs + 8*bw, molDatBucket(8,k) = molDatBucket(8,k) + 1;
    elseif molData(j,k1) < bs + 9*bw, molDatBucket(9,k) = molDatBucket(9,k) + 1;
    elseif molData(j,k1) < bs + 10*bw, molDatBucket(10,k) = molDatBucket(10,k) + 1;
    else % ignore energentum over 100
    end % if molData
  end % for j
end % for k
colour = 'rbgkm';
invokNo,clr = colour(mod(invokNo-1,5)+1)	% display test only
j=1:10;
subplot(3,2,1)
if firstInvocation, hold off, end
plot(bw*(j-5),molDatBucket(:,1),clr)
title('momentum x')
hold on
subplot(3,2,2)
if firstInvocation, hold off, end
plot(bw*(j-5),molDatBucket(:,2),clr)
title('momentum y')
hold on
subplot(3,2,3)
if firstInvocation, hold off, end
plot(bw*j,molDatBucket(:,4),clr)
title('momentum magn.')
hold on
subplot(3,2,4)
if firstInvocation, hold off, end
plot(bw*j,molDatBucket(:,3),clr)
title('energy')
hold on
subplot(3,2,5)
if firstInvocation, hold off, end
scatter(molDatAvg(3),molDatAvg(4),clr)
xlabel('p_x')
ylabel('p_y')
title('average momentum')
hold on
subplot(3,2,6)
if firstInvocation, hold off, end
scatter(molDatAvg(1),molDatAvg(2),clr)
xlabel('c_x')
ylabel('c_y')
title('average position')
hold on

% function molDatAverage = molDataHisto(molData,molDatAverage) THM 101028
% plot histogram of statistics from statSim()
% given
%	molData, N*7 array of [cx,cy,px,py,4,M,E] to which we'll add p
%	molDatAverage, initially [], average of all these
% find
%	molDatAverage = [molDatAverage;molDatAvg]
function molDatAverage = molDataHisto(molData,molDatAverage)
molData(:,8) = sqrt(molData(:,3).^2 + molData(:,4).^2);	% momentum magnitude
sizeMolData = size(molData);
% firstInvocation = size(molDatAverage)==[0,0];% initially [0,0] (or invokNo==1)
sizeMolDataAvg = size(molDatAverage);
invokNo = sizeMolDataAvg(1) + 1;
firstInvocation = invokNo==1;			% initially
molDatAvg = sum(molData)/sizeMolData(1)
molDatAverage = [molDatAverage;molDatAvg];
%molDataStdev = sqrt(sum(molData.^2-molDatAverage.^2)/sizeMolData(1))
molDatBucket = zeros(10,5);
bw = 10;		% bucket width
for k = 1:4				% momentum (1,2), energy (3), |p| (4)
  k1 = condif(k<3,k+2,k+4);		% 1->3 (px), 2->4 (py), 3->7 (E), 4->8
  bs = condif(k<3,-5*bw,0);		% bucket start
  for j = 1: sizeMolData(1)
    if molData(j,k1) < bs + bw, molDatBucket(1,k) = molDatBucket(1,k) + 1;
    elseif molData(j,k1) < bs + 2*bw, molDatBucket(2,k) = molDatBucket(2,k) + 1;
    elseif molData(j,k1) < bs + 3*bw, molDatBucket(3,k) = molDatBucket(3,k) + 1;
    elseif molData(j,k1) < bs + 4*bw, molDatBucket(4,k) = molDatBucket(4,k) + 1;
    elseif molData(j,k1) < bs + 5*bw, molDatBucket(5,k) = molDatBucket(5,k) + 1;
    elseif molData(j,k1) < bs + 6*bw, molDatBucket(6,k) = molDatBucket(6,k) + 1;
    elseif molData(j,k1) < bs + 7*bw, molDatBucket(7,k) = molDatBucket(7,k) + 1;
    elseif molData(j,k1) < bs + 8*bw, molDatBucket(8,k) = molDatBucket(8,k) + 1;
    elseif molData(j,k1) < bs + 9*bw, molDatBucket(9,k) = molDatBucket(9,k) + 1;
    elseif molData(j,k1) < bs + 10*bw, molDatBucket(10,k) = molDatBucket(10,k) + 1;
    else % ignore energentum over 100
    end % if molData
  end % for j
end % for k
colour = 'rbgkm';
invokNo,clr = colour(mod(invokNo-1,5)+1)	% display test only
j=1:10;
subplot(3,2,1)
if firstInvocation, hold off, end
plot(bw*(j-5),molDatBucket(:,1),clr)
title('momentum x')
hold on
subplot(3,2,2)
if firstInvocation, hold off, end
plot(bw*(j-5),molDatBucket(:,2),clr)
title('momentum y')
hold on
subplot(3,2,3)
if firstInvocation, hold off, end
plot(bw*j,molDatBucket(:,4),clr)
title('momentum magn.')
hold on
subplot(3,2,4)
if firstInvocation, hold off, end
plot(bw*j,molDatBucket(:,3),clr)
title('energy')
hold on
subplot(3,2,5)
if firstInvocation, hold off, end
scatter(molDatAvg(3),molDatAvg(4),clr)
xlabel('p_x')
ylabel('p_y')
title('average momentum')
hold on
subplot(3,2,6)
if firstInvocation, hold off, end
scatter(molDatAvg(1),molDatAvg(2),clr)
xlabel('c_x')
ylabel('c_y')
title('average position')
hold on

% function molData = statSim(N,T,r,M,ab)	THM			101026
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = moleculeInit(N,T,startRange,r,M)
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
%	molDatAverage = molDataHisto(molData,[]);
function molData = statSim(N,T,r,M,ab)

% initialize
startRange = [0,0,5,0,0,40,50];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
molData = moleculeInit(N,T,startRange,r,M)	% display test only
molDatAverage = molDataHisto(molData,[]);
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if mod(happIndx,stepChunk)==0
      elapsed = toc;				% elapsed time
      chunkNo = chunkNo + 1;
      if firstChunk
	firstChunk = false;
      end % if firstChunk
      molDatAverage = molDataHisto(molData,molDatAverage);
      beep
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/n '],'s');
      if isempty(reply)
          reply = 'y';
      end
      if reply~='y'	% not 'y'
        break
      else tic					% elapsed time
      end % if reply
    end % if mod()
  end % if collision
end % while(true)

% function [momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg) THM 101103
% accumulated magnitude of momentum and energy from molData into momenCum
% given
%	molData, N*7 array of [cx,cy,px,py,r,M,E] to which we'll add p
%	momenCum, initially zeros(11,2): 10 positions each for p, E
%	momenAvg, initially []: accumulates average momentum magnitude, energy
% find
%	momenCum
%	momenAvg
function [momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg)
molData(:,8) = sqrt(molData(:,3).^2 + molData(:,4).^2); % momentum magnitude
sizeMolData = size(molData);
molDatBucket = zeros(10,5);
bs = 0;				% bucket start
bw = 10;			% bucket width
for k = 1:2			% momentum magnitude (1), energy (2)
  k1 = condif(k==1,8,7);	% molData(8) is p, (7) is E
  for j = 1: sizeMolData(1)
    if molData(j,k1) < bs + bw, momenCum(1,k) = momenCum(1,k) + 1;
    elseif molData(j,k1) < bs + 2*bw, momenCum(2,k) = momenCum(2,k) + 1;
    elseif molData(j,k1) < bs + 3*bw, momenCum(3,k) = momenCum(3,k) + 1;
    elseif molData(j,k1) < bs + 4*bw, momenCum(4,k) = momenCum(4,k) + 1;
    elseif molData(j,k1) < bs + 5*bw, momenCum(5,k) = momenCum(5,k) + 1;
    elseif molData(j,k1) < bs + 6*bw, momenCum(6,k) = momenCum(6,k) + 1;
    elseif molData(j,k1) < bs + 7*bw, momenCum(7,k) = momenCum(7,k) + 1;
    elseif molData(j,k1) < bs + 8*bw, momenCum(8,k) = momenCum(8,k) + 1;
    elseif molData(j,k1) < bs + 9*bw, momenCum(9,k) = momenCum(9,k) + 1;
    elseif molData(j,k1) < bs + 10*bw, momenCum(10,k) = momenCum(10,k) + 1;
    else % ignore energentum over 100
    end % if molData
  end % for j
end % for k
momenAverage(1) = sum(molData(:,8))/10;		% average momentum magnitude
momenAverage(2) = sum(molData(:,7))/10;		% average energy
momenAvg = [momenAvg;momenAverage];

% function molDatAverage =
%	 molDataDistr(molData,molDatAverage,momenCum,momenAvg)	THM	101028
% plot histogram of statistics from statSim()
% given
%	molData, N*7 array of [cx,cy,px,py,r,M,E] to which we'll add p
%	molDatAverage, initially [], average of all these
%	momenCum, array of 10 buckets for each of momentum magnitude, energy
%	momenAvg, array of variable rows holding average mom magn, energy
%		(both these last supplied by function momenCumHisto)
% find
%	molDatAverage = [molDatAverage;molDatAvg]
function molDatAverage = molDataDistr(molData,molDatAverage,momenCum,momenAvg)
molData(:,8) = sqrt(molData(:,3).^2 + molData(:,4).^2);	% momentum magnitude
sizeMolData = size(molData);
% firstInvocation = size(molDatAverage)==[0,0];% initially [0,0] (or invokNo==1)
sizeMolDataAvg = size(molDatAverage);
invokNo = sizeMolDataAvg(1) + 1;
firstInvocation = invokNo==1;			% initially
secndInvocation = invokNo==2;	% or size(momenAvg)==[0,0]
momenAvg			% display test only
molDatAvg = sum(molData)/sizeMolData(1)		% display test only
molDatAverage = [molDatAverage;molDatAvg];
%molDataStdev = sqrt(sum(molData.^2-molDatAverage.^2)/sizeMolData(1))
molDatBucket = zeros(10,5);
bw = 10;		% bucket width
for k = 1:4	% kept, tho need only 2	% momentum (1,2), energy (3), |p| (4)
  k1 = condif(k<3,k+2,k+4);		% 1->3 (px), 2->4 (py), 3->7 (E), 4->8
  bs = condif(k<3,-5*bw,0);		% bucket start
  for j = 1: sizeMolData(1)
    if molData(j,k1) < bs + bw, molDatBucket(1,k) = molDatBucket(1,k) + 1;
    elseif molData(j,k1) < bs + 2*bw, molDatBucket(2,k) = molDatBucket(2,k) + 1;
    elseif molData(j,k1) < bs + 3*bw, molDatBucket(3,k) = molDatBucket(3,k) + 1;
    elseif molData(j,k1) < bs + 4*bw, molDatBucket(4,k) = molDatBucket(4,k) + 1;
    elseif molData(j,k1) < bs + 5*bw, molDatBucket(5,k) = molDatBucket(5,k) + 1;
    elseif molData(j,k1) < bs + 6*bw, molDatBucket(6,k) = molDatBucket(6,k) + 1;
    elseif molData(j,k1) < bs + 7*bw, molDatBucket(7,k) = molDatBucket(7,k) + 1;
    elseif molData(j,k1) < bs + 8*bw, molDatBucket(8,k) = molDatBucket(8,k) + 1;
    elseif molData(j,k1) < bs + 9*bw, molDatBucket(9,k) = molDatBucket(9,k) + 1;
    elseif molData(j,k1) < bs + 10*bw, molDatBucket(10,k) = molDatBucket(10,k) + 1;
    else % ignore energentum over 100
    end % if molData
  end % for j
end % for k
for j = 1:2					% normalize histo to distr
  momenCumDistr(:,j) = momenCum(:,j)/sum(momenCum(:,j));
end % for j
colour = 'bgkm';
invokNo,clr = colour(mod(invokNo-1,4)+1)	% display test only
j=1:10;
subplot(3,2,1)
if firstInvocation, hold off, end
plot(bw*j,molDatBucket(:,4),clr)
title('momentum magn.')
hold on
subplot(3,2,2)
if firstInvocation, hold off, end
plot(bw*j,molDatBucket(:,3),clr)
title('energy')
hold on
subplot(3,2,3)
if firstInvocation, hold off, end
plot(bw*j,momenCumDistr(:,1),clr)
title('cumulative momentum magn.')
hold on
subplot(3,2,4)
if firstInvocation, hold off, end
plot(bw*j,momenCumDistr(:,2),clr)
title('cumulative energy')
hold on
subplot(3,2,5)
if firstInvocation, hold off, end
plot(momenAvg(:,1))
title('average momentum magn.')
xlabel('time (# collisions)')
ylabel('meV-ps/nm')
hold on
subplot(3,2,6)
if firstInvocation, hold off, end
plot(momenAvg(:,2))
title('average energy')
xlabel('time (# collisions)')
ylabel('meV')
hold on

% function molData = cumuSim(N,T,r,M,ab)	THM			101102
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = moleculeInit(N,T,startRange,r,M)
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
%	[momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg)
%	molDatAverage = molDataDistr(molData,molDatAverage,momenCum,momenAvg);
% Differs from statSim by introducing cumulative momentum magnitude and energy;
%  modifies histogram building and plotting routine, replacing momx, momy by
%  these new distributions, and scatter plots of mom, pos by fluctuation plots
%  of average momentum magnitude and energy.
function molData = cumuSim(N,T,r,M,ab)

% initialize
startRange = [0,0,5,0,0,40,50];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
molData = moleculeInit(N,T,startRange,r,M)	% display test only
molDatAverage = [];	% initialize for molDataDistr: no longer display initial
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
momenCum = zeros(10,2);	% initialize parameters for momenCumHist
momenAvg = [];		%  ditto
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if ~firstChunk
    % start accumulating, skipping the distortion of the initial motions
      if mod(happIndx,N)==0		% accumulate for every N collisions
	[momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg);
      end % if mod(happIndx,N)==0
    end % if ~firstChunk
    if mod(happIndx,stepChunk)==0
      elapsed = toc;				% elapsed time
      chunkNo = chunkNo + 1;
      if firstChunk
	firstChunk = false;
      else
	molDatAverage = molDataDistr(molData,molDatAverage,momenCum,momenAvg);
      end % if firstChunk
      beep
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/n '],'s');
      if isempty(reply)
          reply = 'y';
      end
      if reply~='y'	% not 'y'
        break
      else tic					% elapsed time
      end % if reply
    end % if mod()
  end % if collision
end % while(true)

% function [momAvgFT,mAapprox] = FTfluct		THM		101109
% Find discrete Fourier transform of hardwu=ired sequence momAvg: momAvgFT
% Manually select strongest frequency components (7 with largest magnitudes)
% and DFT back to mAapprox
function [momAvgFT,mAapprox] = FTfluct
momAvg = [42.2964;41.1153;42.0875;41.5578;42.5435;41.9132;40.2426;42.4980;41.3177;42.8827;41.9394;40.3250;40.0651;37.5326;37.4996;39.1918;40.4985;40.1718;39.0108;39.7756;39.3788;37.1139;39.9479;40.1349;38.4678;36.7253;34.7730;40.8269;39.9655;39.5493;40.2808;39.6370;38.1273;41.9592;41.8759;44.2589;42.9165;43.0852;42.8527;42.0199];
F40 = makeDFT(40);
momAvgFT = F40*momAvg;
mAapproxFT = zeros(size(momAvgFT));
momAvgFT = [momAvgFT, abs(momAvgFT)];
mAapproxFT(1,1) = momAvgFT(1,1);
mAapproxFT(2,1) = momAvgFT(2,1);
mAapproxFT(7,1) = momAvgFT(7,1);
mAapproxFT(4,1) = momAvgFT(4,1);
mAapproxFT(8,1) = momAvgFT(8,1);
mAapproxFT(17,1) = momAvgFT(17,1);
mAapproxFT(5,1) = momAvgFT(5,1);
for j = 2:19	% put in symmetric components
  mAapproxFT(42-j,1) = momAvgFT(j,1);
end % for j

mAapprox = conj(F40')*mAapproxFT;
mAapprox = [mAapprox,abs(mAapprox)];
plot(1:40,momAvg,'r',1:40,mAapprox(:,2),'b',1:40,momAvgFT(1,1)/sqrt(40),'k')
title('Fourier analysis of average momentum magnitude')
xlabel('time (# collisions/10)')
ylabel('average momentum magnitude')
legend('avg mom','largest 7 freq.','time average','Location','North')

% function momAvgFT = FTfluctAgg(aggSize,clr)		THM	101111
% Group the 40 values hardwired in momAvg into 40/aggSize (must be integer)
% aggregate averages. So aggSize may be 1,2,4,5,8,10,20,40. Then, with these,
% Take Fourier transform to find the principal frequencies.
%	clr is the plot colour
function momAvgFT = FTfluctAgg(aggSize,clr)
momAvg = [42.2964;41.1153;42.0875;41.5578;42.5435;41.9132;40.2426;42.4980;41.3177;42.8827;41.9394;40.3250;40.0651;37.5326;37.4996;39.1918;40.4985;40.1718;39.0108;39.7756;39.3788;37.1139;39.9479;40.1349;38.4678;36.7253;34.7730;40.8269;39.9655;39.5493;40.2808;39.6370;38.1273;41.9592;41.8759;44.2589;42.9165;43.0852;42.8527;42.0199];
sizeMomAvg = size(momAvg);
k = 0;
for j = 1:sizeMomAvg
  if aggSize==1 | mod(j,aggSize)==1
    if k>0, cumMomAvg(k) = cumMomAvg(k)/aggSize; end
    k = k + 1;
    cumMomAvg(k) = 0;
  end % if mod
  cumMomAvg(k) = cumMomAvg(k) + momAvg(j);
end % for k
cumMomAvg(k) = cumMomAvg(k)/aggSize;	% last term
noTerms = sizeMomAvg/aggSize;		% must be integer
DFT = makeDFT(noTerms(1));
momAvgFT = DFT*cumMomAvg';
momAvgFT = [momAvgFT, abs(momAvgFT)];
figure(3)			% hardwired: will wipe previous figure 3
plot(aggSize:aggSize:40,cumMomAvg,clr)
hold on				% must hold off before first invocation
title('momentum magnitude averages fluctuate')
xlabel('time (# collisions/10)')
ylabel('average momentum magnitude')

% function mimStats = momAvgHisto(aggSize,clr)		THM		101111
% Find and plot 10-element histogram for hardwired momentum magnitude averages:
% group the 40 values hardwired in momAvg into 40/aggSize (must be integer)
% aggregate averages. So aggSize may be 1,2,4,5,8,10,20,40.
%	clr is colour of plot
function momStats = momAvgHisto(aggSize,clr)
% this one is really the sequence of average momenta
% plotFig = 4;
% momAvg = [42.2964;41.1153;42.0875;41.5578;42.5435;41.9132;40.2426;42.4980;41.3177;42.8827;41.9394;40.3250;40.0651;37.5326;37.4996;39.1918;40.4985;40.1718;39.0108;39.7756;39.3788;37.1139;39.9479;40.1349;38.4678;36.7253;34.7730;40.8269;39.9655;39.5493;40.2808;39.6370;38.1273;41.9592;41.8759;44.2589;42.9165;43.0852;42.8527;42.0199];
% this one is really a set of average energies (from ../newsim1/manyInits)
plotFig = 5;
momAvg = [25.4606;26.7735;26.6858;26.1923;26.5929;26.2111;25.7814;25.8228;25.9716;26.2799;25.4807;25.9138;25.8934;26.5559;26.9727;25.4016;26.6654;26.2584;26.4078;26.2142;25.2606;26.1378;25.2521;26.5894;25.4875;26.5820;25.3646;25.6381;25.3599;26.5136;26.7209;24.2845;25.7262;25.5478;25.3614; 24.1373;26.6495;26.1698;25.7589;25.8821];
sizeMomAvg = size(momAvg);
k = 0;
for j = 1:sizeMomAvg
  if aggSize==1 | mod(j,aggSize)==1
    if k>0, cumMomAvg(k) = cumMomAvg(k)/aggSize; end
    k = k + 1;
    cumMomAvg(k) = 0;
  end % if mod
  cumMomAvg(k) = cumMomAvg(k) + momAvg(j);
end % for k
cumMomAvg(k) = cumMomAvg(k)/aggSize;	% last term
noTerms = sizeMomAvg/aggSize;		% must be integer
sizeCumMomAvg = size(cumMomAvg);
frstMom = 0;
scndMom = 0;
for j = 1:noTerms
  frstMom = frstMom + cumMomAvg(j);
  scndMom = scndMom + cumMomAvg(j)^2;
end % for j
momStats(1) = frstMom/noTerms(1);			% mean
momStats(2) = scndMom/noTerms(1) - momStats(1)^2;	% variance
momStats(3) = sqrt(momStats(2));		% standard deviation
minCumMomAvg = min(cumMomAvg);
maxCumMomAvg = max(cumMomAvg);
noBuck = condif(aggSize==1,10,condif(aggSize==2,5,2));
stepBuck = (maxCumMomAvg - minCumMomAvg)/(noBuck - 1);
histo = zeros(1,noBuck);
for j = 1:sizeCumMomAvg(2);
  if cumMomAvg(j) < minCumMomAvg + stepBuck, histo(1) = histo(1) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 2*stepBuck, histo(2) = histo(2) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 3*stepBuck, histo(3) = histo(3) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 4*stepBuck, histo(4) = histo(4) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 5*stepBuck, histo(5) = histo(5) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 6*stepBuck, histo(6) = histo(6) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 7*stepBuck, histo(7) = histo(7) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 8*stepBuck, histo(8) = histo(8) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 9*stepBuck, histo(9) = histo(9) + 1;
  elseif cumMomAvg(j) < minCumMomAvg + 10*stepBuck, histo(10) = histo(10) + 1;
  else % ignore anything over maxMomAvg
  end % if
end % for j
x = minCumMomAvg:stepBuck:maxCumMomAvg;
figure(plotFig)			% hardwired: will erase previous figure 4 or 5
bar(x,histo,clr)
hold on				% must preface 1st invocation by hold off
% title('Histogram of momentum magnitude averages')
% xlabel('momentum, meV-ps/ns')
title('Histogram of energy averages')
xlabel('energy, meV')
ylabel('# occurrences')

% function variability = FTfluctVar		THM			101111
% Hardwire the Fourier transform magnitudes from FTflcutAgg and add up total
% magnitude as vector length. Each set gives an element of  variability  output
function variability = FTfluctVar
momVar39=[0.0642;0.0148;0.0228;0.0171;0.0085;0.0282;0.0201;0.0106;0.0047;0.0152;0.0157;0.0069;0.0045;0.0111;0.0118;0.0174;0.0061;0.0114;0.0079;0.0098;0.0079;0.0114;0.0061;0.0174;0.0118;0.0111;0.0045;0.0069;0.0157;0.0152;0.0047;0.0106;0.0201;0.0282;0.0085;0.0171;0.0228;0.0148;0.0642]*100/sqrt(40);
momVar19=[0.0448;0.0095;0.0154;0.0081;0.0024;0.0168;0.0113;0.0033;0.0057;0.0035;0.0057;0.0033;0.0113;0.0168;0.0024;0.0081;0.0154;0.0095;0.0448]*100/sqrt(20);
momVar9=[0.0307;0.0068;0.0122;0.0098;0.0013;0.0098;0.0122;0.0068; 0.0307]*100/sqrt(10);;
momVar4=[1.8550;0.8812;0.8812;1.8550]/sqrt(5);
variability(1) = vecLength(momVar39);
variability(2) = vecLength(momVar19);
variability(3) = vecLength(momVar9);
variability(4) = vecLength(momVar4);

% function avgEnergy = manyInits		THM			101111
% call moleculeInit 40 times, 10 molecules eqch, and record average energies
function avgEnergy = manyInits
for j = 1:40
  molData = moleculeInit(10,300,[0,0,5,0,0,40,50],0.05,4);
  sizeMolData = size(molData);
  totEn = 0;
  for k = 1:sizeMolData(1)
    totEn = totEn + molData(k,7);
  end % for k
  avgEnergy(j) = totEn/sizeMolData(1);
end % for j

% function enerIgn = energyHisDis(momenCum,enerIgn)		THM	101116
% finds and plots energy ignorance from momenCum(2)
% in and out: enerIgn, column vector of successive values of energy ignorance
function enerIgn = energyHisDis(momenCum,enerIgn)
enerProb = momenCum(:,2)/sum(momenCum(:,2));	% display test only
sizeEnerProb = size(enerProb);		% display test only
enerIg = 0;
for j = 1: sizeEnerProb(1)
  enerP = enerProb(j);
  enerIg = enerIg - enerP.*condif(enerP==0,1,log(enerP));
end % for j
enerIgn = [enerIgn;enerIg];
sizeEnerIgn = size(enerIgn);		% display test only
figure(6)	% hardwired figure 6 will overwrite any previous figure 6
plot(0:sizeEnerIgn(1)-1,enerIgn)
title('Energy ignorance vs time')
xlabel('time (# collisions)')
ylabel('energy ignorance (nats)')

% function molData = entrSim(N,T,r,M,ab)	THM			101102
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = molecuReInit(N,T,startRange,r,M)	fixed initial config.
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
%	[momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg)
%	enerIgn = energyHisDis(momenCum,enerIgn);
% Differs from cumuSim by finding and plotting cumulative energy ignorance
% (entropy).
function molData = entrSim(N,T,r,M,ab)

% initialize
startRange = [0,0,5,0,0,40,50];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
% use one or other of the following two lines:
%molData = molecuReInit(N,T,startRange,r,M)	% display test only
molData = moleculeInit(N,T,startRange,r,M)	% display test only
momenCum = zeros(10,2);	% initialize parameters for momenCumHist
momenAvg = [];		%  ditto
[momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg);
enerIgn = energyHisDis(momenCum,[]); % capture and plot initial energy ignorance
momenCum = zeros(10,2);	% reset parameters for momenCumHist
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if mod(happIndx,1)==0		% accumulate for every N collisions
      [momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg);
      enerIgn = energyHisDis(momenCum,enerIgn);
      momenCum = zeros(10,2);		% reset so accumulates over only 10 coll
    end % if mod(happIndx,N)==0
    if mod(happIndx,stepChunk)==0
      j=1:1000;
      beta = 1/25.0221;			% reciprocal of mean energy
      Zterms = exp(-beta*10*j);		% bucket width bw = 10
      Z = sum(Zterms);
      BoltzIgn = 1 + log(Z)		% display
      hold on
      plot(1:happIndx,1.92*ones(happIndx,1),'r')	% use theor max
      hold off
      enerIgn'				% display
      elapsed = toc;			% elapsed time
      chunkNo = chunkNo + 1;
      beep
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/n '],'s');
      if isempty(reply)
          reply = 'y';
      end
      if reply~='y'	% not 'y'
        break
      else tic					% elapsed time
      end % if reply
    end % if mod()
  end % if collision
end % while(true)

% function [t,c] = molCollWallWhen(cb,p,M,r,abcd)	THM		101117
% find
%   time of collision, t
%   centre after collision:
%     [a-r;y] for right, [r-a;y] for left, [x;b-r] for top, [x,r-b] for bottom
%     - could use sum of appropriate for corner, but let's default to wall
%       if a tie, then schedule an immediate hit on top or bottom
% given
%   centre of molecule at time 0, column vector cb
%   momentum, column vector p
%   mass (amu), M
%   vector abcd describing box (abcd(1),abcd(2))*(abcd(3),abcd(4))
% uses	v = condif(cond,then,els);
%	z = implies(x,y);
%	found = nonnegrows(matrix);
%	found = minrows(matrix);
% Differs from corresponding routine in newsim1 by testing for diaphragm at y=0.
function [t,c] = molCollWallWhen(cb,p,M,r,abcd)
 m = M*10.36;		% convert from amu to milli-eV (ps/nm)^2
 v = p/m;		% velocity
% distances from location  cb  to wall or middle horizontal diaphragm
 dap = abcd(2) - r - cb(1);	% positive distance to vertical wall
 dav = abcd(1) + r - cb(1);	% negative distance to vertical wall
 dbp = abcd(4) - r - cb(2);	% positive distance to horizontal wall
 dbv = abcd(3) + r - cb(2);	% negative distance to horizontal wall
 db0 = -cb(2);		% NB -ve if above y=0 so needs -ve v(2) to cross
% if distance==0 then depend on sign of  v ; time is dist/v; new molecule pos, c
 record(1,:) = condif(implies(dap==0,v(1)>0),[dap/v(1),abcd(2)-r,cb(2)+v(2)*dap/v(1)],[inf,0,0]);
 record(2,:) = condif(implies(dav==0,v(1)<0),[dav/v(1),abcd(1)+r,cb(2)+v(2)*dav/v(1)],[inf,0,0]);
 record(3,:) = condif(implies(dbp==0,v(2)>0),[dbp/v(2),cb(1)+v(1)*dbp/v(2),abcd(4)-r],[inf,0,0]);
 record(4,:) = condif(implies(dbv==0,v(2)<0),[dbv/v(2),cb(1)+v(1)*dbv/v(2),abcd(3)+r],[inf,0,0]);
 record(5,:) = condif(db0~=0,[db0/v(2),cb(1)+v(1)*db0/v(2),0],[inf,0,0]);
        % if at diapraghm (db0==0) then it has already been processed
 next = minrows(nonnegrows(record));	% 2 identical if corner, else 1
 sizeNext = size(next);			% display test only
 t = next(1,1);			% time is smallest of the 5 times
 c = [next(1,2);next(1,3)];	% new position of molecule 4 corresponding event

% function [p,cntB] = molCollWallCountBot(pb,tmult,bmult,abcd,r,c,cntB,mix)
%								THM	101117
% (initial version, extracted from final)
% given vectors momentum pb, box boundaries abcd, position c; scalars radius r,
%  integer count cntB, boolean mix (0 if y=0 diaphragm closed and gases mix;
%  1 if diaphragm open and gases don't mix.
%  The multipliers tmult, bmult change top, bottom momenta due to moving wall,
%	heat, respectively. (Set to 1 for no heat transfer, fixed volume.)
% find new momentum p
% Differs from same in newsim1 by interaction with y=0 diaphragm, by counting
%  number of collisions with bottom wall, and by allowing momentum multipliers
%  at top and bottom.
function p = molCollWallCountBot(pb,tmult,bmult,abcd,r,c,cntB,mix)
% tmult = 1; bmult = 1; cntB = 0;		% for future refinement
abr = abcd - [-r,r,-r,r];	% extend wall to _centre_ of molecule
p = pb;				% initialize: momentum out = momentum in
if c(1) < abr(1) | c(1) > abr(2) | c(2) < abr(3) | c(2) > abr(4)
  display(['molCollWallCountBot: ' num2str(c') ' out of box ' num2str(abr)])
elseif (c(1)==abr(1) | c(1)==abr(2)) & (c(2)==abr(3) | c(2)==abr(4)), p = -pb;
					% 4 corners: reverse momentum
elseif (c(1)==abr(1) | c(1)==abr(2)) & c(2) == 0, p = -pb; % diaphragm ends:   "
elseif c(1)==abr(1) | c(1)==abr(2), p(1) = -pb(1);%side wall: reverse x-momentum
elseif c(2) == abr(4), p(2) = -tmult^1.5*pb(2);%top: (approx) reverse y-momentum
			p(1) = tmult^.5*pb(1);
elseif c(2) == abr(3), p(2) = -bmult*pb(2);	%bot: default reverse y-momentum
			p(1) = bmult*pb(1);
			cntB = cntB + 1;
elseif c(2) == 0
  if mix, p = pb;		% diaphragm open: default keep momentum
  else p(2) = -pb(2);		% diaphragm closed: reverse y-momentum
  end % if mix
else display('molCollWallCountBot: no wall collision')
end %if

% 101117
% >> p = molCollWallCountBot([0;1],1,1,[-20,20,-10,10],.05,[0;0],0,1) 
% p =  0
%      1 
% >> p = molCollWallCountBot([0;1],1,1,[-20,20,-10,10],.05,[0;0],0,0) 
% p =  0
%     -1

% error on l. 127? : if executes only when happIndx==1	110223
% function molData = tmptureSim(N,T,r,M,ab)	THM			101118
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% N,T,r,M 2*1 column vectors to accommodate 2 gases
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = molecuReInit(N,T,startRange,r,M)	fixed initial config.
%    or	molData1 = moleculeInit(N(1),T(1),startRange1,r(1),M(1))
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0,mix);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
%	[momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg)
%	enerIgn = energyHisDis(momenCum,enerIgn);
% Differs from entrSim by allowing 2 gases in upper and lower half of box,
%  and possibility of removing diaphragm separating them (say after initial
%  equilibrium reached); uses entrSim's entropy routine and. like it plots
%  entropy in steps of 1. Instead of using cumuSim's molDataDistr(), the
%  momentum magnitude histograms are plotted directly from the main program.
% (A subsequent comparison is done of the results with theory by
%	[pAvg,pStdev] = plotMaxwell(T,M) )
function molData = tmptureSim(N,T,r,M,ab)

% initialize
startRange1 = [0,0,5,0,0,40,50];
startRange2 = [0,0,5,0,0,-50,-40];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
% use one or other of the following two lines:
%molData = molecuReInit(N,T,startRange,r,M)	% display test only
molData1 = moleculeInit(N(1),T(1),startRange1,r(1),M(1))   % display test only
molData2 = moleculeInit(N(2),T(2),startRange2,r(2),M(2))   % display test only
molData = [molData1;molData2];
momenCum = zeros(10,2);	% initialize parameters for momenCumHist
momenAvg = [];		%  ditto
[momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg);
enerIgn = energyHisDis(momenCum,[]); % capture and plot initial energy ignorance
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
mix = 0;		% initially don't mix the gases in top and bottom
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0,mix);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
%   if mod(happIndx,max(N))==0		% accumulate for every N collisions
    if mod(happIndx,1==0)		% snapshot every collisions
      [momenCum,momenAvg] = momenCumHisto(molData,momenCum,momenAvg);
      enerIgn = energyHisDis(momenCum,enerIgn);
      momenCumDistr(:,1) = momenCum(:,1)/sum(momenCum(:,1));	% momMagn distr
      figure(8)		% hardwired: will overwrite any existing figure 8
      j=1:10;
      bw = 10;
      plot(bw*j,momenCumDistr(:,1))
      title('cumulative momentum magn.')
      xlabel('momentum, meV-ps/nm')
      ylabel('probability')
      momenCum = zeros(10,2);		% reset so accumulates over only 10 coll
    end % if mod(happIndx,N)==0
    if mod(happIndx,stepChunk)==0
      enerIgn'				% display
      elapsed = toc;			% elapsed time
      chunkNo = chunkNo + 1;
      beep
      avgEnergy1 = sum(molData(1:N(1),7))/N(1);
      avgEnergy2 = sum(molData(N(1)+1:N(1)+N(2),7))/N(2);
      energyX = (molData(:,3).^2./(2*10.36*molData(:,6)));
      avgEnergyX = sum(energyX)/(N(1)+N(2));	% check equipartition
      energyY = (molData(:,4).^2./(2*10.36*molData(:,6)));
      avgEnergyY = sum(energyY)/(N(1)+N(2));
      disp(['Average energies: gas 1 ' num2str(avgEnergy1) ' gas 2 ' ...
	num2str(avgEnergy2) ' x partition ' num2str(avgEnergyX) ...
	' y partition ' num2str(avgEnergyY)])
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/m[mix]/n '],'s');
      if isempty(reply)
	reply = 'y';
      end % if isempty
      if reply=='m'			% toggle mixing gases
	mix = ~mix;
	reply = 'y';
      end % if reply==
      if reply~='y'	% not 'y'
        break
      else tic					% elapsed time
      end % if reply~=
    end % if mod()
  end % if collision
end % while(true)

% function smooth2enerIgn(runlen)		THM			101118
% ruunig average smoothing of energy ignorance found by ../newsim2/tmptureSim()
% input  runlen  length of smoothing runs
function smooth2enerIgn(runlen)
enerIgn1 =...	% 72 lines
[0.9797; 0.9797; 0.9896; 0.9797; 0.9687; 0.9687; 0.9687;...
0.9797; 1.0287; 1.0287; 1.0287; 1.0287; 1.0287; 1.0177;...
1.0177; 1.0177; 1.0056; 1.0056; 1.0056; 1.0056; 1.0056;...
1.0056; 1.0056; 1.0056; 1.0056; 1.0056; 1.0056; 1.0056;...
1.0056; 1.0056; 1.0056; 1.0056; 1.0056; 1.0056; 1.0056;...
1.0056; 1.0056; 1.0056; 1.0056; 1.0406; 1.0406; 1.0406;...
1.0406; 1.0406; 1.0406; 1.0406; 1.0406; 1.0406; 1.0406;...
1.0406; 1.0406; 1.0406; 1.0406; 1.0406; 1.0406; 1.0406;...
1.0406; 1.0406; 1.0406; 1.0273; 1.0273; 1.0273; 1.0273;...
1.0273; 1.0273; 1.0273; 1.0871; 1.0871; 1.1602; 1.2374;...
1.2929; 1.2929; 1.3236; 1.3700; 1.3700; 1.3647; 1.3507;...
1.3740; 1.3902; 1.3816; 1.4147; 1.3979; 1.4034; 1.3694;...
1.3577; 1.3330; 1.3330; 1.3635; 1.3724; 1.3753; 1.3753;...
1.3648; 1.3878; 1.3799; 1.3925; 1.3920; 1.3611; 1.3733;...
1.3560; 1.3918; 1.4063; 1.4231; 1.4286; 1.4477; 1.4477;...
1.4697; 1.4925; 1.4832; 1.4936; 1.4680; 1.4680; 1.4642;...
1.4642; 1.4642; 1.4733; 1.4783; 1.4822; 1.4609; 1.4609;...
1.4609; 1.4426; 1.4426; 1.4347; 1.4347; 1.4347; 1.4563;...
1.4563; 1.4531; 1.4498; 1.4267; 1.4699; 1.4926; 1.4789;...
1.4873; 1.4873; 1.4873; 1.4699; 1.4699; 1.4827; 1.4845;...
1.4827; 1.4842; 1.4640; 1.4535; 1.4523; 1.4523; 1.4275;...
1.4462; 1.4462; 1.4462; 1.4456; 1.4657; 1.4614; 1.4614;...
1.4614; 1.4718; 1.4718; 1.4963; 1.4963; 1.4963; 1.4991;...
1.5018; 1.4965; 1.4937; 1.4937; 1.4965; 1.4408; 1.4408;...
1.4408; 1.4489; 1.4330; 1.4315; 1.4315; 1.4366; 1.4344;...
1.4309; 1.4363; 1.4152; 1.4179; 1.4179; 1.4179; 1.4179;...
1.3943; 1.3943; 1.3825; 1.3825; 1.3599; 1.3632; 1.3540;...
1.3605; 1.3522; 1.3285; 1.3522; 1.3522; 1.3522; 1.3757;...
1.3757; 1.3757; 1.3757; 1.3757; 1.3757; 1.3757; 1.3757;...
1.3757; 1.3975; 1.4178; 1.4368; 1.4545; 1.4580; 1.4806;...
1.4806; 1.4806; 1.4979; 1.5141; 1.5141; 1.5088; 1.5088;...
1.5230; 1.5364; 1.5488; 1.5565; 1.5482; 1.5482; 1.5482;...
1.5482; 1.5482; 1.5353; 1.5353; 1.5353; 1.5353; 1.5511;...
1.5511; 1.5596; 1.5596; 1.5547; 1.5547; 1.5547; 1.5547;...
1.5547; 1.5547; 1.5694; 1.5694; 1.5694; 1.5547; 1.5457;...
1.5604; 1.5569; 1.5569; 1.5422; 1.5542; 1.5671; 1.5671;...
1.5671; 1.5671; 1.5671; 1.5792; 1.5792; 1.5792; 1.5903;...
1.5903; 1.5513; 1.5513; 1.5513; 1.5611; 1.5702; 1.5702;...
1.5702; 1.5702; 1.5702; 1.5702; 1.5702; 1.5611; 1.5650;...
1.5650; 1.5650; 1.5650; 1.5650; 1.5657; 1.5727; 1.5721;...
1.5721; 1.5630; 1.5593; 1.5593; 1.5593; 1.5593; 1.5593;...
1.5593; 1.5649; 1.5664; 1.5664; 1.5937; 1.5937; 1.5937;...
1.5918; 1.5889; 1.5889; 1.5738; 1.5738; 1.5738; 1.5738;...
1.5738; 1.5498; 1.5392; 1.5229; 1.5229; 1.5229; 1.5285;...
1.5285; 1.5285; 1.5377; 1.5438; 1.5438; 1.5285; 1.5269;...
1.5199; 1.5199; 1.5199; 1.5199; 1.5173; 1.5173; 1.5173;...
1.5173; 1.5173; 1.5280; 1.5335; 1.5335; 1.5433; 1.5433;...
1.5727; 1.5727; 1.5727; 1.5727; 1.5774; 1.5952; 1.5949;...
1.5949; 1.5546; 1.5546; 1.5546; 1.5546; 1.5546; 1.6011;...
1.6011; 1.6011; 1.6011; 1.5931; 1.5973; 1.5973; 1.5893;...
1.5893; 1.5468; 1.5468; 1.5650; 1.5738; 1.5588; 1.5588;...
1.5588; 1.5674; 1.5625; 1.5452; 1.5452; 1.5452; 1.5684;...
1.5684; 1.5789; 1.5938; 1.5903; 1.5903; 1.5994; 1.6078;...
1.6078; 1.6537; 1.6537; 1.6449; 1.6449; 1.6353; 1.6449;...
1.6578; 1.6685; 1.6685; 1.6685; 1.6657; 1.6525; 1.6525;...
1.6571; 1.6571; 1.6571; 1.6468; 1.6468; 1.6468; 1.6647;...
1.6647; 1.6637; 1.6625; 1.6652; 1.6652; 1.6193; 1.6297;...
1.6297; 1.6297; 1.6377; 1.6390; 1.6368; 1.6368; 1.5987;...
1.5873; 1.5873; 1.5948; 1.5990; 1.6206; 1.6206; 1.6489;...
1.6365; 1.6365; 1.6365; 1.6365; 1.6365; 1.6365; 1.6136;...
1.6011; 1.5934; 1.5934; 1.6214; 1.6214; 1.6214; 1.6214;...
1.6342; 1.6382; 1.6198; 1.5935; 1.5952; 1.5891; 1.5891;...
1.5891; 1.5769; 1.6152; 1.6152; 1.6021; 1.6021; 1.6483;...
1.5966; 1.5945; 1.5945; 1.5945; 1.5945; 1.6001; 1.6047;...
1.6278; 1.6638; 1.6638; 1.6513; 1.6101; 1.6307; 1.6455;...
1.6440; 1.6440; 1.6564; 1.6564; 1.6491; 1.6591; 1.6411;...
1.6022; 1.5788; 1.5504; 1.5628; 1.5785; 1.5934; 1.5934;...
1.5563; 1.5594; 1.5632; 1.5720; 1.5673; 1.5673; 1.5700;...
1.5700; 1.5700; 1.5623; 1.5797; 1.5951; 1.6013; 1.6013;...
1.6013; 1.6013; 1.6013; 1.6013; 1.5903; 1.5903; 1.5903;...
1.5903; 1.5903; 1.5903; 1.5868; 1.5947; 1.5718; 1.5718;...
1.5718; 1.5862; 1.6073; 1.6073];
top = 501 - 1 - runlen;
for j = 1:top
  accum = 0;
  for k = 1:runlen
  % cleverer code would initialize accum array then replace 1 value at a time
    accum = accum + enerIgn1(j+k-1);
  end % for k
  enerApp1(j) = accum/runlen;
end % for j
figure(13)
plot(1:501,enerIgn1,'b',1+runlen/2:top+runlen/2,enerApp1,'r')
title('smooth2EnerIgn1')
xlabel('time (# collisions)')
ylabel('energy ignorance before mixing starts')
enerIgn2=...	% 72 lines
[1.6073; 1.6073; 1.6054;...
1.5763; 1.5763; 1.5796; 1.5796; 1.5701; 1.5701; 1.5701;...
1.5426; 1.5517; 1.5517; 1.5601; 1.6010; 1.6010; 1.6110;...
1.6110; 1.6110; 1.6571; 1.6571; 1.6571; 1.7058; 1.6921;...
1.6921; 1.7245; 1.7200; 1.7519; 1.7545; 1.7483; 1.7396;...
1.6843; 1.6843; 1.7002; 1.6783; 1.6883; 1.6783; 1.6745;...
1.6511; 1.6661; 1.6661; 1.6661; 1.6564; 1.6564; 1.6564;...
1.6564; 1.6564; 1.6564; 1.6754; 1.6754; 1.6754; 1.6866;...
1.6784; 1.6774; 1.6942; 1.7034; 1.7034; 1.7122; 1.7069;...
1.7364; 1.7364; 1.7364; 1.7477; 1.7477; 1.7477; 1.7477;...
1.7384; 1.7301; 1.7321; 1.7321; 1.6935; 1.6996; 1.7078;...
1.7197; 1.7188; 1.7188; 1.7122; 1.7122; 1.7122; 1.7076;...
1.7076; 1.7390; 1.7361; 1.7416; 1.7374; 1.7463; 1.7374;...
1.7374; 1.7299; 1.7220; 1.7153; 1.6959; 1.6959; 1.6959;...
1.6959; 1.6959; 1.7102; 1.7070; 1.7070; 1.7228; 1.7319;...
1.7345; 1.7395; 1.7395; 1.7310; 1.7310; 1.7398; 1.7429;...
1.7817; 1.7980; 1.8116; 1.8116; 1.8163; 1.8393; 1.8483;...
1.8393; 1.8288; 1.8288; 1.8197; 1.8293; 1.8293; 1.8255;...
1.8255; 1.8255; 1.8295; 1.8295; 1.8304; 1.8304; 1.8388;...
1.8388; 1.8408; 1.8408; 1.8463; 1.8463; 1.8463; 1.8463;...
1.8408; 1.8481; 1.8481; 1.8456; 1.8465; 1.8465; 1.8465;...
1.8465; 1.8428; 1.8353; 1.8361; 1.8361; 1.8361; 1.8354;...
1.8073; 1.8327; 1.8333; 1.8333; 1.8308; 1.8369; 1.8389;...
1.7934; 1.7900; 1.7950; 1.7942; 1.7615; 1.7615; 1.8145;...
1.8145; 1.8118; 1.8118; 1.8118; 1.8099; 1.8099; 1.8238;...
1.8238; 1.8152; 1.8152; 1.8190; 1.8120; 1.8514; 1.8533;...
1.8527; 1.8527; 1.8541; 1.8541; 1.8527; 1.8540; 1.8484;...
1.8519; 1.8448; 1.8376; 1.8410; 1.8410; 1.8363; 1.8363;...
1.8363; 1.8298; 1.8298; 1.8224; 1.7900; 1.7769; 1.7811;...
1.7811; 1.7697; 1.7557; 1.7636; 1.7591; 1.7591; 1.7904;...
1.7904; 1.7904; 1.7910; 1.7786; 1.7786; 1.7515; 1.7555;...
1.7736; 1.7682; 1.7682; 1.7691; 1.7691; 1.7691; 1.7587;...
1.7587; 1.7560; 1.7671; 1.7772; 1.7801; 1.7788; 1.7284;...
1.7293; 1.7198; 1.7158; 1.7146; 1.6848; 1.6956; 1.6963;...
1.6865; 1.7001; 1.7333; 1.7340; 1.7394; 1.7789; 1.7937;...
1.8022; 1.8098; 1.8046; 1.7867; 1.7867; 1.7960; 1.7960;...
1.7960; 1.7919; 1.7907; 1.7907; 1.7907; 1.7907; 1.7986;...
1.7958; 1.7958; 1.7931; 1.7931; 1.7924; 1.7924; 1.7852;...
1.7852; 1.7898; 1.7819; 1.7886; 1.7863; 1.7776; 1.8022;...
1.8097; 1.8175; 1.8175; 1.8048; 1.8078; 1.8076; 1.8076;...
1.8076; 1.8178; 1.8184; 1.8184; 1.8184; 1.8080; 1.8089;...
1.8056; 1.8033; 1.7978; 1.7978; 1.8136; 1.8136; 1.8148;...
1.8212; 1.8212; 1.8212; 1.8190; 1.8225; 1.8158; 1.8069;...
1.8134; 1.8128; 1.8128; 1.8128; 1.8128; 1.8049; 1.8049;...
1.8198; 1.8188; 1.8188; 1.8188; 1.8209; 1.8209; 1.8209;...
1.8209; 1.8206; 1.8206; 1.8209; 1.8108; 1.8174; 1.8121;...
1.8121; 1.7868; 1.7636; 1.7690; 1.7636; 1.7703; 1.7703;...
1.7688; 1.7740; 1.7740; 1.7740; 1.8207; 1.8207; 1.8224;...
1.8154; 1.8052; 1.8052; 1.7985; 1.7996; 1.7986; 1.8085;...
1.8085; 1.8085; 1.8098; 1.8180; 1.8210; 1.8225; 1.8172;...
1.8181; 1.8181; 1.8188; 1.8269; 1.8202; 1.8249; 1.8233;...
1.8233; 1.8233; 1.8068; 1.8141; 1.8141; 1.8141; 1.8068;...
1.8074; 1.8097; 1.8097; 1.8082; 1.8096; 1.8003; 1.8003;...
1.8003; 1.7966; 1.8145; 1.7950; 1.7950; 1.7950; 1.7950;...
1.8049; 1.7950; 1.7950; 1.7950; 1.7950; 1.8010; 1.8124;...
1.8085; 1.8085; 1.8240; 1.8206; 1.8206; 1.8200; 1.8116;...
1.8116; 1.8473; 1.8800; 1.8800; 1.8800; 1.8800; 1.8800;...
1.8874; 1.8824; 1.8824; 1.9067; 1.9067; 1.9067; 1.9067;...
1.9019; 1.8891; 1.8891; 1.8674; 1.8763; 1.8763; 1.8727;...
1.8760; 1.8760; 1.8850; 1.8850; 1.8876; 1.8968; 1.9043;...
1.9043; 1.9067; 1.9126; 1.9126; 1.9093; 1.9103; 1.9103;...
1.9326; 1.9326; 1.9318; 1.9318; 1.9302; 1.9302; 1.9293;...
1.9293; 1.9313; 1.9313; 1.9313; 1.9313; 1.9319; 1.9327;...
1.9267; 1.9335; 1.9314; 1.9350; 1.9164; 1.9164; 1.9164;...
1.8789; 1.8781; 1.8791; 1.8768; 1.8782; 1.8782; 1.8759;...
1.8759; 1.8742; 1.8729; 1.8729; 1.8729; 1.8729; 1.8772;...
1.8230; 1.8230; 1.8096; 1.8043; 1.8043; 1.8069; 1.8083;...
1.8083; 1.8060; 1.8060; 1.7975; 1.8151; 1.8038; 1.8182;...
1.8247; 1.8069; 1.8069; 1.8069; 1.8069; 1.8069; 1.7985;...
1.7964; 1.8224; 1.8084; 1.8172; 1.8172; 1.8148; 1.8115;...
1.8115; 1.8362; 1.8690; 1.8675; 1.8602; 1.8612; 1.8641;...
1.8675; 1.8713; 1.8724; 1.8750; 1.8750; 1.8750; 1.8764];
top = 500 - 1 - runlen;
for j = 1:top
  accum = 0;
  for k = 1:runlen
  % cleverer code would initialize accum array then replace 1 value at a time
    accum = accum + enerIgn2(j+k-1);
  end % for k
  enerApp2(j) = accum/runlen;
end % for j
figure(14)
plot(1:500,enerIgn2,'b',1+runlen/2:top+runlen/2,enerApp2,'r')
title('smooth2EnerIgn2')
xlabel('time (# collisions)')
ylabel('energy ignorance after mixing starts')
figure(15)
plot(1:1001,[enerIgn1;enerIgn2],'b',1+runlen/2:500-runlen/2, enerApp1,'r',502+runlen/2:1000-runlen/2,enerApp2,'r')

% function molData = presSim(N,T,r,M,ab)	THM			101124
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = moleculeInit(N,T,startRange,r,M)
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
%	molDatAverage = molDataHisto(molData,[]);
% Differs from statSim() by plotting rates of change of momentum per unit side;
% also avoids first 100 collisions.
function molData = presSim(N,T,r,M,ab)

% initialize
startRange = [0,0,5,0,0,40,50];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
molData = moleculeInit(N,T,startRange,r,M)	% display test only
% molDatAverage = molDataHisto(molData,[]);
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
oldTime = 0;
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
totTop = 0; numTop = 0;	% initialize momentum counting accumulators
totBot = 0; numBot = 0;	% initialize momentum counting accumulators
totRgt = 0; numRgt = 0;	% initialize momentum counting accumulators
totLft = 0; numLft = 0;	% initialize momentum counting accumulators
pressure = [];
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % distinguish collision on top, bottom, right, left sides and accumulate
    if deltaPx1==0
      if deltaPy1<0
	totTop = totTop + deltaPy1;
	numTop = numTop + 1;
      else % deltaPy1>=0 and the == doesn't matter
	totBot = totBot + deltaPy1;
	numBot = numBot + 1;
      end % if deltaPy1<0
    end %if deltaPx1==0
    if deltaPy1==0
      if deltaPx1<0
	totRgt = totRgt + deltaPx1;
	numRgt = numRgt + 1;
      else % deltaPx1>=0 and the == doesn't matter
	totLft = totLft + deltaPx1;
	numLft = numLft + 1;
      end % if deltaPx1<0
    end %if deltaPy1==0
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if mod(happIndx,stepChunk)==0
      elapsed = toc;				% elapsed time
      chunkNo = chunkNo + 1;
      if firstChunk
	firstChunk = false;
      end % if firstChunk
      % molDatAverage = molDataHisto(molData,molDatAverage);
      beep
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/n '],'s');
      if isempty(reply)
          reply = 'y';
      end
      if reply~='y'	% not 'y'
	pressure = [pressure;[totTop/a/t,totBot/a/t,totRgt/b/t,totLft/b/t]];
	pressure				% display pressure
	sizePressure = size(pressure);
	P = sizePressure(1);
	figure(16)			% hardwired: will overwrite any fig 16
	plot(1:P,pressure(:,1),'b',1:P,pressure(:,2),'g',1:P,pressure(:,3),'r',1:P,pressure(:,4),'m')
	title('pressures on four sides')
	xlabel(['chunks of ' num2str(stepChunk) ' intermol coll'])
	ylabel('pressure (meV/nm^2)')
	% NB 3-D pressure is meV/nm^3 = 1.6 atm (1 atm = 101325 Pa (J/m^3) and
	%  conversion meV = 1.602 176 53*10^-22 J)
	legend('top','bottom','right','left','Location','East')
        break
      else tic					% elapsed time
	t = time - oldTime;
	a = ab(1);
	b = ab(2);
	pressure = [pressure;[totTop/a/t,totBot/a/t,totRgt/b/t,totLft/b/t]];
	totTop = 0; numTop = 0;	% reset momentum counting accumulators
	totBot = 0; numBot = 0;	% reset momentum counting accumulators
	totRgt = 0; numRgt = 0;	% reset momentum counting accumulators
	totLft = 0; numLft = 0;	% reset momentum counting accumulators
	oldTime = time;
      end % if reply
    end % if mod()
  end % if collision
end % while(true)

% function molData = press2Sim(N,T,r,M,ab)	THM			101118
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% N,T,r,M 2*1 column vectors to accommodate 2 gases
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = molecuReInit(N,T,startRange,r,M)	fixed initial config.
%    or	molData1 = moleculeInit(N(1),T(1),startRange1,r(1),M(1))
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0,mix);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
% Differs from tmptureSim by tracking the pressures of the two gases before and
% after mixing: see ../newsim1/presSim(). Does not record energy stats.
function molData = press2Sim(N,T,r,M,ab)

% initialize
startRange1 = [0,0,5,0,0,40,50];
startRange2 = [0,0,5,0,0,-50,-40];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
% use one or other of the following two lines:
%molData = molecuReInit(N,T,startRange,r,M)	% display test only
molData1 = moleculeInit(N(1),T(1),startRange1,r(1),M(1))   % display test only
molData2 = moleculeInit(N(2),T(2),startRange2,r(2),M(2))   % display test only
molData = [molData1;molData2];
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
oldTime = 0.0;
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
mix = 0;		% initially don't mix the gases in top and bottom
totTop = 0; numTop = 0;	% initialize momentum counting accumulators
totBot = 0; numBot = 0;	% initialize momentum counting accumulators
pressure = [];
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0,mix);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % distinguish collision on top for gas 1, bottom for gas 2 and accumulate
    if deltaPx1==0
      if deltaPy1<0
	if v1<=N(1)			% molecule belongs to upper gas (1)
	  totTop = totTop + deltaPy1;
	  numTop = numTop + 1;
	end % if v1<=N(1)
      else % deltaPy1>=0 and the == doesn't matter
	if v1>N(1)			% molecule belongs to lower gas (2)
	  totBot = totBot + deltaPy1;
	  numBot = numBot + 1;
	end % if v1>N(1)
      end % if deltaPy1<0
    end %if deltaPx1==0
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if mod(happIndx,stepChunk)==0
      elapsed = toc;			% elapsed time
      chunkNo = chunkNo + 1;
      beep
      avgEnergy1 = sum(molData(1:N(1),7))/N(1);
      avgEnergy2 = sum(molData(N(1)+1:N(1)+N(2),7))/N(2);
      energyX = (molData(:,3).^2./(2*10.36*molData(:,6)));
      avgEnergyX = sum(energyX)/(N(1)+N(2));	% check equipartition
      energyY = (molData(:,4).^2./(2*10.36*molData(:,6)));
      avgEnergyY = sum(energyY)/(N(1)+N(2));
      disp(['Average energies: gas 1 ' num2str(avgEnergy1) ', gas 2 ' ...
	num2str(avgEnergy2) ', x partition ' num2str(avgEnergyX) ...
	', y partition ' num2str(avgEnergyY)])
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/m[mix]/n '],'s');
      if isempty(reply)
	reply = 'y';
      end % if isempty
      if reply=='m'			% toggle mixing gases
	mix = ~mix;
	reply = 'y';
      end % if reply==
      if reply~='y'	% not 'y'
	pressure = [pressure;[-totTop/a/t,totBot/a/t]];
	pressure				% display pressure
	sizePressure = size(pressure);
	P = sizePressure(1);
        break
      else tic					% elapsed time
	t = time - oldTime;
	a = ab(1);
	b = ab(2);
	pressure = [pressure;[-totTop/a/t,totBot/a/t]];
	sizePressure = size(pressure);
	P = sizePressure(1);
	figure(17)			% hardwired: will overwrite any fig 17
	plot(1:P,pressure(:,1),'b',1:P,pressure(:,2),'g')
	title('pressures of upper and lower gases')
	xlabel(['chunks of ' num2str(stepChunk) ' intermol coll'])
	ylabel('pressure (meV/nm^2)')
	% NB 3-D pressure is meV/nm^3 = 1.6 atm (1 atm = 101325 Pa (J/m^3) and
	%  conversion meV = 1.602 176 53*10^-22 J)
	legend('upper','lower','Location','West')
	totTop = 0; numTop = 0;	% reset momentum counting accumulators
	totBot = 0; numBot = 0;	% reset momentum counting accumulators
	oldTime = time;
      end % if reply~=
    end % if mod()
  end % if collision
end % while(true)

% function molData = chmpotSim(N,T,r,M,ab)	THM			101118
% N molecules each of radius r (nm), mass M (AMU) at temperature T (degK)
% N,T,r,M 2*1 column vectors to accommodate 2 gases
% Helium 0.049nm (1), 4AMU; Neon 0.051 m (1), 20AMU; Argon 0.088nm (1.8), 40AMU)
% find; molecules in box [-ab(1),ab(1)]x[-ab(2),ab(2)] nanometers
%	(initially symmetrical, stored abcd [abcd(1),abcd(2)]x[abcd(3),abcd(4)])
% Return molData
% uses	molData = molecuReInit(N,T,startRange,r,M)	fixed initial config.
%    or	molData1 = moleculeInit(N(1),T(1),startRange1,r(1),M(1))
%	[delTime,centre] = molCollWallWhen(cb,p,M,r,ab);
%	[delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M,M,r,r);
%	p1= molCollWallCountBot(p1b,1,1,abcd,r1,c1,0,mix);
%	[p1,p2] = moleculeCollide(M,M,p1b,p2b,0,C12);
% Differs from tmptureSim, press2Sim by tracking the centres of masses of the
% two gases before and after mixing
function molData = chmpotSim(N,T,r,M,ab)

% initialize
startRange1 = [0,0,5,0,0,40,50];
startRange2 = [0,0,5,0,0,-50,-40];
%	      l h s x y  l h	s is stdev in nm of spread along line of angle
%	      ang   cntr	from ang_l to ang_h centred at x,y
%	     ----pos---  mom	mom_l to mom_h is angle of motion direction
% for box -20:20*-10:10 try exactly 45 degrees from exactly (0,0)
% generates molData[N * cx,cy,px,py,r,M,E]
%			 1  2  3  4 5 6 7
% L: nm   T: ps  M: meV(ps/nm)^2   L/T: nm/ps   ML/T: meV(ps/nm)
% use one or other of the following two lines:
%molData = molecuReInit(N,T,startRange,r,M)	% display test only
molData1 = moleculeInit(N(1),T(1),startRange1,r(1),M(1))   % display test only
molData2 = moleculeInit(N(2),T(2),startRange2,r(2),M(2))   % display test only
molData = [molData1;molData2];
oldMolData = molData;
sizeMolData = size(molData);
noVals = sizeMolData(1);
abcd(1) = -ab(1); abcd(2) = ab(1); abcd(3) = -ab(2); abcd(4) = ab(2);
happIndx = 0;		% keep track of number of steps
time = 0.0;		% keep track of time
stepChunk = 100;	% iterate this many steps then check if more wanted
firstChunk = true;	% to detect first  collisions at resume time
%stepChunk = 1;		% test only
chunkNo = 0;
mix = 0;		% initially don't mix the gases in top and bottom
CofM = [];
tic					% elapsed time

% Loop starts here
while(true)
  collision = false;
  % find time of next collision
  record = [];				% record of collision times, results
  for v1 = 1:noVals
    c1b = [molData(v1,1); molData(v1,2)];
    p1b = [molData(v1,3); molData(v1,4)];
    M1 = molData(v1,6);
    r1 = molData(v1,5);
    [delTime,centre] = molCollWallWhen(c1b,p1b,M1,r1,abcd);  % display test only
    record = [record;[delTime,v1,v1,centre(1),centre(2),centre(1),centre(2)]];
    for v2 = v1+1:noVals
      c2b = [molData(v2,1); molData(v2,2)];
      p2b = [molData(v2,3); molData(v2,4)];
      M2 = molData(v2,6);
      r2 = molData(v2,5);
      [delTime,cntr1,cntr2] = moleculeCollWhen(c1b,c2b,p1b,p2b,M1,M2,r1,r2);
      record = [record;[delTime,v1,v2,cntr1(1),cntr1(2),cntr2(1),cntr2(2)]];
    end % for v2
  end % for v1

  %sqrt choose and execute next collision
  next = minrows(record);	% smallest delTime wins  display test only
  time = time + next(1);
  v1 = next(2);
  v2 = next(3);
  c1 = [next(4);next(5)];
  c2 = [next(6);next(7)];
  deltaX1 = c1(1) - molData(v1,1);
  deltaY1 = c1(2) - molData(v1,2);
  deltaX2 = c2(1) - molData(v2,1);
  deltaY2 = c2(2) - molData(v2,2);
  r1 = molData(v1,5);
  r2 = molData(v2,5);
  M1 = molData(v1,6);
  M2 = molData(v2,6);
  p1b = [molData(v1,3); molData(v1,4)];
  p2b = [molData(v2,3); molData(v2,4)];
  if v1==v2				% collision is with wall
    p1 = molCollWallCountBot(p1b,1,1,abcd,r1,c1,0,mix);
    deltaPx1 = p1(1) - p1b(1);		% change of momentum at wall
    deltaPy1 = p1(2) - p1b(2);
    % molData for v1 is common to both these branches: done afterwards
  else % v1~=v2
    collision = true;			% test only: used to update afterColl
    C12 = c2 - c1;
    [p1,p2] = moleculeCollide(M1,M2,p1b,p2b,0,C12);	% display test only
    molData(v2,1) = c2(1);
    molData(v2,2) = c2(2);
    molData(v2,3) = p2(1);
    molData(v2,4) = p2(2);
    molData(v2,7) = (p2(1).^2 + p2(2).^2)/(2*M2*10.36);		% energy
  end % if v1==v2
  molData(v1,1) = c1(1);
  molData(v1,2) = c1(2);
  molData(v1,3) = p1(1);
  molData(v1,4) = p1(2);
  molData(v1,7) = (p1(1).^2 + p1(2).^2)/(2*M1*10.36);		% energy
  for v = 1:noVals			% update all other molecules' positions
    if v==v1 | v==v2, continue, end % if
    deltaX = molData(v,3)/(molData(v,6)*10.36)*next(1);
    deltaY = molData(v,4)/(molData(v,6)*10.36)*next(1);
    molData(v,1) = molData(v,1) + deltaX;	% update position
    molData(v,2) = molData(v,2) + deltaY;
  end % for v
  if collision
    happIndx = happIndx + 1;
    if mod(happIndx,stepChunk)==0
      elapsed = toc;			% elapsed time
      chunkNo = chunkNo + 1;
      beep
      avgEnergy1 = sum(molData(1:N(1),7))/N(1);
      avgEnergy2 = sum(molData(N(1)+1:N(1)+N(2),7))/N(2);
      energyX = (molData(:,3).^2./(2*10.36*molData(:,6)));
      avgEnergyX = sum(energyX)/(N(1)+N(2));	% check equipartition
      energyY = (molData(:,4).^2./(2*10.36*molData(:,6)));
      avgEnergyY = sum(energyY)/(N(1)+N(2));
      disp(['Average energies: gas 1 ' num2str(avgEnergy1) ', gas 2 ' ...
	num2str(avgEnergy2) ', x partition ' num2str(avgEnergyX) ...
	', y partition ' num2str(avgEnergyY)])
      reply = input([num2str(elapsed) ' elapsed, ' num2str(time) ' ps.: another ' num2str(stepChunk) ' collisions? y/m[mix]/n '],'s');
      if isempty(reply)
	reply = 'y';
      end % if isempty
      if reply=='m'			% toggle mixing gases
	mix = ~mix;
	reply = 'y';
      end % if reply==
      if reply~='y'	% not 'y'
        break
      else tic					% elapsed time
	N1 = N(1);
	N2 = N(2);
	CofM1 = [sum(molData(1:N1,1)),sum(molData(1:N1,2))]/N1;
	CofM2 = [sum(molData(N1+1:N1+N2,1)),sum(molData(N1+1:N1+N2,2))]/N2;
	CofM = [CofM;[CofM1,CofM2]];
	sizeCofM = size(CofM);
	C = sizeCofM(1);
	figure(18)			% hardwired: will overwrite any fig 18
%	plot(1:C,CofM(:,1),'b',1:C,CofM(:,2),'g',1:C,CofM(:,3),'r',1:C,CofM(:,4),'m')
	plot(1:C,CofM(:,2),'b',1:C,CofM(:,4),'r')
	title('centres of mass of upper and lower gases')
	xlabel(['chunks of ' num2str(stepChunk) ' intermol coll'])
	ylabel('centres of mass (nm)')
	% NB 3-D pressure is meV/nm^3 = 1.6 atm (1 atm = 101325 Pa (J/m^3) and
	%  conversion meV = 1.602 176 53*10^-22 J)
%	legend('upper x','upper y','lower x','lower y')
	legend('upper y','lower y','Location','Best')
      end % if reply~=
    end % if mod()
  end % if collision
end % while(true)