====== Matlab script ======

% caculates the cavity spectrum and the frequency distance between TEM_00
% and higher order modes m+n. Two curved mirrors possible.
%
% 1. Vers.: 14.03.12 , T.Legero
% Revised : 19.11.13 , T. Legero , Two curved mirrors possible.

clear all;
clc;

Len = 0.48;       % cavity length in [m]
Curv1 = 1000000000;        % mirror no. 1 curvature in [m]
Curv2 = 1;        % mirror no. 2 curvature in [m]
Refl = 0.99998;   % reflectivity of mirrors
c = 299792458;    % speed of light [m/s]

FSR = c /(2*Len);   % Free-Spectral-Range
Gouy = FSR/pi*acos(sqrt((1-Len/Curv1)*(1-Len/Curv2)));  % Gouy-Phase

% Finesse = pi*sqrt(Refl)/(1-Refl);
% linewidth = FSR/Finesse;

fprintf('Cavity length       L = %g m  =>  FSR = %g MHz\n', Len, FSR/1E6);
fprintf('Radius of Curvature R1 = %g m and R2 = %g m => Gouy = %g MHz\n\n', Curv1, Curv2, Gouy/1E6);

%fprintf('Reflectivity R = %g  => Finesse  F = %d \n => Linewidth  d = %g MHz \n\n' , Refl, round(Finesse),linewidth/1E6);

mn_max = 100;
w_max = 100;

delta = zeros(mn_max,w_max);

for w = 1:w_max
    for mn = 1:mn_max
        delta(mn,w) = w*FSR - mn*Gouy; % Freq.Diff. von (q-w)ter höherer 
                                       % Mode mn zu q-ter TEM_00-Mode 
        if (abs(delta(mn,w)) < 10E6)
           fprintf('Freq. diff. of m+n = %g higher order mode to 00-mode is %g MHz\n', mn, delta(mn,w)/1E6); 
        end
                                           
    end
end

%% plot

mn = 4;                  % Modenumber of Higher-Order-Mode to plot
R = 0.95:0.01:1.05;      % Vary Mirror-Curvature

G = FSR./pi*acos(sqrt(1-Len./R));  % Gouy-Phase
%G = FSR./pi*acos(sqrt((1-Len./R).*(1-Len./R)));  % Gouy-Phase
d = FSR - mn*G;                    % Difference between TEM_00 and TEM_mn

plot(R,G);