%---------------------------------------------------------------------  
% function [q,psi] = FT_ABCD_trace_beam(M,nr, qin)
% 
% Traces a Guassian beam parameter (qin) through a system
% of ABCD matrices (M is a list of ABCD matrices). It returns the
% beam parameters after each matrix and the overall Gouy phase (psi) 
%
% M: vector of ABCD matrices (nx2x2)
% nr: vector of indec of refraction with n(1) being the index
%     of refraction for qin and n(2) for after M(1), ...
% qin: input Guassian parameter for first ABCD matrix M(1,:,:)
%
% q: vector of Gaussian beam parameters for each ABCD matrix
%    e.g. q(3)=M(3)*M(2)*M(1)*qin
% psi: overall Guoy phase
%
% Part of the SimTools package
% Andreas Freise 29.08.08 afreise@googlemail.com
%---------------------------------------------------------------------  

function [q,psi] = trace_beam(M, nr, qin)
  
  [n,x,y]=size(M)
  
  Mtmp=zeros(2,2);
  q=zeros(1,n);
 
  Mtmp=reshape(M(1,:,:),2,2);
  q(1)=qMq(Mtmp,qin,nr(1),nr(2));
  psi=0;
  if (Mtmp(1,2))
    psi=abs(q2Psi(q(1))-q2Psi(qin));
  end
  
  for i=2:n

    Mtmp=reshape(M(i,:,:),2,2);
    q(i)=qMq(Mtmp,q(i-1),nr(i),nr(i+1));
    if (Mtmp(1,2))
      psitmp=abs(q2Psi(q(i))-q2Psi(q(i-1)));
      psi=psi+psitmp;
    end
  
  end
  
  