%--------------------------------------------------------------------------
% FTE_zernike_scattering.m
%
% An example Matlab script which plots the coupling from one LG mode into
% another according to an approximation using Zernike polynomials.  The
% example calculates the coupling from a surface made up of several Zernike
% polynomials with random amplitudes  and uses Finesse Tools to approximate
% the coupling from one mode to another via each individual polynomial.  
% The approximate coupling is then plotted as we include more polynomial 
% orders (n) in our approximation.  This value is compared with the 
% coupling value generated numerically from the complete map.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    03.07.2012
%--------------------------------------------------------------------------

clear all;

% Beam parameters
w = 3.2e-2;             
lambda = 1064e-9;
k = 2*pi/lambda;
gp = FT_init_gauss_param(lambda,1,'z',0,'w0',w);

% Mirror parameters
R = 17e-2;
gridsize = 512;

% Incident and coupled modes
p = 3;
l = 3;
p_ = 2;
l_ = 5;

% Number of terms to be used in approximation.  N.B. can only use 2 terms
% (i.e. exp(2ikZ) = 1 + 2ikZ is expanded up to first order) if we are
% adding the coupling coefficients to give coupling from entire map, the 2
% order terms can not simple be added together (see
% FT_zernike_coupling_of_LG_modes.m)
no_terms = 2;

% 3 Zernike polynomials to add to make up surface.  All couple
% significantly from p,l into p_,l_
n = [2 4 6 8 10 12 14 16 18 20];
m = 2*(2*randi(2,1,10)-3);
nmax = 20;
% Amplitude in nm.  The gaussian weighting is used to produce a more
% 'realistic' distribution of 'random' amplitudes for different order
% polynomials, i.e. the higher order polynomials (higher spatial
% frequencies) have greater chance of having a lower amplitude.
amps = 4*rand(1,10).*exp(-((1:10)/3).^2);

figure
    plot(n,amps,'b','LineWidth',2)
    xlabel('n')
    ylabel('amplitude')
    %set(gca,'XScale','Log')
    set(gca,'XLim',[min(n) max(n)])
    set(gca,'YScale','Log')
    grid on;
    drawnow;

% Add polynomials to zernike coefficients structure
zc = FT_init_zernike_coefficients(nmax);
for i=1:length(n)
    zc = FT_update_zernike_coefficient(zc,n(i),m(i),amps(i));
end

% Create map from zernike coefficients
zmap = FT_create_synthetic_zernike_map(zc,'original',2*R,gridsize,R,'ztest');

% Plot map
    FT_plot_mirror_map(zmap,0)
    drawnow;
    
[rows,cols] = size(zmap.data);

% x and y to match map with the LG mode
x = zmap.xstep * linspace(zmap.x0-cols,zmap.x0,cols);
y = zmap.ystep * linspace(zmap.y0-rows,zmap.y0,rows);

% Set area outside the map as 0
zmap.data(zmap.nan)=0;
Z = zmap.scaling.*zmap.data;

% Reflected field and coupled mode
field1 = FT_LG_field(gp,p,l,x,y,[0,0,0]) .* exp(2*1i*k*Z);
field2 = FT_LG_field(gp,p_,l_,x,y,[0,0,0]);

% Plot field
figure
    FT_FFT_plot_field(field1,1)
    drawnow;

nn=0:nmax;

% If we're looking at the coupling back into the incident mode look at the
% coupling out of this mode for simplicity
if p==p_ && l==l_
    kplpl = 1;
else
    kplpl = 0;
end

% Coupling from complete map
k_map(1:length(nn))=FT_conv_fields(field2,field1) - kplpl;

kk = 0;

% Coupling from individual Zernike polynomials
for n=0:nmax    
    
    for m=-n:2:n
        
        A = FT_get_zernike_coefficient(zc,n,m);
        cc = FT_zernike_coupling_of_LG_modes(p,l,p_,l_,lambda,w,R,A,n,m,no_terms);
        
        kk = kk + (cc-kplpl);
    end
    
    k_zernikes(n+1) = kk;
end

% Plot coupling 
figure
    hold on;
    plot(nn,abs(k_zernikes),'b','LineWidth',2)
    plot(nn,abs(k_map),'r--','LineWidth',2)
    % Labels
    title('Zernike approximate coupling')
    xlabel('Zernike order')
    ylabel('Coupling coefficient')
    legend('k from map','k from sum of individual Zernike ks','Location','SouthEast')
    
        
        