% 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

%---------------------------------------------------------
% This is a function to compute harmonic measure 
% a  a X 1 rectangle using sparse linear algebra. We 
% assume a is a positive integer and the grid size
% is 1/m where m is also a positive integer.

function out = hm_linalg(a,m)
   tic;
   L1=(2*m +1)
   L2=(2*a*m+1)
   L=L1*L2
   M=eye(L);
   Y=zeros(L,1);
   % left side 
   for k=1:L1
     Y(k)=1;
   end
   % interior
   for k=2:L1-1
     for j=1:L2-2
        M(L1*j+k,L1*(j+1)+(k))=-.25;
        M(L1*j+k,L1*(j-1)+(k))=-.25;
        M(L1*j+k,L1*(j)+(k+1))=-.25;
        M(L1*j+k,L1*(j)+(k-1))=-.25;
     end
   end
   figure
   spy(M)
   X = M\Y;
   prob=X((L+1)/2)
   Z=reshape(X,[L1,L2]);
   depth=50;
   figure
   image(depth*Z);
   axis equal;
       %colormap(flipud(jet(depth)));
       colormap(jet(depth));
       %colormap(hsv(depth));
       %colormap(gray(depth));
       %colormap(flipud(winter(depth)));
       % add a legend saying how colors give iterates
       colorbar;
   toc;

return
pwd