function[COR, Del1, Del2, Qcoef2_shift]=co_fbcfcor(Qcoef1, Qcoef2, ngrid) 
% DAMOCO Toolbox, function CO_FBCFCOR, version 17.01.11
%
% This function computes the correlation between two coupling functions
% given in terms of Fourier coefficients Qcoef1, Qcoef2. The correlation
% is defined as the maximum over all possible relative shifts of two coupling
% functions. The resolution of tested phase shifts is defined by ngrid.  
% In the computation, the coupling function of the first oscillator is taken
% as a reference, while the second one is shifted with respect to the
% first.
% 
% Form of call: 
%               COR = co_fbcfcor(Qcoef1, Qcoef2) 
%               COR = co_fbcfcor(Qcoef1, Qcoef2, ngrid) 
%               [COR, Delta1, Delta2]=co_fbcfcor(Qcoef1, Qcoef2, ngrid) 
%               [COR, Delta1, Delta2, Qcoef2_shift]=co_fbcfcor(Qcoef1, Qcoef2, ngrid) 
% Input:
%       Qcoef1, Qcoef2: The Fourier coefficients of both coupling functions
%       ngrid:          grid size, defaul is ngrid =100
% Output:
%       COR:            Maxmimal correlation of the coupling functions
%       Del1, Del2      Shifts of the own (phi2) and external (phi1) 
%                       phases which maximize the correlation between
%                       the shifted second coupling function and
%                       the first one (reference)  
%       Qcoef2_shift:   Fourier coefficients of the shifted coupling
%                       function for the second oscillator. Correlation
%                       between this function and the first one is maximal
%                       for zero shifts. 
if nargin<3; ngrid=100; end;
S = size(Qcoef1); S = S(1);
or = (S-1)/2;
Qcoef1(or+1, or+1)=0; % Autonomous frequencies are deleted if they are still present
Qcoef2(or+1, or+1)=0;
[X,Y] = meshgrid(2*pi*(0:ngrid-1)/ngrid, 2*pi*(0:ngrid-1)/ngrid);
M = zeros(size(X));
autCor1 = real(sum(sum(Qcoef1.*conj(Qcoef1)))); % Autocorrelation
autCor2 = real(sum(sum(Qcoef2.*conj(Qcoef2)))); 
for n = - or : or;
    for m = - or : or;     % Computing correlation for combinations of shifts
       M = M + (exp( 1i * (-n*X - m*Y) ) * Qcoef2(-n+or+1, -m+or+1) * Qcoef1(or+1+n, or+1+m)); 
    end;
end;
M=real(M);
Cor12 = M ./ sqrt( autCor2 .* autCor1 );
COR  = max(max(Cor12)); % computing maximal  correlation 
if nargout > 1 
    [m1, m2] = find(Cor12==COR);
    Del1 = X(m1,m2);
    Del2 = Y(m1,m2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargout >3;
    Qcoef2_shift = zeros(S,S);
    for n = - or : or;
        for m = - or : or;       % Computing coefficients of the shifted coupling function 
            Qcoef2_shift(n+or+1, m+or+1) = exp( 1i * ( n*Del1 + m*Del2 )) * Qcoef2(n+or+1, m+or+1); 
        end;
    end;
end;
end