%--------------------------------------------------------------------------
% function [K] = FT_knm(gp_in, gp_out, A, n_in,n_out,m_in, m_out, err)
%
% Calculates the coupling coefficient for a beam interacting with a
% surface described by the function A(x,y). This function can be complex.
% Calclates coupling coefficient for mode [n_in,m_in]->[n_out,m_out].
%
% gp_in:  Incoming beam parameters
% gp_out: Outgoing beam parameters
% A:      A function of x and y that describes the surface
% n_in:   Incoming x mode
% n_out:  outgoing x mode
% m_in:   incoming y mode
% m_out:  outgoing y mode
% err:    Absolute error for integration routine
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Daniel Brown    25.11.2012
%--------------------------------------------------------------------------
%
function [K] = FT_knm(gp_in, gp_out, A, n_in,n_out,m_in, m_out, err)

%set bounds to incluce all but 1ppm of beam
xrange=5*max(gp_in.w(1),gp_out.w(1))*sqrt(max(n_in,n_out)+1);
yrange=5*max(gp_in.w(2),gp_out.w(2))*sqrt(max(m_in,m_out)+1);

in =  @(x,y)(FT_unm(gp_in.lambda./gp_in.nr, gp_in.w0(1),  gp_in.w0(2),  x, y, gp_in.z(1), gp_in.z(2), n_in, m_in));
out = @(x,y)(FT_unm(gp_out.lambda./gp_out.nr,gp_out.w0(1), gp_out.w0(2), x, y, gp_out.z(1), gp_out.z(2), n_out, m_out));

Func = @(x,y) in(x,y).* A(x,y).*conj(out(x,y));

K = integral2(Func, -xrange,xrange,-yrange,yrange,'AbsTol',err);

end
