%--------------------------------------------------------------------------
% function [K,nm_in,nm_out] = FT_map_knm( map, gp_in, gp_out, maxtem,err)
%
% Calculates the coupling coefficient for a beam interacting with a
% surface defined by a typical FINESSE map file
%
% map:    Map object
% gp_in:  Incoming beam parameters
% gp_out: Outgoing beam parameters
% maxtem: Maximum mode number to compute upto
% err:    Absolute error for integration routine
% 
% Output
% K:            Matrix of coupling coefficents
% nm_in,nm_out: Incoming and outgoing modes that match upto matrix K
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
%
% Daniel Brown    25.11.2012
%--------------------------------------------------------------------------
%
function [K,nm,nm_] = FT_map_knm( map, gp_in, gp_out, maxtem,err)
    MAX =  base2dec(sprintf('%i%i',maxtem,maxtem),maxtem+1);
    K = zeros(MAX,MAX);
    curr = 0;
    
    X = ((0:(size(map.data,1) -1 ))- map.x0).*map.xstep;
    Y = ((0:(size(map.data,2) -1 ))- map.y0).*map.ystep;
    
    A = @(x,y) exp(1i.*2.*pi.*interp2(X,Y,map.data'.*map.scaling,x,y)./gp_in.lambda);
    
    currnm = 0;
    
    nm = cell((MAX+1),1);
    nm_ = cell((MAX+1),1);
    
    for n_in = 0:maxtem
        for m_in = 0:maxtem
        
            currnm = currnm + 1;
            nm(currnm) = cellstr(sprintf('%i%i',n_in,m_in));
            currnm_ = 0;
            
            for n_out = 0:maxtem
                
                for m_out = 0:maxtem
                    currnm_ = currnm_ + 1;
                    
                    if currnm == 1
                        nm_(currnm_) = cellstr(sprintf('%i%i',n_out,m_out));                    
                    end
                    
                    fprintf('Calculating %s->%s = ',char(nm(currnm)),char(nm_(currnm_)));
                    
                    i = base2dec(char(nm(currnm)),maxtem+1);
                    j = base2dec(char(nm_(currnm_)),maxtem+1);
                    
                    K(i+1,j+1) = FT_knm(gp_in,gp_out,A,n_in,n_out,m_in,m_out,err);
                    fprintf('%.6e %.6fi\n',real(K(i+1,j+1)),imag(K(i+1,j+1)));
                    curr = curr+1;
                end
            end
        end
    end
end

