%--------------------------------------------------------------------------
% function [ann] = FT_Winkler_coupling(lambda,w,n,n_,tilt,h0,slambda,phi)
%
% A function for Matlab which calculate the coupling into higher order
% modes from a distorted surface, using Winklers approximation (Winkler
% 1994).  The approximation from this paper has been modified to include 
% incident modes other than n=0.  The coupling is calculated from a surface
% distorted at one particular spatial frequency.
%
% lambda:   wavelength of the beam incident on the surface [m]
% w:        beam size at the surface [m]
% n:        mode number (incident mode)
% n_:       mode number (coupled mode)
% tilt:     angle of tilt of the surface [degrees]
% h0:       amplitude of surface distortion [m]
% slambda:  spatial wavelength of distortion [m]
% phi:      spatial phase [degrees]
%
% ann:      coupling coefficient
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    13.05.2012
%--------------------------------------------------------------------------
%

function [ann] = FT_Winkler_coupling(lambda,w,n,n_,tilt,h0,slambda,phi)

    baseid = 'Winkler_coupling';

    % Error messages
    % Mode/ zernike indices must be integers
    if (round(n)~=n || round(n_)~=n_)
        result='Invalid mode indices: n and n_ must be integers';
        msgid=[baseid,':checkarguments']
        error(msgid,result)
    end


    % Putting angles into radians
    phi = phi*pi/180;
    tilt = tilt*pi/180;
    
    % Calculate normaised parameters
    k = 2*pi/lambda;
    W = sqrt(2)*pi*w./slambda;
    I0 = 2*k*h0;
    alpha = sqrt(2)*k*w*tilt;

    % Define variables
    n1 = min(n,n_);
    n2 = max(n,n_);
    dn = n2-n1;
    
    % The tilt only couples between modes with dn = 1 for this
    % approximation, set kronecker delta
    if dn==1
        d_dn1 = 1;
    else
        d_dn1 = 0;
    end
    
    % Adapted ormula, taking into account different incident modes
    % General factor
    a0 = 1i * sqrt(factorial(n1)/factorial(n2)) * (1/sqrt(2^dn));
    
    % Tilt term
    A1 = alpha * factorial(n2)/factorial(n1) * d_dn1;
    
    % Spatial wavelengths term
    A2 = I0 * exp(-(W.^2)/4) .* W.^dn .* FT_LaguerrePol(n1,dn,W.^2/2) * cos(phi + dn*pi/2);
    
    % Final coupling coefficient
    ann = a0 * (A1+A2); 

end

