%--------------------------------------------------------------------------
% FTE_zernike_coupling_test.m
% 
% An example script for Matlab which plots the results for coupling between
% two LG modes incident on a surface described by a Zernike polynomial.
% This example compares the coupling calculated via an analytical
% approximation (FT_zernike_coupling_of_LG_modes.m) and calculated 
% numerically.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    18.05.2012
%--------------------------------------------------------------------------
%

clear all;

% aligo cavity mirror radius
R = 17e-2;      

% 1nm amplitude
A = 1;

% Input mode
lambda = 1064e-9;
k = 2*pi/lambda;
p = 3;
l = 3;

% Coupled mode
p_ = 2;
l_ = 5;

% Zernike polynomial with required m (astigmatic polynomial Z2,2)
n = 2;
m = 2;

% Calculate coupling over range of beam sizes
w = linspace(0.0001,R,100);

% Number of grid points and optical centre
npoints = 1000;
x0 = (npoints+1)/2;
y0 = (npoints+1)/2;

% Timing variables
total_runs = length(w);
elapsed_time = 0;

for i=1:length(w)
    
    tic;
    
    % Grid coordinates for numerical calculation
    r = linspace(-12*w(i),12*w(i),npoints);
    x = r;
    y = r;
    
    % Gaussian parameter
    gp = FT_init_gauss_param(lambda,1,'z',0,'w0',w(i));
    
    % Zernike radius in grid points 
    Rz = R/abs(r(1)-r(2));
    
    % Calculate zernike surface
    Z = A * 1e-9 * FT_zernike_polynomial_for_grid(npoints,npoints,x0,y0,Rz,Rz,n,m);    
        
    % Make points outside polynomial 0
    idx = find(isnan(Z)==1);
    Z(idx) = 0;
    
    % Calculate reflected field
    field1 = FT_LG_field(gp,p,l,x,y,[0,0,0]) .* exp(2*1i*k*Z);
    
    % Coupled field
    field2 = FT_LG_field(gp,p_,l_,x,y,[0,0,0]);
    
    % Calculate numerical and analytical coupling
    k_num(i) = FT_conv_fields(field2,field1);
    k_analytic(i) = FT_zernike_coupling_of_LG_modes(p,l,p_,l_,lambda,w(i),R,A,n,m,3);

    elapsed_time = FT_print_progress(elapsed_time,i,total_runs,0);
    
end

P_num = k_num.*conj(k_num);
P_analytic = k_analytic.*conj(k_analytic);

% Plot coupling results
figure
    % Amplitude 
    hold on;
    plot(w/R,abs(k_num),'b','LineWidth',3)
    
    plot(w/R,abs(k_analytic),'r--','LineWidth',3)
   
    legend('numerical','analytic')
    
    % Beam size compared to mirror/zernike radius
    xlabel('Beam radius [R] ')
    ylabel('Coupling coefficient ')
    

    
    
    
    
    
    
    
