% example mixing problem
% by Justin Dunmyre, 6/2012
% modified by Gavin LaRose, 6/28/2012
% this demonstration models the number of black balls in a box as
%    a mixing problem and carries out the simulation of the number
%    of balls in the box, showing the number of balls in the box
%    at any time and a graph of the number of balls and the value
%    of the differential equation model that advances at the same
%    rate.
%
% (c)2012 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.

function Example_Mixing_Problem

%----------------------------------------------------------------------
% configuration variables

    % the number of black balls added each time unit
    blackInPerFrame=3;
    % the size of the hole through which balls leave the box, in
    %     ball units
    holeSize = 4;
    % then we put in invisible white balls to fill out the volume
    whiteInPerFrame = holeSize - blackInPerFrame;
    % length of simulation, in time units
    simulationLength=200;
    % we plot every frameskip steps in the simulation
    frameskip=2;

    % initial number of black dots in the box
    numBlack=0;
    % "volume" of the box, in number of balls; should be divisible
    %    by the boxWidth, below.  good values are 50, 100 and 200
    volume=100;
    % initial number of white balls (that is, empty spaces) in box
    initWhite=volume-numBlack;
    
    boxWidth=10;
    boxHeight=volume/boxWidth;

%----------------------------------------------------------------------
% generate the demo
   
    % solution to the modeling differential equation; this gives
    %    the number of black balls in the box
    [t,theorySol]=ode45(@mixingODE,1:simulationLength,numBlack);
  
    % this figures out the size of the screen and puts the figure
    %    at the left side of the screen, one fourth of the way up
    %    the screen, with width 1/4 of the screen width and height
    %    1/2 of the screen height.
    scrsz = get(0,'ScreenSize');
    figh=figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    
    txtbox1=annotation(figh,'textbox',...
      [0.629 0.9 0.358 0.0792],...
      'String',{'Removed from system'},...
      'FontSize',16,...
      'FitBoxToText','on');

    txtbox2=annotation(figh,'textbox',... 
      [0.15 0.51 0.332 0.0792],...
      'String',{'Added to system'},...
      'FontSize',16,...
      'FitBoxToText','on');

    % contColor is a cell array giving the colors of each ball
    %    (unit volume) in the box
    contColor={};
    % initially we have initWhite white balls and numBlack black
    for j=1:initWhite
        contColor{end+1}='w';
    end
    for j=1:numBlack
        contColor{end+1}='k';
    end
    
    % plot the initial state
    subplot(2,1,1), plotContainer(contColor);
    title_handle=title('Balls in a Box');
    set( title_handle, 'FontSize', 18, 'Position', [0,boxHeight+1] );
    subplot(2,1,2), plot(1,100*numBlack/length(contColor),'k-',...
        t(1:1),100*theorySol(1:1)/volume,'b-','LineWidth',2);
    legend_handle = legend( '% of black balls', ...
        'solution to ODE', 'Location', 'SouthEast');
    set( legend_handle, 'FontSize', 14 );
    
     % simulation loop
    for j=2:simulationLength
        % shake the box (randomly relocate all of the dots)
        contColor=contColor(randperm(length(contColor)));
        % then initially we have the same number of black balls
        %    as we did at the last time interval
        numBlack(j)=numBlack(j-1);
        % the loop below adds blackInPerFrame black balls, bumping
        %    a dot off the end of the array for each one we add and
        %    keeping track of the number of black balls in the
        %    box
        for k=1:blackInPerFrame
            contColor=cat(2,'k',contColor);
            cout=contColor{end};
            if(cout~='k')
                numBlack(j)=numBlack(j)+1;
            end
            contColor=contColor(1:end-1);
        end
        % similar for the number of white balls added
        for k=1:whiteInPerFrame
            contColor=cat(2,'w',contColor);
            cout=contColor{end};
            if(cout=='k')
                numBlack(j)=numBlack(j)-1;
            end
            contColor=contColor(1:end-1);
        end
    
        if(j/frameskip==floor(j/frameskip))
            subplot(2,1,1),...
              plotContainer(contColor);
            title_handle=title('Balls in a Box');
            set( title_handle, 'FontSize', 18, ...
                'Position', [0,boxHeight+1] );
            
            subplot(2,1,2),...
              plot(1:length(numBlack),100*numBlack/length(contColor),...
              'k-',t(1:j),100*theorySol(1:j)/volume,'b-','LineWidth',2);
            legend_handle = legend( '% of black balls', ...
                'solution to ODE', 'Location', 'SouthEast');
            set( legend_handle, 'FontSize', 14 );
            
            axis([1 simulationLength 0 volume]);
            
            getframe();
        end
    end
  
    % this function plots the box full of white and black dots
    function plotContainer(contColor)
        plot(1,boxWidth,'w.');
        inboxULX=.5; inboxULY=1.5;
        inboxLLX=.5; inboxLLY=.5;
        inboxURX=blackInPerFrame+whiteInPerFrame+.5;inboxURY=1.5;
        inboxLRX=blackInPerFrame+whiteInPerFrame+.5;inboxLRY=.5;
        inboxX=[inboxULX inboxURX inboxLRX inboxLLX inboxULX];
        inboxY=[inboxULY inboxURY inboxLRY inboxLLY inboxULY];
        line(inboxX,inboxY,'Color',[.25 .25 .25],'LineWidth',2,...
            'LineStyle',':')
        axis([.3 (boxWidth+0.7) .3 (boxHeight+0.7)])
        outboxULX=boxWidth-blackInPerFrame-whiteInPerFrame+.5;
        outboxULY=boxHeight+.5;
        outboxLLX=boxWidth-blackInPerFrame-whiteInPerFrame+.5;
        outboxLLY=boxHeight-.5;
        outboxURX=boxWidth+.5; outboxURY=boxHeight+.5;
        outboxLRX=boxWidth+.5; outboxLRY=boxHeight-.5;
        outboxX=[outboxULX outboxURX outboxLRX outboxLLX outboxULX];
        outboxY=[outboxULY outboxURY outboxLRY outboxLLY outboxULY];
        line(outboxX,outboxY,'Color',[.25 .25 .25],'LineWidth',2,...
            'LineStyle',':')
        set(gca,'XTick',[]);
        set(gca,'YTick',[]);
        dummyA=zeros(boxWidth,boxHeight);
        hold on;
        for k=1:boxWidth
            for ell=1:boxHeight
                index=sub2ind(size(dummyA),k,ell);
                plot(k,ell,[contColor{index} '.'],'MarkerSize',24);
            end
        end
        hold off;
    end

    function yprime=mixingODE(t,y)
       yprime = blackInPerFrame - holeSize*y/volume;
    end
end

% all done
%----------------------------------------------------------------------
