% ### EXcoffee.m ###    11.16.14    {C. Bergevin}
% [modified version of EXrandomWalk2D.m motivated by the problem shown in
% Giordano (1997) Fig.7.18ff]
% **NOTE**: There is a minor bug in this version such that it is possible for
% some 'cream' to leave the 'cup' (despite the specifed boundary conditions)
clear;
% -------------
N= 10;      % one plus sqrt of total # of (independent) walkers (each starts at unique x,y point about origin)
M= 300;      % Total # of steps for each walker
method= 2;  % see comments above
BND= 10;    % bounding limits for initial grid of walkers at t=0
axisLim= 100;   % size of coffee cup
diffC= 1;   % diffusion const. (i.e., scaling factor for step size)
framerate= 1/30;    % pause length [s] for animation
Sgrid= 8;           % grid spacing for entropy calculation
% -------------
% +++
space= (2*BND)/N;
[X,Y]= meshgrid(-BND:space:BND,-BND:space:BND);
E= size(X,1); % # of elements
SgridX= linspace(-axisLim,axisLim,Sgrid);   % set grid bounds for entropy calc.
SgridY= linspace(-axisLim,axisLim,Sgrid);

figure(1); clf; grid on; xlabel('x-postion'); ylabel('y-postion');
% visualize before onset?
if (1==1), plot(X,Y,'ko','MarkerSize',5); axis([-axisLim axisLim -axisLim axisLim]);   end
% +++
% To do
% - apply boundary condition (i.e., ensure no steps past walls)
% - fix entropy calc. (i.e., if prob.=0??)
for r= 1:M
    
    if method==1
        % random L/R and U/D step with equal probability
        tempX= rand(E,E);   tempY= rand(E,E);  % determine random vals.
        temp2X= tempX<0.5;  temp2Y= tempY<0.5;  % determine L vs R and U vs D
        X(temp2X)= X(temp2X)+1; X(~temp2X)= X(~temp2X)-1;
        Y(temp2Y)= Y(temp2Y)+1; Y(~temp2Y)= Y(~temp2Y)-1;
    else
        % sample step from normal distribution
        stepX= randn(E,E); stepY= randn(E,E);
        X= X+ diffC*stepX;  Y= Y+ diffC*stepY;
        % verify step is not past walls; if so, bounce back in opposite direction
        [aa,bb]= find(abs(X)>axisLim); [cc,dd]= find(abs(Y)>axisLim); 

        % +++ --> correct for points that have moved past the walls
        % not quite right, but kinda works
        X(aa,bb)= X(aa,bb)-2*diffC*stepX(aa,bb); Y(cc,dd)= Y(cc,dd)-2*diffC*stepY(cc,dd);
        
        % more right (I think), but doesn't work
        %X(aa,bb)= sign(X(aa,bb))*2*axisLim-X(aa,bb); Y(cc,dd)= Y(cc,dd)-2*diffC*stepY(cc,dd);
        
        % uncomment to allow for flagging when 'cream' leaves the cup
        if(max(abs(X(:))>axisLim)), return;  end
    end
    % visualize
    figure(1)
    plot(X,Y,'ko','MarkerSize',5); axis([-axisLim axisLim -axisLim axisLim]); pause(framerate);
    % do binning to determine 'probability' distribution 
    histS= hist2(X(:),Y(:),SgridX,SgridY)/E^2; % use external function hist2.m; and normalize to a probability
    histS= histS(:);    % convert to a single column vector
    zeroI= ~histS==0; % need to filter out states with zero elements so to avoid computational error (since 0*log(0)= NaN)
    S(r)= -sum(histS(zeroI).*log(histS(zeroI))); % calculate entropy (S) 
end;

figure(2)
plot(S,'LineWidth',2); hold on; grid on;
xlabel('time step'); ylabel('entropy');