function [qrefl,qtrans,qcav,prefl,ptrans,pcav]=cav(qin,l,m,lambda,L,R1,R2,r1,r2,l1,l2,n_iter,varargin)
% Propagates a Gaussian beam through a cavity. 
% The function returns the q-factors of the reflected beam, the transmitted 
% beam and the cavity beam at the first mirror. Also returns the reflected, 
% transmitted and intracavity field amplitudes. (Notation: prefl,ptrans,pcav
% refer to the complex field amplitudes from each bounce, not power. Thus the 
% total transmitted power, for example, is Ptransmitted = epsilon0*abs(sum(ptrans))^2.
%
% AUTHOR: Andri M. Gretarsson, 2003.
% LAST MODIFIED: 2007 by AMG.
%
% SYNTAX: [qrefl,qtrans,qcav,prefl,ptrans,pcav]=...
%         cav(qin,l,m,lambda,L,R1,R2,r1,r2,l1,l2,n_iter <,n,Lmirr,pin>);
%
% INPUT VARIABLES
% ---------------
% R1,R2     = Radii of curvature of the end mirrors.  Positive if mirror is 
%             concave as seen from inside the cavity, negative otherwise. 
%             (Sides facing outwards are assumed to be flat.) R1 corresponds
%             to the input mirror and R2 to the end mirror.
%
%           FIGURE:
%           
%           Input beam   Input mirror     End Mirror    Transmitted beam
%           -----------------|(==============)|- - - - - - - - - - - 
%                                Cavity beam
%
%           Note: Reflected beam is in same location as the Input beam but 
%           travels in the opposite direction (away from the cavity).
%     
% L         = Cavity length
% l,m       = TEM_lm Gaussian mode incidnet on the cavity.
% lambda    = wavelength
% qin       = q of the incoming beam immediately before the input optic 
%             outward facing side (assumed flat).
% r1        = Reflectoin coefficient (fraction of field _amplitude_ reflected) of input mirror.
% r2        = Reflection coefficient of end mirror.
% l1        = Loss coefficient for a single pass through the input mirror.
%             (Fraction of field _amplitude_ lost.)
% l2        = Loss coefficient for a single pass through the end mirror.
%             (Fraction of field _amplitude_ lost.)
% n_iter    = Number of iterations (cavity traversals) to calculate.
% n         = Optional: 1x3 vector of indices of ref. of: mirror substrates,
%             cavity medium, external medium. 
%             Default is n=[1.46,1,1].
% Lmirr     = Optional: 1x2 vector of mirror thicknesses (on the optic
%             axis). First value corresponds to the input mirror, second to
%             the end mirror. Default: Lmirr=[0,0].
% pin       = Amplitude factor of the beam entering the cavity
%
% NOTES:  In the current version, antireflective coatings on the outside
% faces (flat faces) of the mirrors are assumed to be perfect.  This could
% be improved in future versions. Also, the reflectivity and transmissivity
% of the coatings are assumed to be the same for light incident from
% either the substrate side or cavity side of the coating.  Clearly, this
% could be improved also.
%
%--------------------------------------------------------------------------
% SYNTAX: [qrefl,qtrans,qcav,prefl,ptrans,pcav]=...
%         cav(qin,l,m,lambda,L,R1,R2,r1,r2,l1,l2,n_iter <,n,Lmirr,pin>);
%--------------------------------------------------------------------------

if ( nargin>=13 && ~isempty(varargin{1}) ), n=varargin{1}; else n=[1.46,1,1]; end
if ( nargin>=14 && ~isempty(varargin{2}) ), Lmirr=varargin{2}; else Lmirr=[0,0]; end
if ( nargin>=15 && ~isempty(varargin{3}) ), pin=varargin{3}; else pin=1; end

nsubs=n(1);
ncav=n(2);
noutside=n(3);
L1=Lmirr(1);
L2=Lmirr(2);

% Variable initialization
t1 = sqrt(1-r1^2-l1^2);
t2 = sqrt(1-r2^2-l2^2);
qrefl = zeros(n_iter,1); qcav = zeros(n_iter,1); qtrans = zeros(n_iter,1);
prefl = zeros(n_iter,1); pcav = zeros(n_iter,1); ptrans = zeros(n_iter,1);
% #cav  ->  immediately inside input mirror travelling towards end mirror
% #refl ->  immediately outside input mirror (outside cavity) travelling towards laser
% #trans->  immediately outside end mirror (outside cavity) travelling away from cavity

% Note the poor notation for field amplitudes throughout:  prefl, pcav,
% etc.  The letter p does NOT mean power in this case!

% Commonly used quantitites
freeL  = free(L,ncav);
freeL1 = free(L1,nsubs);
freeL2 = free(L2,nsubs);
mirrR1 = mirr(R1);              % input mirror
mirrR2 = mirr(R2);              % end mirror
% mirrR1phase=exp(i*2*pi*(nsubs*L1/lambda-floor(nsubs*L1/lambda)));
% mirrR2phase=exp(i*2*pi*(nsubs*L2/lambda-floor(nsubs*L2/lambda)));
% tripphase=exp(i*2*pi*(ncav*L/lambda-floor(ncav*L/lambda)));
mirrR1phase=exp(i*2*pi*(nsubs*L1/lambda));
mirrR2phase=exp(i*2*pi*(nsubs*L2/lambda));
tripphase=exp(i*2*pi*(ncav*L/lambda));

% Prompt reflection
[qrefl(1),prefl(1)]= prop(qin,...
    fdie(nsubs,noutside)*freeL1*mirr(-R1)*freeL1*fdie(noutside,nsubs),...
    [l,m],-r1*mirrR1phase^2*pin);

qcav(1)=i*1; pcav(1)=0;
qtrans(1)=i*1; pcav(1)=0;

% First entry into cavity
[qcav(2),pcav(2)] = prop(qin,...
    sdie(-R1,nsubs,ncav)*freeL1*fdie(noutside,nsubs),...
    [l,m],t1*mirrR1phase*pin);

% First transmission through end mirror
[qtrans(2),ptrans(2)] = prop(qcav(2),...
    fdie(nsubs,noutside)*freeL2*sdie(R2,ncav,nsubs)*freeL,...
    [l,m],t2*tripphase*mirrR2phase*pcav(2));

% First round trip leakage through input mirror
[qrefl(2),prefl(2)] = prop(qcav(2),...
    fdie(nsubs,noutside)*freeL1*sdie(R1,ncav,nsubs)*freeL*mirrR2*freeL,...
    [l,m],t1*r2*mirrR1phase*tripphase^2*pcav(2));

for s = 3:n_iter+2
    
    [qcav(s),pcav(s)] = prop(qcav(s-1),...
        mirrR1*freeL*mirrR2*freeL,...
        [l,m],r1*r2*tripphase^2*pcav(s-1));

    [qtrans(s),ptrans(s)] = prop(qcav(s),...
        fdie(nsubs,noutside)*freeL2*sdie(R2,ncav,nsubs)*freeL,...
        [l,m], t2*tripphase*mirrR2phase*pcav(s));
      
    [qrefl(s),prefl(s)] = prop(qcav(s),...
        fdie(nsubs,noutside)*freeL1*sdie(R1,ncav,nsubs)*freeL*mirrR2*freeL,...
        [l,m],t1*r2*mirrR1phase*tripphase^2*pcav(s));
    

end