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

---

%MATLABpak09cIII:
%
%	oneMolecAvg.m				oneMolecExpAvg.m
%

% function [sumAC,MTBC] = oneMolecAvg(timeStep,noSteps,noIts) THM	110523
%(adapted 110606 from autocorrgas/oneMolecAvg.m, simplified almost back to orig)
% input: integer  timeStep  gives a fixed time between "collisions",
% i.e., before the "velocity" randomly changes, but, depending on internal
% parameter  mult , step size will be multiplied by  mult  1/mult of the time
%  integer  noSteps  gives number of these steps.
% output: MTBC is mean time between collisions, from TTBC/noColl
%  sumAC is sum of all terms of autocorrelation: use ac = autocorrel(velseq)
%   noIts  times and accumulate into  acc
% Also plots the two autocorrelations
function [sumAC,MTBC] = oneMolecAvg(timeStep,noSteps,noIts)
mult = 1;		% internal parameter: good results only for mult = 1
TTBC = 0;
for it = 1:noIts
  velIndx = 1;
  for j = 1:noSteps-1		% calculate random velocities after each "coll."
    vs = randn(1);			% new vel after each coll: normal
    if mod(j,mult)==0			% step size * mult, 1/mult of the time
      for k = 1:mult*timeStep		% mult*timestep  for one j
	velseq(:,velIndx) = vs;		% repeat for size of mult*timeStep
	velIndx = velIndx + 1;
      end % for k
      TTBC = TTBC + mult*timeStep;
    else				% doubling, added 110523
      for k = 1:timeStep		% timeStep  for all other j
	velseq(:,velIndx) = vs;		% repeat for size of timeStep
	velIndx = velIndx + 1;
      end % for k
      TTBC = TTBC + timeStep;
    end % if mod()			% doubling, added 110523
  end % for j
  ac = autocorrelND(velseq);	% autocorrel normalized to 1 for lag 0
  if it==1
    acc = ac;			% initialize
  else
    acc = acc + ac;		% sum of autocorrelations
  end % if it
end % for it
acc = acc/noIts;		% average of autocorrelations
MTBC = TTBC/(noSteps-1)/noIts;	% goes with  mult
sizeV = size(velseq);
sizeV2 = sizeV(2);
plot(1-sizeV2:sizeV2-1,acc)
sumAC = sum(acc);

% function [sumAC,MTBC] = oneMolecExpAvg(intercolldist,noIts) THM	110607
%(adapted 110607 from oneMolecAvg.m, to use exponential distr of time intervals)
%(now these notes differ from the code: intercolldist generated internally)
% input: intercolldist is array of exponentially distributed time intervals
%	 whose size is the number of timesteps: must be generated externally
%	 to be reused noIts times---intercolldist = intercoll(1,noSteps,tau);
%	 otherwise get different sized velseq and hence ac
%  noIts: is number of iterations for averaging
% output: MTBC is mean time between collisions, from TTBC/noColl
%  sumAC is sum of all terms of autocorrelation: use ac = autocorrel(velseq)
%   noIts  times and accumulate into  acc
% Also plots the autocorrelation
function [sumAC,MTBC] = oneMolecExpAvg(tau,noSteps,noIts)
% sizeIntercolldist = size(intercolldist);
% noSteps = sizeIntercolldist(2);		% intercolldist is row vector
TTBC = 0;
totSumAC1 = 0;
totSumAC2 = 0;
for it = 1:noIts
  velIndx = 1;
  intercolldist = intercoll(1,noSteps,tau);
  for j = 1:noSteps		% calculate random velocities after each "coll."
    vs = randn(1);			% new vel after each coll: normal
    timeStep = intercolldist(j);	% exponentially distributed
    for k = 1:timeStep;		% timeStep  for all other j
      velseq(:,velIndx) = vs;		% repeat for size of timeStep
      velIndx = velIndx + 1;
    end % for k
    TTBC = TTBC + timeStep;
  end % for j
  ac = autocorrelND(velseq);	% autocorrel normalized to 1 for lag 0
% display(['oneMolecExpAvg: size(velseq) = ' num2str(size(velseq)) ', size(ac) = ' num2str(size(ac))])	% test only
  sizeV = size(velseq);
  sizeV2 = sizeV(2);
  sumAC1 = sum(ac(sizeV2:2*sizeV2-1));
% display(['oneMolecExpAvg: sumAC = ' num2str(sumAC1)])
  totSumAC1 = totSumAC1 + sumAC1;
  totSumAC2 = totSumAC2 + sumAC1^2;
% if it==1
%   acc = ac;			% initialize
% else
%   acc = acc + ac;		% sum of autocorrelations
% end % if it
end % for it
% acc = acc/noIts;		% average of autocorrelations
MTBC = TTBC/noSteps/noIts;
sumAC = totSumAC1/noIts;
stdAC = sqrt((totSumAC2 - totSumAC1^2/noIts)/noIts);
display(['oneMolecExpAvg: mean sumAC = ' num2str(sumAC) ', stDev sumAC = ' num2str(stdAC)])	% test only
% sizeV = size(velseq);
% sizeV2 = sizeV(2);
% plot(1-sizeV2:sizeV2-1,acc)
% sumAC = sum(acc);
% sumAC = sum(acc(sizeV2:2*sizeV2-1));	% non-negative terms only
% display(['oneMolecExpAvg: size(velseq) = ' num2str(sizeV2) ', size(ac) = ' num2str(size(acc))])	% test only