eta=-1; k=0.01;


ro=i*eta/2/(1-i*eta/2);
to=1/(1-i*eta/2);

for in=0:99
  k=2*pi/100*in;
  p=exp(i*k);
  Rl=ro;
  Rr=ro;
  Tl=to;
  Tr=to;
  dat=[];
  for j=1:400
    t=to+(rand(1)+i*rand(1))/10.;
    while (t*t' > 1)
      t=to+(rand(1)+i*rand(1))/10.;
    end
    r=-sqrt(-t/t')*sqrt(1-t'*t);

    Tln=Tl*p*t/(1-p^2*r*Rr); 
    Trn=t*p*Tr/(1-p^2*r*Rr);

    Rln=Rl+p^2*Tl*Tr*r/(1-p^2*r*Rr);
    Rrn=r+p^2*t*t*Rr/(1-p^2*r*Rr);

  
    Rl=Rln; Rr=Rrn; Tl=Tln; Tr=Trn;
%    P=[abs(Rl)^2+abs(Tl)^2 abs(Rr)^2+abs(Tr)^2]
    dat=[dat ; Tl];
  end
  abs(Tl)
  hold off
  plot(dat,'b-')
  axis([-1 1 -1 1])
  hold on
  plot(real(dat(1)),imag(dat(1)),'r*')
  in
  pause
%  if (in<10)
%    eval(['print -dgif8 out0' num2str(in) '.gif']);
%  else
%    eval(['print -dgif8 out' num2str(in) '.gif']);
%  end
end