function [sigfc1, sigfc2, sigma1, sigma2, arg] = co_fbsolv(Fcoef1, Fcoef2, ngrid)
% DAMOCO Toolbox, function CO_FBSOLV, version 17.01.11
% Given the coefficients Fcoef1, Fcoef2 of the Fourier expansion of the
% protophase dynamics, computed by co_fexp2, this function computes 
% Fourier coefficients sigfc1, sigfc2 of the phase transformation. 
% Optionally it returns the transformation functions, computed on the grid
% of size ngrid.
% 
% Form of Call:
%           [sigfc1,sigfc2] =  co_fbsolv(Fcoef1, Fcoef2)
%           [sigfc1,sigfc2,sigma1,sigma2,arg] = co_fbsolv(Fcoef1,Fcoef2)
%           [sigfc1,sigfc2,sigma1,sigma2,arg] = co_fbsolv(Fcoef1,Fcoef2,ngrid)
%
% INPUT
%  Fcoef1:   Coefficients of Fourier expansion of the protophase dynamics 
%            of the 1st oscillator, produced by co_fexp2. 
%  Fcoef2:   Coefficients of protophase dynamics of the 2nd oscillator,  
%            produced by co_fexp2.
%  ngrid:    grid size is required only if sigma1, sigma 2 are computed
%            Default value is ngrid = 100.
%   
% OUTPUT
%  sigfc1:   Fourier coefficients of the 1st phase transformation function
%  sigfc2:   Fourier coefficients of the 2st phase transformation function
%  sigma1:   Transformation function for the 1st oscillator
%  sigma2:   Transformation function for the 2nd oscillator
%            These function can be plotted as  plot(arg, sigma2)
%  arg:      argument for plotting sigma1 / sigma2
%
or = size(Fcoef1); or = ((or(1)-1)/2);                   % Compute order, or, from matrix Fcoef1;
[S] = fsolve(@nls, zeros(4*or+2,1), [], Fcoef1, Fcoef2); % Solving non-linear-system for coefficients in vector S
sigfc1 = S(1:2*or);                                      % Reorganzing the coefficients of the solution S
sigfc1 = [sigfc1(1:or); 1; sigfc1(or+1:end)];            
sigfc2 = S(2*or+1:4*or);                                 
sigfc2 = [sigfc2(1:or); 1; sigfc2(or+1:end)];            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%     Computing transformation functions sigma1, sigma2
if nargout == 5;
    if nargin ==2; ngrid=100; end;
    or=(length(sigfc1)-1)/2;
    arg=2*pi*(0:ngrid-1) / (ngrid-1);
    sigma1=zeros(1,ngrid); sigma2=sigma1;
    for n = -or:-1;
        sigma1 = sigma1 + 2*real(sigfc1(n+or+1)*exp(1i*n*arg)); 
        sigma2 = sigma2 + 2*real(sigfc2(n+or+1)*exp(1i*n*arg)); 
    end;
    sigma1 = sigma1+1;
    sigma2 = sigma2+1;
end;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x] = nls(X, Fcoef1, Fcoef2)           % Defining non-linear system from coefficients Fcoef1, Fcoef2
or = size(Fcoef1); or = ((or(1)-1)/2);

S = X(1:2*or);
S = [S(1:or);1;S(or+1:end)];

R = X(2*or+1:4*or);
R = [R(1:or);1;R(or+1:end)];
omeg1 = X(end-1);
omeg2 = X(end);

M = S*R.';

Fcoef1 = conj(Fcoef1); Fcoef2 = Fcoef2';
x(1) = sum(sum(M.*Fcoef1))-omeg1;
x(2) = sum(sum(M.*Fcoef2))-omeg2;

J = Fcoef1; K = Fcoef2;
for n = 1:or;
    J(1,:) = []; K(:,1) = [];
    J = [J;zeros(1,2*or+1)]; K = [K zeros(2*or+1,1)]; 
    x = [x sum(sum(M.*J))];
    x = [x sum(sum(M.*K))];
end;

J = Fcoef1; K = Fcoef2;
for n = 1:or;
    J(end,:) = []; K(:,end) = [];
    J = [zeros(1,2*or+1); J]; K = [zeros(2*or+1,1) K]; 
    x = [x sum(sum(M.*J))];
    x = [x sum(sum(M.*K))];
end;
end