% scripts for Dec 5, 2017   MAT i331 


% scripts for tuesday Nov 28 2017

%------------------------------------------



% this is a function to compute the harmonic measure of one side
% of a rectable from its center. We assume the box is of the 
% form [-a,a] times [-1,1] and the grid size in 1/m, where a and m
% are inputs.  A random walk is run N times and we count the number 
% of times it first hits the boundary on the right hand side.

function out=hm_box(a,m,N)
   tic;
   delta=1/m;
   ep=.000001;
   total=0;      
   for n=1:N
       flag=0;
       x=0;
       y=0;
       points=[x,y];
       while flag==0;
          step=randi(4);
          switch step
              case 1
                  x=x+delta;
              case 2
                  x=x-delta;
              case 3 
                 close all
                 y=y+delta;
              case 4
                  y=y-delta;
          end % switch
          points=[points;x,y];
          if  (abs(x)>a-ep)||(abs(y)>=1-ep)
              if x>a-ep
                  total=total+1;
              end 
              save(n,:)=[x,y];
              flag=1;
          end   
      end % while
   end % for n
   prob=total/N
   toc
   figure;
   axis equal;
   hold on
   x=[-a,a,a,-a,-a];
   y=[-1,-1,1,1,-1];
   plot(x,y,'LineWidth',2)
   plot(points(:,1),points(:,2),'b');
   scatter(save(:,1),save(:,2),'r','filled')
return

%----------------------------------------

% this is a function to compute the harmonic measure of one side
% of a rectable from its center. We assume the box is of the 
% form [-a,a] times [-1,1] and the grid size in 1/m, where a and m
% are inputs.  A kakutani walk is run N times and we count the number 
% of times it first hits the boundary on the right hand side.

function out=kakutani(a,N)
   tic;
   ep=.000001;
   total=0;      
   for n=1:N
       flag=0;
       x=0;
       y=0;
       points=[x,y];
       while flag==0;
          radius=min(a-abs(x),1-abs(y));
          z=radius*exp(2*pi*1i*rand);
          x=x+real(z);
          y=y+imag(z);
          points=[points;x,y];
          if  (abs(x)>a-ep)||(abs(y)>=1-ep)
              if x>a-ep
                  total=total+1;
              end 
              save(n,:)=[x,y];
              flag=1;
          end   
      end % while
   end % for n
   prob=total/N
   toc
   figure;
   axis equal;
   hold on
   x=[-a,a,a,-a,-a];
   y=[-1,-1,1,1,-1];
   plot(x,y,'LineWidth',2)
   plot(points(:,1),points(:,2),'b');
   scatter(save(:,1),save(:,2),'r','filled')
return


%-----------------------------------------------


% Draw diffusion limited aggregation

function   out=dla(n)
   clear centers;
   clear lograd;
   tic;
   centers=[0];
   counter=1;
   attempts=1;
   %centers=zeros(1,n);
   hold on;
   while counter < n
      maxradius=max(abs(centers));
      z=(maxradius+2)* exp(2 * pi * i * random('unif',0,1));
      step=min(abs(z-centers))-2;
      while((step>.01)&&(step<5000))
         z=z+ step * exp( 2 * pi * i * random('unif',0,1));
         step=min(abs(z-centers))-2;
      end
      if (step<.2)
          attempts=attempts+1;
              counter=counter+1;
              centers(counter)=z;
      end % if
      if counter==1000
         figure;
         axis equal;
         hold on;
      end
      if mod(counter,1000)==0
        title(['DLA, n = ',num2str(counter)])
       %  title('DLA');
        scatter(real(centers),imag(centers) ,'k','.');
        drawnow;
        [counter/1000,toc]
      end % if
   end % for k
   if  n< 1000
       figure;
       hold on; grid on; axis equal;
       plotcirc(centers,zeros(1,length(centers)));
   else
      figure
       hold on;
       axis equal;
       title(['DLA, n = ',num2str(length(centers))])
       colormap jet;
       scatter(real(centers),imag(centers),[],[1:length(centers)],'.');
   end
   
   return 
   figure;   hold on;   grid on;
   title('Cummulative maximum of DLA centers')
   cm=cummax(abs(centers));
   plot(cm);
   figure;   hold on;   grid on;
   title('Log of cummax of DLA centers')
   x=[1:length(cm)];
   plot( log(cm(100:end))./log(x(100:end)) );
  toc;
  
  return


%--------------------------------------------------------

% Draw DLA, but keep tract of the vertices of the convex 
% hull at each stage. Draw DLA with two colors, depending 
% on whether a new disk is on the convex hull boundary when 
% it is added.

function   dla2(n)
   tic
   clear centers;
   clear lograd;
   tic;
   centers=[0,2];
   hullx=[0,2];hully=[0,0];
   pointer=[0,0;1,2];
   counter=2;
   attempts=2;
 %  centers=zeros(1,10^7);
   for k=1:n
      maxradius=max(abs(centers))+2;
      z=(maxradius+1)* exp(2 * pi * i * random('unif',0,1));
      step=min(abs(z-centers))-2;
      while((step>.01)&&(step<5000))
         z=z+ step * exp( 2 * pi * i * random('unif',0,1));
         step=min(abs(z-centers))-2;
      end % while
      if  (step > 4999)  % use Poissonkernel to retart 
      end % if 
      if (step<.02)
          [m,index]=min(abs(z-centers));
          pointer=[pointer;index,z];
          attempts=attempts+1;
          counter=counter+1;
          centers(counter)=z;
          newhullx=[hullx,real(z)];
          newhully=[hully,imag(z)];
          K=convhull(newhullx,newhully);
          hullx=newhullx(K);
          hully=newhully(K);
          saveK(counter)=length(K);
          if max(K)==length(newhullx)
             hull(counter)=1;
          else 
             hull(counter)=0;
          end % if
          hull;          
      end % if 
      if mod(k,10000)==0
         [k/10000,toc]
      end % if
   end % for k
   
  figure;  hold on; axis equal;
  title('Red disks on convex hull when added');
  list=find(hull==0);
  L1=length(list);
  scatter(real(centers(list)),imag(centers(list)),'.','k');
  list=find(hull==1);
  L2=length(list);
  percentred=L2/(L1+L2)
  scatter(real(centers(list)),imag(centers(list)),'.','r');
  % plot dla 
  plotcirc(centers,hull);
  figure; hold on; grid on;
  title('Number of points on convex hull at time t');
  plot(saveK)
  %figure
  %smooth=conv(saveK,ones(1,10))/10;
  %smooth=smooth(10:10:length(saveK-10));
  %plot(smooth);
  toc

return




%---------------------------------------

% function to plot unit circles in the plane.
% this is called by dla2 to draw the dla disks
% for small aggregates.  The disks are colored
% red or black, depending on the value in 'list'

function plotcirc(z,list)

    figure;
    hold on;
    axis equal;
    t=[0:.1:2.1*pi];
    rad=1;
    for k=1:length(z)
       x=real(z(k));
       y=imag(z(k));
       if list(k)==0
          plot([x+rad*cos(t)],[y+rad*sin(t)],'k');
       else
          plot([x+rad*cos(t)],[y+rad*sin(t)],'r');
       end
    end % for k

return;

%------------------------------------------------
%-----------------------------------------------



function paths5(N)
tic;
% set number of steps in the random walk
%N=1000000;
t=randi(4,1,N);
u=round(cos(.5*pi*t));
v=round(sin(.5*pi*t));
% define the random walk
x=cumsum(u)-u(1);
y=cumsum(v)-v(1);

% find min and max value that path reaches
minx=min(x);
maxx=max(x);
miny=min(y);
maxy=max(y);
% compute size of smallest grid containing path
xrange=maxx-minx;
yrange=maxy-miny;
% define an array big enough to hold path
% plus a empty `frame' around path  
visited=zeros(xrange+4,yrange+4);

% redefine x and y vectors to lie in grid 
% but not touch outer most sites 
newx=x-minx+2;
newy=y-miny+2;

for k=1:N
     visited(newx(k),newy(k))=1;
end % for k 


% define distance to starting point
% if distance is already define do nothing, 
% otherwise distance is one more that 
% the minimum of four surrounding sites
% maxpath keeps track of how far from 
% origin we can get in k steps

grid=zeros(xrange+4,yrange+4);
checked=zeros(xrange+4,yrange+4);
upper=1000000;
maxpath=0;
grid=grid+upper;
x1=newx(1);
y1=newy(1);
grid(x1,y1)=0;
checked(x1,y1)=1;
% listx and listy contain coordinates of set of points 
% distance k from starting point 
clear listx;
clear listy; 
clear newlistx;
clear newlisty;
newlistx(1)=x1;
newlisty(1)=y1;
listlength=1;
flag=1;
dist=0;
while flag==1
   dist=dist+1;
   counter=0;
   listlength=length(newlistx);
   listx=newlistx;
   listy=newlisty;
   clear newlistx;
   clear newlisty;
   for k=1:listlength
       x1=listx(k);
       y1=listy(k);
       % check each neighbor if in trace and not checked
       nextx=x1+1;
       nexty=y1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1; 
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if 
       nextx=x1-1; nexty=y1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1; 
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if 
       nextx=x1; nexty=y1+1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1; 
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if 
       nextx=x1; nexty=y1-1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1; 
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if 
   end % for k
   if counter==0
      flag=0;
   end % if 
end % while 


maxpath=0;
for k=1:N
    if grid(newx(k),newy(k)) > maxpath
       maxpath=grid(newx(k),newy(k));
       maxx=newx(k);
       maxy=newy(k);
    end  % if grid 
end % for k

% find path between origin and maximal distance path

   clear savex;
   clear savey;
   oldx=maxx;
   oldy=maxy;
   savex(1)=oldx;
   savey(1)=oldy;
   for k=1:maxpath
      dist=maxpath-k;
      if   grid(oldx-1,oldy)==dist
           nextx=oldx-1;
           nexty=oldy; 
      end % if
      if   grid(oldx+1,oldy)==dist
           nextx=oldx+1;
           nexty=oldy; 
      end % if
      if   grid(oldx,oldy-1)==dist
           nextx=oldx;
           nexty=oldy-1; 
      end % if
      if   grid(oldx,oldy+1)==dist
           nextx=oldx;
           nexty=oldy+1; 
      end % if
      savex(k+1)=nextx;
      savey(k+1)=nexty;
      oldx=nextx;
      oldy=nexty; 
   end % for k 
   figure;
   axis equal;
   hold on; 
   set(gca,'visible','off');
   plot(newx,newy);

   figure;
   axis equal;
   hold on; 
   set(gca,'visible','off');
   plot(newx,newy);
   plot(savex,savey,'k','LineWidth',3);
   toc;
return