%--------------------------------------------------------------------------
% FTE_Winkler_test.m
%
% A Matlab example script which plots the coupling from one mode to another
% due to reflection of a distorted surface.  The coupling is calculated 
% over a range of beam sizes, both analytically (Winkler approximation) and
% numerically.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    18.05.2012
%--------------------------------------------------------------------------

clear all;

% aligo ETM radius
R=17e-2;

% 1nm amplitude distortion
A = 1e-9;
phi = 0;        % Symmetric distortion couples into modes with n = n0+2*l;

% Surface defined by one spatial wavelength
slambda = R/2;

% No tilt of the surface
tilt = 0;

% Mode parameters (coupling from n0 to n)
n0 = 3;
n = 7;
lambda = 1064e-9;
k = 2*pi/lambda;

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

for i=1:length(w)

    % Calculate beam for numerical coupling calculation
    r=linspace(-6*w(i),6*w(i),1000);
    y = 0;
    x = r;
    gp = FT_init_gauss_param(lambda,1,'z',0,'w0',w(i));
    
    % Calculate distorted surface
    z = A*FT_harmonic_surface_distortion(r,slambda,phi);
    
    % Calculate reflected beam
    field0 = FT_HG_field(gp,n0,0,x,y,[0,0,0]);
    field0 = field0.*exp(2*1i*k*z);
    
    % Coupled beam
    field1 = FT_HG_field(gp,n,0,x,y,[0,0,0]);
    
    % Calculate numerical and analytic coefficients (fields should be this
    % way round to give right phase)
    C_num(i) = FT_conv_fields(field1,field0);
    C_winkler(i) = FT_Winkler_coupling(lambda,w(i),n0,n,tilt,A,slambda,phi);
    
end

% Calculate power coupled in ppm
P_num = 1e6 * abs(C_num).^2;
P_winkler = 1e6 * abs(C_winkler).^2;

% Power difference 
dP = abs(P_num - P_winkler);

% Plot result
figure
    hold on;
    plot(w/slambda,P_num,'b','LineWidth',3)
    plot(w/slambda,P_winkler,'r--','LineWidth',3)

    grid on;
    
    % Beam size in multiples of the spatial wavelength
    xlabel('Beam radius [\lambda_s] ')
    ylabel('Power coupled [ppm] ')
    legend('numerical','analytical')



