%--------------------------------------------------------------------------
% function [f,phi] = FT_HOM_cavity_resonance(L,Rc1,Rc2,n,m)
%
% A function for Matlab which calculates the resonace frequency (and mirror 
% tuning) for a Hermite-Gauss higher order mode in a 2 mirror cavity.  The
% resonance is calculated with respect to the 00 mode, i.e. the difference
% in resonance frequency from the 00 mode.  N.B. for astigmatic mirrors the
% resonance frequency will be slightly different for different HGnm modes
% of the same order.
%
% L:    Cavity length [m]
% Rc1:  Radius of curvature of the input mirror [m].  If the mirror is
%       astigmatic give Rc1 as a vector [Rc1x Rc1y] which refer to the 
%       radii of curvature in the x/y direction.  For flat mirrors we have
%       Rc1 = 0;
% Rc2:  Radius of curvature of the end mirror [m]. [Rc2x Rc2y] for
%       astigmatic mirror.  For flat mirror Rc2 = 0;
% n,m:  HG mode indices
%
% f:    Resonance frequency (offset from 00 mode) [Hz]
% phi:  Resonance tuning (or additional accumulated Gouy phase) from 00
%       [deg].
%
% Charlotte Bond    26.11.2012
%--------------------------------------------------------------------------
%

function [f,phi] = FT_HOM_cavity_resonance(L,Rc1,Rc2,n,m)

    baseid = 'HOM_cavity_resonance';
    
    % Error messages
    % Checking number of Rc1/Rc2 values (should be 1 or 2)
    if (length(Rc1)>2 || length(Rc2)>2)
        result='Invalid Rc vectors: Rc1/Rc2 should have length 1 (non-astigmatic) \nor 2 (astigmatic)';
        msgid=[baseid,':checkarguments']
        error(msgid,result);
    end

    % Speed of light
    c = 299792458;
    
    % Set astigmatism in x/y direction depending on values given
    Rc1x = Rc1(1);
    Rc1y = Rc1(length(Rc1));
    Rc2x = Rc2(1);
    Rc2y = Rc2(length(Rc2));
    
    % Initial g-factors (i.e. for flat mirror)
    g1x = 1;
    g1y = 1;
    g2x = 1;
    g2y = 1;
    
    % Set g-factors according to Rcs
    if (Rc1x~=0)
        g1x = 1 - L/Rc1x;
    end
    if (Rc1y~=0)
        g1y = 1 - L/Rc1y;
    end
    if (Rc2x~=0)
        g2x = 1 - L/Rc2x;
    end
    if (Rc2y~=0)
        g2y = 1 - L/Rc2y;
    end

    % Calculate additional Gouy phase to 00 mode
    phi = n*acos(sqrt(g1x*g2x)) + m*acos(sqrt(g1y*g2y));
    
    % calculate resonance frequency
    f = c/(2*L) * (1/pi) * phi;
    
    % Gouy phase in degrees
    phi = phi*180/pi;


end

