%--------------------------------------------------------------------------
% FTE_scattering_from_syrfaces.m
%
% Example Matlab script which illustrates how the Finesse Tools functions 
% can be used to analyse coupling into higher order modes caused by
% scattering from a distorted surface.  
%
% This example plots the coupling from a distorted surface in 3 different 
% ways, from a surface made of 3 different spatial frequencies.  The
% coupling is calculated from the total surface and from the individual
% spatial frequencies (both numerically and analytically) which are then
% added to get the total.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    13.05.2012
%--------------------------------------------------------------------------

clear all;

profile_plot = 'surface_profile';
power_plot = 'comparing_coupling_methods';
diff_plot = 'approximation_accuracy';

% Different maximum amplitudes for the overall distorted surface
A = linspace(0.01e-9,100e-9,1000);

% Beam parameters (use aligo cavity parameters)
w0 = 6.2e-2;         % Beam size at the ETM
z = 0;
lambda = 1064e-9;
k = 2*pi/lambda;
nr=1;
gp = FT_init_gauss_param(1064e-9,nr,'w0',w0,'z',z);

Nin=0;

% Set range 
R = 17e-2;          % Radius of aligo cavity mirrors
r = linspace(-2*R,2*R,1000);
rstep = abs(r(1)-r(2));
scaling = 1e-9;
r0 = (length(r)+1)/2;

% Calculate incident field
y = 0;
x = r;
field = FT_HG_field(gp,Nin,0,x,y,[0,0,0]);

% Setup 3 individual surface distortions with different spatial
% frequencies, relative amplitudes and phases
l1 = 2*R;
l2 = 3*R/2;
l3 = R;

a1 = 1;
a2 = 1;
a3 = 1;

phi1 = 360*rand();
phi2 = 360*rand();
phi3 = 360*rand();

% Calculate individual distortions
z1 = a1*FT_harmonic_surface_distortion(r,l1,phi1);
z2 = a2*FT_harmonic_surface_distortion(r,l2,phi2);
z3 = a3*FT_harmonic_surface_distortion(r,l3,phi3);

% Calculate total surface distortion
ztotal = z1 + z2 + z3;

% Normalising
rms = sqrt(sum(ztotal.^2)/length(ztotal));
z1 = z1/rms;
z2 = z2/rms;
z3 = z3/rms;
ztotal = ztotal/rms;

% Produce and plot surface profile
profile = FT_new_surface_profile(ztotal,r0,rstep,scaling);
figure
    FT_plot_surface_profile(profile,0);

% Calculate coupling for range of amplitudes
% Numeric -> calculate coupling numerically from seperate spatial 
% frequencies (C_123)
% Numeric -> calculate coupling from combined surface (C_total)
for i=1:length(A)
    
    % Scale distortions
    z1_ = A(i)*z1;
    z2_ = A(i)*z2;
    z3_ = A(i)*z3;
    ztotal_ = A(i)*ztotal;
    
    % Calculate reflected field
    field1 = field.*exp(2*1i*k*z1_);
    field2 = field.*exp(2*1i*k*z2_);
    field3 = field.*exp(2*1i*k*z3_);
    field_total = field.*exp(2*1i*k*ztotal_);
    
    % Calculate the coupling into modes with different n
    for n=0:5
        
        field_ref = FT_HG_field(gp,n,0,x,y,[0,0,0]);
        
        c1 = FT_conv_fields(field_ref,field1);
        c2 = FT_conv_fields(field_ref,field2);
        c3 = FT_conv_fields(field_ref,field3);
        
        C_123(n+1,i) = c1+c2+c3;
        
        C_total(n+1,i) = FT_conv_fields(field_ref,field_total);
        
    end
    
end


% Calculate anlytic coupling
% Analytic -> seperate spatial frequencies
for j=1:length(A)    
    for n=0:5
    
        % Scale amplitudes
        A1 = A(j)*a1/rms;
        A2 = A(j)*a2/rms;
        A3 = A(j)*a3/rms;
        
        % Calculate analytic coupling using Winkler approximation
        ca1 = FT_Winkler_coupling(lambda,w0,Nin,n,0,A1,l1,phi1);
        ca2 = FT_Winkler_coupling(lambda,w0,Nin,n,0,A2,l2,phi2);
        ca3 = FT_Winkler_coupling(lambda,w0,Nin,n,0,A3,l3,phi3);
        
        C_winkler(n+1,j) = ca1+ca2+ca3;
        
    end
    
end

% Calculate coupled power in ppm
P_123 = 1e6*abs(C_123).^2;
P_winkler = 1e6*abs(C_winkler).^2;
P_total = 1e6*abs(C_total).^2;

% Put amplitude in nm
Anm = A*1e9;

% Different colours for different plots
colour = ['r' 'b' 'g' 'c'];

% Plot results for first 3 order couplings to compare the different methods
figure
    hold on;
    for n=0:3
        i=n+1;
        if n>=Nin
            i = n+2;
        end
        plot(Anm,P_winkler(i,:),colour(n+1),'LineWidth',2)
        st(n+1,:) = sprintf('n = %g',i-1);
    end

    for n=0:3
        i=n+1;
        if n>=Nin
            i = n+2;
        end
        plot(Anm,P_123(i,:),'k-.','LineWidth',2)
        plot(Anm,P_total(i,:),'k--','LineWidth',1)
    end

    legend(st(1,:),st(2,:),st(3,:),st(4,:),'numerical sum','numerical total','Location','NorthWest')

    set(gca,'XScale','Log')
    set(gca,'YScale','Log')

    xlabel('rms [nm] ')
    ylabel('Power coupled [ppm] ')
    grid on;
    hold off;

% Calculate relative difference in power scattered
dP = abs((P_total - P_winkler)./(P_total));

% Plot difference
colour2 = ['r' 'b' 'g' 'c'];

figure
    hold on;
    for n=0:3
        i=n+1;
        if n>=Nin
            i = n+2;
        end
        plot(Anm,dP(i,:),colour2(n+1),'LineWidth',2)
        st(n+1,:) = sprintf('n = %g',i-1);
    end
    legend(st(1,:),st(2,:),st(3,:),st(4,:))

    set(gca,'XScale','Log')
    set(gca,'YScale','Log')

    xlabel('rms [nm] ')
    ylabel('Fractional power difference')
    grid on;
