function [phi1,phi2,omeg1,omeg2,q1,q2,Nrmq1,Nrmq2,dirin] =...
    co_ittransf2(theta1,theta2,or,fsamp,ngrid)
% DAMOCO Toolbox, function CO_ITTRANSF2, version 17.01.11
% This high-level function performs all steps of the phase transformation at once,
% using the iteration method.
% Given two protophases theta1 and theta2, the order of the Fourier expansion, or, 
% the sampling frequency of the data, fsamp,  and the grid size ngrid, this
% function computes: 
%      the true phases phi1 and phi2, 
%      the true coupling functions q1 and q2 on a grid of size ngrid, 
%      the frequencies (constant terms of coupling functions) omeg1 and omeg2,
%      the norms of the coupling functions, norm1 and norm2, 
%      and the directionality index dirin.
% Form of call: 
%      [phi1,phi2,omeg1,omeg2,q1,q2,Nrmq1,Nrmq2,dirin] = co_ittransf2(theta1,theta2,or,fsamp,ngrid)
%      [phi1,phi2,omeg1,omeg2,q1,q2,Nrmq1,Nrmq2,dirin] = co_ittransf2(theta1,theta2,or,fsamp)
% Input:
%       theta1:     protophase of the 1st oscillator
%       theta2:     protopahse of the 2nd oscillator
%       or:         order of the Fourier expansion: or = 10 is a good first choice.
%       fsamp:      sampling frequency
%       ngrid:      size of the grid to compute q1, q2; default value is ngrid=100
% Output: 
%       phi1:       true phase of the 1st oscillator
%       phi2:       true phase of the 2nd oscillator
%       omeg1:      autonomous frequency (constant term of the coupling function) of the 1st oscillator
%       omeg2:      autonomous frequency (constant term of the coupling function) of the 2nd oscillator
%       q1:         coupling function of the 1st oscillator, computed on an a grid 
%                   Note: It does not contain the constant term, i.e. the
%                   phase dynamics is given by dphi1/dt = omeg1 + q1(phi1,phi2).
%       q2:         coupling function of the 2nd oscillator, computed on a
%                   grid, without constant term.
%       Nrmq1:      the norm of q1.
%       Nrmq2:      the norm of q2.
%       dirin:      the index of directionality of interaction between both oscillators: 
%                   dirin = 1: unidirectional coupling from 1 to 2. 
%                   dirin = -1: unidirectional coupling from 2 to 1.
%                   dirin = 0: equally strong bidirectional coupling  
%
if nargin < 5       % default value for the grid size
    ngrid=100; 
end         
co_testproto(theta1, theta2);        % Testing protophases
pr_theta1 = co_fbtransf1(theta1);    % Preprocessing with univariate Fourier based transform
pr_theta2 = co_fbtransf1(theta2);    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%       WARNING,  if protophases are too strongly correlated
[dummy, maxind]=co_maxsync(pr_theta1, pr_theta2, 2*or); clear dummy;
if maxind > 0.6;
    disp(['Synchronisation index is ' num2str(maxind)]);
    disp('Protophases are strongly correlated, the results may not be reliable!');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%           Iteration Method
[Fcoef1,Fcoef2,f1,f2]=co_fexp2(pr_theta1,pr_theta2,or,fsamp,ngrid); 
niter=4;       % number of iterations; 4 is usually sufficient
showiter=1;    % intermediate results are shown if showiter >0
[q1,q2,omeg1,omeg2,Nrmq1,Nrmq2,sigma1,sigma2]=co_itersolv(f1,f2,niter,showiter);
q1=q1-omeg1; q2=q2-omeg2;            % coupling functions do not contain constant terms
phi1=co_gth2phi(pr_theta1,sigma1);   % final transform of 1D-transformed protophases
phi2=co_gth2phi(pr_theta2,sigma2);
dirin = co_dirin(Nrmq1, Nrmq2, omeg1,omeg2);        % Direction of coupling
end