% m30_circles.m - Sinai circular obstacles
%
% Set up circular obstacles and draw them:
  pert = input('perturbation? ');
  hold off, shg, set(gcf,'doublebuffer','on');
  [ii,jj] = meshgrid(-4:4); B = [ii(:)'; jj(:)']; N = length(B);
  radius = 1/3; circle = radius*exp(1i*(0:50)*pi/25);
  figure(1); hold off
  for i = 1:N
    cr = B(1,i) + real(circle); ci = B(2,i) + imag(circle);
    fill(cr,ci,[.5 .5 1]), hold on, plot(cr,ci,'-b','linewidth',.6) 
  end
  axis(1.85*[-1 1 -1 1]), grid on, axis square
% set(gca,'xtick',-2:2,'ytick',-2:2)
  for run = 2:-1:1
    x0 = [.5;.1]; v = [1;0]; t = 0;                              
    if run==2, x0 = x0 + [0;pert]; end
    path = [t [1 1i]*x0];          
    while t < 15
      X0 = repmat(x0,1,N); S = v'*(B-X0)/(v'*v);             
      Y = X0 + v*S; D = sqrt(sum((Y-B).^2));          
      I = find((D<radius) & (S>0));     
      if length(I)==0, break, end
      y = Y(:,I); s = S(:,I); b = B(:,I);
      I = find(s==min(s)); I = I(1);    
      y = y(:,I); s = s(:,I); b = b(:,I);
      x1 = y - sqrt(radius^2-norm(y-b)^2)*v/(v'*v);  
      diam = x1-b; diam = diam/norm(diam);
      v = -v/norm(v); zdiam = diam(1)+1i*diam(2);
      zv = v(1) + 1i*v(2); zv = zdiam*(zdiam/zv);
      v = [real(zv); imag(zv)]; t = t + norm(x1-x0);
      path = [path; [t [1 1i]*x1]]; x0 = x1;
    end
    lastpt = interp1(path(:,1),path(:,2),15,'linear','extrap');
    path(end,:) = [15 lastpt]; tt = 0:.075:15;
    zpath = interp1(path(:,1),path(:,2),tt,'linear','extrap');
    if run==1
      xpath1 = real(zpath); ypath1 = imag(zpath); path1 = path;
    else
      xpath2 = real(zpath); ypath2 = imag(zpath); path2 = path;
    end
  end
  h2 = plot(xpath2(1),ypath2(1),'.r','markersize',28);
  h1 = plot(xpath1(1),ypath1(1),'.k','markersize',28);
  tnext = 0;
  for i = 2:length(xpath1) 
    t = tt(i);
    if t>=tnext
      title(['t = ' int2str(tnext)],'fontsize',20)
      tnext = tnext+1;
    end
    set(h2,'xdata',xpath2(i),'ydata',ypath2(i))
    set(h1,'xdata',xpath1(i),'ydata',ypath1(i))
    drawnow
  end
  plot(path2(:,2),'-r','linewidth',1.2)
  plot(path1(:,2),'-k','linewidth',1.2)
  plot(interp1(path2(:,1),path2(:,2),0:15,'linear','extrap'),'.r','markersize',17)
  plot(interp1(path1(:,1),path1(:,2),0:15,'linear','extrap'),'.k','markersize',17)
  title(['at t=15: (x,y) = (' num2str(real(lastpt)) ...
     ',' num2str(imag(lastpt)) ')'],'fontsize',18)
  disp('press return to see error plot')
  pause
  figure(2)
  err = sqrt((xpath1-xpath2).^2+(ypath1-ypath2).^2);
  semilogy(tt,err,'-','linewidth',2)
  xlabel('t','fontsize',18), ylabel('separation','fontsize',18)
  set(gca,'xtick',0:15), grid on
  title('exponential separation of trajectories','fontsize',20)