% first pass at a Monte Carlo simulation of deposition and movement of
% atoms on a surface.

% notation:  Here's an example of atom locations:
%          
% 3        X   X
% 2      X   X   X   X
% 1    X   X   X   X   X  X
%
%      1 2 3 4 5 6 7 8 9
%
% top(1)=1, top(2)=2, top(3)=3, top(4)=2, top(5)=3, top(6)=2, top(7)=1

% other m-files needed:
% plotatom      - plots a single atom
% plotallatoms  - plots all atoms
% moveatom      - plots movement of an atom
%
% Only certain motions are allowed in this version:  only isolated atom can
% move, for example

% probabilities of action at each location, each step:
pmove = .1;       % probability of movement in each direction (must be < 0.5!)
pdeposit = .005;   % probability of deposition

% region to plot:
i1 = 80;
i2 = 120;
% we work over a larger region 60:140.

% initialize even larger region so we don't need to worry about boundary 
% conditions 
% Two layers to start:
top = zeros(1,200);
top(1:2:200) = 1;
top(2:2:200) = 2;

% a little hump to start:
top(101) = 3;
top(102) = 4;
top(103) = 5;
top(104) = 4;
top(105) = 3;


% plot initial layers:
plotallatoms
pause(.5)


for step = 1:200

   % deposition of new atoms:
   for i=60:140
      if (top(i-1)>= top(i) & top(i+1)>=top(i) & rand(1)<pdeposit)
	 top(i) = top(i) + 2;
	 moveatom(i,top(i)+6,i,top(i));
	 plotallatoms
	 end
      end % loop on i for deposition

  % motion of existing atoms:
  for i=60:140
  if (top(i)>top(i-2) & top(i)>top(i+2))
      % single atom, could possible move
      random_num = rand(1);
      if (random_num < pmove)
	  % jump left
	  topi = top(i);
          top(i) = top(i) - 2;
	  top(i-2) = top(i-2) + 2;
	  moveatom(i,topi,i-2,top(i-2));
	  plotallatoms

	  j = i-2;
	  while (top(j-1) < top(j+1))
	     % if its not on a flat surface, move left again:
	     topj = top(j);
	     top(j) = top(j)-2;
	     top(j-1) = top(j-1) + 2;
	     moveatom(j,topj,j-1,top(j-1));
	     plotallatoms
	     j = j-1;
	     end

      else if (random_num > 1-pmove)
	  % jump right
	  topi = top(i);
          top(i) = top(i) - 2;
	  top(i+2) = top(i+2) + 2;
	  moveatom(i,topi,i+2,top(i+2));
	  plotallatoms

	  j = i+2;
	  while (top(j+1) < top(j-1))
	     % if its not on a flat surface, move right again:
	     topj = top(j);
	     top(j) = top(j)-2;
	     top(j+1) = top(j+1) + 2;
	     moveatom(j,topj,j+1,top(j+1));
	     plotallatoms
	     j = j+1;
	     end
	  end %jump right
      end
      end
    end % loop on i

  end % loop on step