function[Qcoef1, Qcoef2, Qcoef3]=co_fcpltri(phi1, phi2, phi3, Dphi1, Dphi2, Dphi3, N)
% DAMOCO Toolbox, function CO_FCPLTRI, version 06.03.14
% Given phases of 3 oscillators and their derivatives, the function yields the coupling functions
% q1(phi1, phi2, phi3),
% q2(phi2, phi3, phi1), and
% q3(phi3, phi1, phi2)
% via fitting a Fourier series: exp( i(n*phi1 + m*phi2 +k*phi3) ).
% The results for each oscillator are given by 3-dimensional matrix of coefficients
% Qcoef1(n,m,k)
% Qcoef2(m,k,n)
% Qcoef3(k,n,m)
% where the indices are associated with the phases as indicated above:
% n to ph1,
% m to phi2
% k to phi3
% Hence, the first dimension of each matrix always corresponds to the phase, which dynamics
% is represented by the coefficients.
% Form of call:
% [Qcoef1, Qcoef2, Qcoef3] = co_fcpltri(phi1, phi2, phi3, or, safr)
%
% Input: phi1: phase of the 1st oscillator,
% phi2: phase of the 2nd oscillator,
% phi3: phase of the 3rd oscillator,
% Dphi1: derivative of phi1
% Dphi2: derivative of phi2
% Dphi3: derivative of phi3
% N: order of the Fourier series
% Output: Qcoef1, Qcoef2, Qcoef3 are Fourier coefficients of the coupling functions
phi1 = unwrap(phi1); % unwrapped phases are needed here
phi2 = unwrap(phi2);
phi3 = unwrap(phi3);
A=zeros(4*N+1,4*N+1,4*N+1); % This matrix contains the coefficients A(n+p),(m+q),(k+r)
% for the linear system of equations to
% obtain the coefficients Qcoef(n,m,k)
ncf=2*N+1; ncf2=ncf*ncf; % number of coefficients in each dimension
for n = -2*N : 2*N;
for m = - 2*N : 2*N;
for k = -2*N : m;
A( n+ncf, m+ncf, k+ncf ) = mean( exp( 1i*( n*phi1 + m*phi2 + k*phi3 ) ) );
A( -n+ncf, -m+ncf, -k+ncf ) = conj (A( n+ncf, m+ncf, k+ncf ));
end;
end;
end;
B1 = zeros(ncf^3, 1); % This vector contains the coefficients B1(n,m,k) for the linear equation system to
% obtain the coefficients Qcoef1(n,m,k) of the phase dynamics of phi1
B2 = B1; % This vector contains the coefficients B2(n,m,k) for the linear equation system to
% obtain the coefficients Qcoef2(m,k,n) of
% the phase dynamics of phi2
B3 = B1; % This vector contains the coefficients B3(n,m,k) for the linear equation system to
% obtain the coefficients Qcoef3(k,n,m) of the phase dynamics of phi2
C = zeros(ncf^3,ncf^3); % The elements of the matrix A are reorganized in C to match the requirements of the
% MATLAB function to solve systems of linear equations
ind = 1;
for r = -N : N;
for s = -N : N;
for q = -N : N;
EXP=exp( -1i*( r*phi1 + s*phi2 + q*phi3) );
B1(ind) = mean( Dphi1 .* EXP );
B2(ind) = mean( Dphi2 .* EXP );
B3(ind) = mean( Dphi3 .* EXP );
clear EXP;
ind = ind+1;
for n = -N : N;
for m = -N : N;
for k = -N : N;
C((r+N)*ncf2 + (s+N)*ncf + (q+N+1) , (n+N)*ncf2 + (m+N)*ncf + (k+N+1) ) = A( (n-r)+ncf, (m-s)+ncf, (k-q)+ncf );
% Elements of the matrix A are reorganized in C
end;
end;
end;
end;
end;
end;
clear A phi* Dphi*
coeff1 = C\B1; % solving linear system to obtain coefficients
coeff2 = C\B2; % solving linear system to obtain coefficients
coeff3 = C\B3; % solving linear system to obtain coefficients
clear C B*
%%% rewriting the coefficients in NxMxK matrix-form
Qcoef1 = zeros(ncf,ncf,ncf); Qcoef2 = Qcoef1; Qcoef3 = Qcoef1;
for n = 1 : ncf;
for m = 1 : ncf;
for k = 1 : ncf;
ind=(n-1)*ncf2 + (m-1)*ncf + k;
Qcoef1(n, m, k) = coeff1(ind);
Qcoef2(m, k, n) = coeff2(ind);
Qcoef3(k, n, m) = coeff3(ind);
end;
end;
end;
end