function [q,p]=prop(q1,abcd,varargin)
% Propagates a Gaussian beam. [q,p]=prop(q1,abcd <,[l,m],p1>);
%--------------------------------------------------------------------
% Propagates a Gaussian beam with complex radius of curvature q1
% and amplitude factor p1 (optional), according to the ABCD matrix
% supplied.
%
% Returns the new complex beam radius q=(A*q1+B)/(C*q1+D) and the 
% new amplitude factor p = 1/(A+B/q1)^(1+l+m) by which the field is
% multiplied. (See eqn 4.7.30 in Principles of Lasers by O. Svelto.)
% If q1 is a vector q and p will be vectors of the same size.
%
% For a Hermite Gaussian l,m are the mode designators.
% For a Laguerre Gaussian l=2p and m is the usual m.
%
% SYNTAX: [q,p]=prop(q1,abcd <,[l,m],p1>);
%           <...> indicates optional arguments
%--------------------------------------------------------------------

if ( nargin>=3 && ~isempty(varargin{1}) ), mode=varargin{1}; else mode=[0,0]; end
if nargin>=4, p1=varargin{2}; else p1=ones(size(q1)); end

A=abcd(1,1);
B=abcd(1,2);
C=abcd(2,1);
D=abcd(2,2);

ell=mode(1);
emm=mode(2);

q = (A*q1 + B)./(C*q1 + D);
p = p1.*exp(i*angle(1./(A+B./q1).^(1+ell+emm)));