%      SyncKuramotoDemo.M
%   
%  Figures to demonstrate the synchronisation of
%  Kuramoto oscillators
%     partly from KuramotoIMP.m
%  
% J. Lunze
% 
% 17.2.2018
% Version v. 27.12.2018
% for 2nd edition: 29.2.2021
%
echo off
clear
close all


%%   Oscillator behaviour
%
figure(1)
omega=1;
theta0=2;
Tend=15;
Time=0:0.1:Tend;
Theta=theta0+omega*Time;
subplot(6,1,[1 2])
plot(Time, sin(Theta),'b')
axis([0, Tend, -1.3, 1.3]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
yticks([-1  1]);
latextitle('$${\tt SyncKuramotoDemo:}$$ Oscillator behaviour');
hold off;
%
subplot(6,1,[3 4])
plot(Time, Theta,'b');
latexylabel('$$\theta_i$$ in rad');
rotateY;
latexxlabel('$$t$$');
axis([0, Tend, 0, 17]);
yticks([2  10 ]);
xticks([0 5 10 15]);
hold off;

%%  Consensus of uniform Kuramoto oscillators
%
figure(2)
omega=1;
N1=10;
Tend=15;
ks=0.03;
%  complete couplings
L=-ones(N1, N1);
for k1=1:N1
    L(k1,k1)=N1-1;
end
%  Network model
Abar=-ks*L;
Bbar=zeros(N1,1);
Cbar=eye(N1,N1);
Network=ss(Abar, Bbar, Cbar, 0);
%x0=4*rand(N1,1)-2;    %  initial states in [-2, 2]
theta0=[   
    0.0577
    1.5371
    0.3521
   -1.9810
   -1.2005
   -0.3722
    0.9948
    1.9023
    1.1599
   -0.7259];
xbar=sum(theta0)/N1;      %  consensus state
[Y,Time]=initial(Network,theta0,Tend);
subplot(6,1,[1 2]);
plot(Time, sin(Y+omega*Time*ones(1,N1)), 'b');
axis([0 Tend -1.3 1.3]);
latexylabel('$$y_i\;$$');
rotateY;
yticks([-1 1]);
xleer;
latextitle('$${\tt SyncKuramotoDemo:}$$ Uniform oscillator network');
hold off
%
%   output of the Kuramoto oscillators
subplot(6,1,[3 4])
plot(Time,Y,'b');
hold on
latextext(15.3, 0, '$$\bar{\theta}$$');
%latexxlabel('$$t$$');
latexylabel('$$\tilde{\theta}_i$$ in rad');
rotateY;
xticks([0 5 10 15]);
%yticks([-2 0 2]);
axis([0 Tend -2.2 2.2]);
latexxlabel('$$t$$');
xticks([0 5 10 15]);
hold off

%%  Arrow diagram
%
figure(18)
subplot(4,6,[1 2 7 8])
%   unit circle:
k1=0;
for alpha=0:0.1:2*pi+0.09
    k1=k1+1;
    x(k1)=cos(alpha);
    y(k1)=sin(alpha);
end
%
%    results of figure 2, initial state
plot(cos(theta0), sin(theta0),'bo');
hold on
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
latexylabel('Im');
rotateY;
yticks([-1 0 1]);
xticks([-1 0 1]);
latextext(0.3, 1.3, '$$t=0$$');
hold off
%
%    results of figure 2  for t=6
subplot(4,6,[3 4 9 10])
time=6;
ktime=ceil(length(Time)/Tend*time);
plot(cos(Y(ktime,:)), sin(Y(ktime,:)),'bo');
hold on
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
%latexylabel('Im');
%rotateY;
yleer;
yticks([]);
xticks([-1 0 1]);
latextext(0.3, 1.3, ['$$t=\,$$' num2str(time)]);
hold off
%
%    results of figure 2 (Tend)
subplot(4,6,[5 6 11 12])
plot(cos(Y(end,:)), sin(Y(end,:)),'bo');
hold on
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
%latexylabel('Im');
%rotateY;
yleer;
yticks([]);
xticks([-1 0 1]);
latextext(0.2, 1.3, ['$$t=\,$$' num2str(Tend)]);
hold off


%%  Consensus of non-uniform Kuramoto oscillators
%                    KuramotoIMP.m, figure 5
%                  = KuramotoIJC.m, figure 1
figure(3)
N1=10;
omega0=0.25:0.17:1.9;  % only the first N1 elements will be used
Tend=15;
ks=0.1;
%  complete couplings
L=-ones(N1, N1);
for k1=1:N1
    L(k1,k1)=N1-1;
end
%  Network model
%  Oscillators:
A01=[0 1; 0 0];
B01=[1; 0];
C01=[1 0];
n=2;
%  Controller selection to ensure the stability of 
%  the interconnected followers
Abar=-kron(ks*L, B01*C01)+kron(diag(ones(N1,1)), A01);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C01);
Network=ss(Abar, Bbar, Cbar, 0);
%  initial states
x0=zeros(N1*n,1);
for k1=1:N1
    x0((k1-1)*n+1)=theta0(k1);
    x0((k1-1)*n+2)=omega0(k1);
end
[Y,Time]=initial(Network,x0,Tend);
subplot(6,1,[1 2]);
plot(Time, sin(Y), 'b');
axis([0 Tend -1.3 1.3]);
latexylabel('$$y_i\;$$');
rotateY;
yticks([-1 1]);
xleer;
latextitle('$${\tt SyncKuramotoDemo:}$$ Non-uniform oscillator network');
hold off
%
%   output of the Kuramoto oscillators
subplot(6,1,[3 4])
plot(Time,Y,'b');
hold on
%latexxlabel('$$t$$');
latexylabel('$$\theta_i$$ in rad');
rotateY;
xticks([0 5 10 15]);
yticks([0  10]);
axis([0 Tend -3 15]);
latexxlabel('$$t$$');
xticks([0 5 10 15]);
hold off
%
%  Determination of the remaining phase shifts
%  \bar{\theta}(t) - \theta_i(t)  ->  -1/k/N (\bar{omega} - \omega_i)
omegabar=sum(omega0(1:N1))/N1;
for k1=1:N1
    deltatheta(k1)=abs(1/ks/N1*(omegabar-omega0(k1))); 
end
epsilon=max(deltatheta);
%   compare with the oscillator behaviour
Theta=Y(end,:);
thetabar=sum(Theta)/N1;
max(thetabar-Theta);   % to compare with epsilon above


%%  Arrow diagram
%
figure(19)
subplot(4,6,[1 2 7 8])
%   unit circle:
k1=0;
for alpha=0:0.1:2*pi+0.09
    k1=k1+1;
    x(k1)=cos(alpha);
    y(k1)=sin(alpha);
end
%
%    results of figure 3, time=13
time=13;
ktime=ceil(length(Time)/Tend*time);
plot(cos(Y(ktime,:)), sin(Y(ktime,:)),'bo');
hold on
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
latexylabel('Im');
rotateY;
yticks([-1 0 1]);
xticks([-1 0 1]);
latextext(0.2, 1.3, ['$$t=\,$$' num2str(time)]);
hold off


%%  Synchronisation region of extended Kuramoto oscillators
%
figure(4)
ke=1;
kr=1;
kbar=kr/ke;
alphamin=-5;
alphamax=5;
betamin=-5;
betamax=5;
deltaalpha=0.01;
deltabeta=0.05;
k1=0;
Stab=[];
for alpha=0:deltaalpha:kbar-deltaalpha
    k1=k1+1;
    Stab(k1,:)=[alpha, sqrt(alpha*alpha*alpha*ke/(kr*kr-alpha*ke)),... 
                       -sqrt(alpha*alpha*alpha*ke/(kr*kr-alpha*ke))];
end
fill([Stab(:,1); kbar; 5; 5; 0], [Stab(:,2); 10; 10; 0; 0], [0.9 0.9 0.9])
hold on
fill([Stab(:,1); kbar; 5; 5; 0], [Stab(:,3); -10; -10; 0; 0], [0.9 0.9 0.9])
latexxlabel('$$\hspace{2cm}$$ Re');
xticks([0])
%
%latextext(kbar, 5, '$$+$$');
latextext(0.5, -4.6, '$$-k^2_{\rm r}/k_{\rm e}$$'); 
yticks([0]);
latexylabel('Im');
rotateY;
axis([-3 5 -4 4])
axis('square');
plot([0 0],[-5 5],'k:', [-5 5],[0 0],'k:');
hold off


%%   Complete synchronisation of Kuramoto oscillators
%           KuramotoIJC.m   figure 2
%
figure(5)
N1=10;
Tend=15;

%  Controller selection to ensure the stability of 
%  the interconnected followers
ke= 10/4;
kr= 10/4;
A01=[0 1; 0 0];
B01=[ke; kr];
C01=[1 0];
D01=0;
n=2;
%
%  complete couplings
L=-ones(N1, N1);
for k1=1:N1
    L(k1,k1)=N1-1;
end
Lhat=L/(N1-1);
ks=0.4;                % coupling strength
Abar=-kron(ks*Lhat, B01*C01)+kron(diag(ones(N1,1)), A01);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C01);
Dbar=zeros(N1, 1);
OscillatorDyn=ss(Abar, Bbar, Cbar, Dbar);
%  initial states as for static coupling network with different omega_i
[Y, Time]=initial(OscillatorDyn, x0, Tend);
Y5=Y(end,:);
subplot(2,1,1)
plot(Time, sin(Y), 'b');
hold on
plot([0 Tend],[0 0],'k:');
axis([0 Tend -1.1 1.1]);
latexxlabel('$$t$$');
latexylabel('$$y_i\;$$');
rotateY;
yticks([-1 0 1]);
xticks([0 5 10 15]);
latextitle('$${\tt SyncKuramotoDemo:}$$ Non-uniform extended oscillator network');
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Kuramoto with nonlinear couplings
%    Uniform oscillators
%     KuramotoDiff.m   figure 10
%
figure(6)
N1=10;
N=10;       % parameter of the simulink model
omega=3*ones(N1,1);
k=0.3;                    % 0.7 coupling strength
Tend=20;
dt=0.01;
T=0:dt:Tend;
x0=[-1.4234   
    -0.1207    
    3.3975   
    -2.3865   
    -2.3054   
    -1.8221    
    2.8593    
    3.1490
    2.8536    
    1.3605];    % initial state for convergence within 20 time units
Theta=[];
r=[];
sim('KuramotoDiffModel');
%   N oscillator outputs
%   in vector Theta
%
subplot(6,2,[1 3])
plot(T, sin(Theta), 'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
%latextitle('$${\tt SyncKuramotoDemo:}$$ Uniform Kuramoto oscillators');
latextext(2, 1.4, '$$k=0.3$$');
yticks([-1 0 1]);
hold off;
%  alternative representation
subplot(6,2,[5 7])
plot(T, Theta, 'b')
xleer;
latexylabel('$$\theta_i$$');
rotateY;
axis([0 Tend 0 60]);
yticks([0 50 100]);
%  Average model
subplot(6,2,[9 11])
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1, 1], 'k:');
yticks([0 1]);
xticks([0 10 20 30 40]);
hold off
%
%    another coupling strength
k=0.7;                    % 0.7 coupling strength
x0=2*[   
    0.0577
    1.5371
    0.3521
   -1.9810
   -1.2005
   -0.3722
    0.9948
    1.9023
    1.1599
   -0.7259];    % uniformly distributed initial phases
Theta=[];
r=[];
sim('KuramotoDiffModel');
%   N oscillator outputs
%   in vector Theta
%
subplot(6,2,[2 4])
plot(T, sin(Theta), 'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i$$');
rotateY;
latextext(2, 1.4, '$$k=0.7$$');
yticks([-1 0 1]);
hold off;
%  alternative representation
subplot(6,2,[6 8])
plot(T, Theta, 'b')
xleer;
latexylabel('$$\theta_i$$');
rotateY;
axis([0 Tend 0 60]);
yticks([0 50 100]);
hold off;
%  Average model
subplot(6,2,[10 12])
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],':', [0 Tend], [1, 1], ':');
yticks([0 1]);
xticks([0 10 20 30 40]);
hold off

%%    Non-uniform Kuramoto oscillators with nonlinear couplings
%      KuramotoDiff.m    figure 1
%
figure(7)
N=10;
% x0=2*pi*rand(N,1);
x0=[1.4752
    2.2190
    5.1597
    0.0968
    0.2703
    1.0618
    4.0785
    4.5975
    4.0699
    2.8332];
x0old=x0;
%omega=2+2*rand(N,1);
omega=2.1:0.2:4;
omegaold=omega;
dt=0.1;
Tend=15;
T=0:dt:Tend;
%
%   Small coupling coefficient: oscillators follow their own frequency
%
k=0.1;                    % coupling strength
klow=k;
Theta=[];
r=[];
sim('KuramotoDiffModel');
Y7=Theta(end,:);
%   N oscillator outputs
%   in vector Theta
%
subplot(6,2,[1 3])
plot(T, sin(Theta), 'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
%latextitle('$${\tt SyncKuramotoDemo:}$$ Uniform Kuramoto oscillators');
latextext(2, 1.4, '$$k=0.1$$');
yticks([-1 0 1]);
hold off;
%  alternative representation
subplot(6,2,[5 7])
plot(T, Theta, 'b')
xleer;
latexylabel('$$\theta_i$$');
rotateY;
axis([0 Tend 0 60]);
yticks([0 50 100]);
%  Average model
subplot(6,2,[9 11])
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1, 1], 'k:');
yticks([0 1]);
xticks([0 10 20 30 40]);
hold off
%    another coupling strength
k=1.4;                    % 0.7 coupling strength
khigh=k;
Tend=5;
T=0:dt:Tend;
Theta=[];
r=[];
sim('KuramotoDiffModel');
Y8=Theta(end,:);
%   N oscillator outputs
%   in vector Theta
%
subplot(6,2,[2 4])
plot(T, sin(Theta), 'b')
axis([0, 5, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i$$');
rotateY;
latextext(.5, 1.4, '$$k=1.4$$');
yticks([-1 0 1]);
hold off;
%  alternative representation
subplot(6,2,[6 8])
plot(T, Theta, 'b')
xleer;
latexylabel('$$\theta_i$$');
rotateY;
axis([0 5 0 20]);
yticks([0 20 40]);
hold off;
%  Average model
subplot(6,2,[10 12])
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, 5, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1, 1], 'k:');
yticks([0 1]);
xticks([0 1 2 3 4 5]);
hold off



%%   No synchronisation occurs for specific initial states
%    for uniform oscillators
%       from KuramotoNonSync.m   figure 1
%      
figure(8)
N=3;
omega=0.5*ones(N,1);
dt=0.1;
Tend=30;
T=0:dt:Tend;
x0=[0;
    2/3*pi;
    4/3*pi];
k=0.3;                    % coupling strength: k/3
klow=k;
sim('KuramotoDiffModel3TS');
%   3 oscillator outputs
%   in vector Theta
%
subplot(6,1,[1 2])
plot(T, sin(Theta), 'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
latextitle('$${\tt SyncKuramotoDemo:}$$ Non-synchronisable phases');
yticks([-1 0 1]);
hold off;
%
%    Slightly changed initial state
x0=[0.1;
    2/3*pi;
    4/3*pi];
sim('KuramotoDiffModel3TS');
%   3 oscillator outputs
%   in vector Theta
%
subplot(6,1,[3 4])
plot(T, sin(Theta), 'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
latexylabel('$$y_i\;$$');
rotateY;
xleer;
latexxlabel('$$t$$');
yticks([-1 0 1]);
hold off;
%
% Anti-phase synchronisation
x0=[0; pi; pi];
sim('KuramotoDiffModel3TS');
%   3 oscillator outputs
%   in vector Theta
%
subplot(6,1,[5 6])
plot(T, sin(Theta(:,1)), 'r')
hold on
plot(T, sin(Theta(:,2)), 'b--')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
latexylabel('$$y_i\;$$');
rotateY;
latexxlabel('$$t$$');
yticks([-1 0 1]);
xticks([0 10 20 30]);
hold off;


%%   Illustration of the order parameter in the arrow diagram
%     for the oscillator behaviour of figures 5 and 7
%    Y5 from figure 5
%    Y7 from figure 7 (left)
%    Y8 from figure 7 (right)
figure(9)
subplot(4,6,[1 2 7 8])
%   unit circle:
k1=0;
for alpha=0:0.1:2*pi+0.09
    k1=k1+1;
    x(k1)=cos(alpha);
    y(k1)=sin(alpha);
end
%
%    results of figure 5
N1=10;
plot(cos(Y5), sin(Y5),'bo');
hold on
plot([0 real(sum(exp(i*Y5))/N1)],...
     [0 imag(sum(exp(i*Y5))/N1)],'b');
plot(real(sum(exp(i*Y5))/N1), imag(sum(exp(i*Y5))/N1), 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
latexylabel('Im');
rotateY;
yticks([-1 0 1]);
xticks([-1 0 1]);
hold off
%
%    results of figure 7 (left)
subplot(4,6,[3 4 9 10])
N1=10;
plot(cos(Y7), sin(Y7),'bo');
hold on
plot([0 real(sum(exp(i*Y7))/N1)],...
     [0 imag(sum(exp(i*Y7))/N1)],'b');
plot(real(sum(exp(i*Y7))/N1), imag(sum(exp(i*Y7))/N1), 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
%latexylabel('Im');
%rotateY;
yleer;
yticks([]);
xticks([-1 0 1]);
hold off
%
%    results of figure 7 (right)
subplot(4,6,[5 6 11 12])
N1=10;
plot(cos(Y8), sin(Y8),'bo');
hold on
plot([0 real(sum(exp(i*Y8))/N1)],...
     [0 imag(sum(exp(i*Y8))/N1)],'b');
plot(real(sum(exp(i*Y8))/N1), imag(sum(exp(i*Y8))/N1), 'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot([-2 2],[0 0],'k:',[0 0],[-2 2],'k:');
axis([-1.1 1.1 -1.1 1.1]);
axis('square');
plot(x, y,'k');  % unit cicle
latexxlabel('Re');
%latexylabel('Im');
%rotateY;
yleer;
yticks([]);
xticks([-1 0 1]);
hold off


%%    Phase locking for in dependence upon coupling strength k
%        KuramotoDiff.m    figure 8
figure(10)        
subplot(4,1,[1 2])
N=10;
Tend=20;
T=0:dt:Tend;
x0=ones(N,1);  %  r_\infty is independent upon x0;
               %  an initial state for quick convergence is used here
omega=2.1:0.2:4;
k2=0;
Coupl=[];
Rinf=[];
kmax=10;
for k=[0.01:0.02:0.8]
    k2=k2+1;
    Coupl(k2)=k;
    sim('KuramotoDiffModel');
    r=[];                    
    k1min=round(0.8*length(Theta));
    for k1=k1min:length(Theta)
        r(k1-k1min+1)=abs(sum(exp(i*Theta(k1,:)))/N);
    end
    Rinf(k2)=sum(r)/length(r);  %  average value of r
end
Tend=150;
T=0:dt:Tend;
for k=[0.9:0.1:2]
    k2=k2+1;
    Coupl(k2)=k;
    sim('KuramotoDiffModel');
    r=[];                    
    k1min=round(0.8*length(Theta));
    for k1=k1min:length(Theta)
        r(k1-k1min+1)=abs(sum(exp(i*Theta(k1,:)))/N);
    end
    Rinf(k2)=sum(r)/length(r);  %  average value of r
end
Tend=30;
T=0:dt:Tend;
for k=[2.1:0.3:kmax]
    k2=k2+1;
    Coupl(k2)=k;
    sim('KuramotoDiffModel');
    r=[];                    
    k1min=round(0.8*length(Theta));
    for k1=k1min:length(Theta)
        r(k1-k1min+1)=abs(sum(exp(i*Theta(k1,:)))/N);
    end
    Rinf(k2)=sum(r)/length(r);  %  average value of r
end
semilogx(Coupl, Rinf, 'b');
latexxlabel('$$k$$');
latexylabel('$$r_\infty\;$$');
rotateY;
axis([0.01, kmax, -0.1, 1.1]);
hold on 
plot([0.01 kmax],[0 0],'k:', [0.01 kmax], [1 1],'k:');
Tend=50;
T=0:dt:Tend;
k2=0;
Coupl=[];
Rinf=[];
for k=[klow, khigh]
    k2=k2+1;
    Coupl(k2)=k;
    sim('KuramotoDiffModel');
    r=[];                    
    k1min=round(0.8*length(Theta));
    for k1=k1min:length(Theta)
        r(k1-k1min+1)=abs(sum(exp(i*Theta(k1,:)))/N);
    end
    Rinf(k2)=sum(r)/length(r);  %  average value of r
end
semilogx(Coupl, Rinf,'bo');
hold on
semilogx([1.2 1.2],[-0.1 1.1], 'k--');
latextext(1.2, -0.3, '$$k_{\rm c}$$');
%title('{\it r}_\infty for different coupling strength');
xticks([0.01 0.1 1 10]);
yticks([0 0.5 1]);
hold off


%%  Phase differences of the three oscillators for different initial states
%     from KuramotoNonSync.m, figure 5 and 6
figure(11)
N=3;
omega=0.5*ones(N,1);
dt=1;
Tend=20;
T=0:dt:Tend;
k=0.3;                    % coupling strength: k/3
x0=[0;
    2/3*pi;
    4/3*pi];
for x2=-pi:0.6:4*pi        % step size 0.5
    for x3=-pi:0.6:4*pi
        x0=[1.0; x2; x3];
        sim('KuramotoDiffModel3TS'); %   3 oscillator outputs
                                     %   in vector Theta
        plot(Theta(:,2)-Theta(:,1), Theta(:,3)-Theta(:,1),'b')
        hold on
    end
end
Stable1=[0 2*pi 4*pi  0   2*pi 4*pi   0   2*pi 4*pi];
Stable2=[0  0    0   2*pi 2*pi 2*pi  4*pi 4*pi 4*pi];
plot(Stable1, Stable2, 'ko');
Zerosum1=[0   0   pi  pi   pi  pi    pi  2*pi 2*pi 3*pi 3*pi 3*pi 3*pi 3*pi 4*pi 4*pi];
Zerosum2=[pi 3*pi  0  pi  2*pi 3*pi 4*pi  pi  3*pi  0    pi  2*pi 3*pi 4*pi  pi  3*pi];
Antiphase1=[2/3*pi 2/3*pi  4/3*pi 4/3*pi 8/3*pi 8/3*pi 10/3*pi 10/3*pi];
Antiphase2=[4/3*pi 10/3*pi 2/3*pi 8/3*pi 4/3*pi 10/3*pi 2/3*pi  8/3*pi];
plot(Antiphase1, Antiphase2, 'k+', Zerosum1, Zerosum2, 'k*');
latexxlabel('$$\theta_2 - \theta_1$$');
latexylabel('$$\theta_3 -\theta_1$$');
rotateY;
latextitle('Synchronisation behaviour of three uniform oscillators')
set(gca, ...
    'XTickLabelMode', 'manual',...
    'XTick', [0 pi 2*pi 3*pi],...
    'TickLabelInterpreter','latex',...
	'XTickLabel',{'$$0$$', '$$\pi$$', '$$2\pi$$', '$$3\pi$$'},...
    'FontName','times',...
	'FontSize',14);
set(gca, ...
    'YTickLabelMode', 'manual',...
    'YTick', [0 pi 2*pi 3*pi],...
    'TickLabelInterpreter','latex',...
	'YTickLabel',{'$$0$$', '$$\pi$$', '$$2\pi$$', '$$3\pi$$'},...
    'FontName','times',...
	'FontSize',14);
axis([-1 11 -1 11]);
axis('square');
hold off

%
%  state space with stable and unstable equilibrium points
figure(12)
plot(Stable1, Stable2, 'ko');
hold on
plot(Antiphase1, Antiphase2, 'k+', Zerosum1, Zerosum2, 'k*');
latexxlabel('$$\theta_2 - \theta_1$$');
latexylabel('$$\theta_3 -\theta_1$$');
rotateY;
latextitle('Synchronisation behaviour of three uniform oscillators')
set(gca, ...
    'XTickLabelMode', 'manual',...
    'XTick', [0 pi 2*pi 3*pi],...
    'TickLabelInterpreter','latex',...
	'XTickLabel',{'$$0$$', '$$\pi$$', '$$2\pi$$', '$$3\pi$$'},...
    'FontName','times',...
	'FontSize',14);
set(gca, ...
    'YTickLabelMode', 'manual',...
    'YTick', [0 pi 2*pi 3*pi],...
    'TickLabelInterpreter','latex',...
	'YTickLabel',{'$$0$$', '$$\pi$$', '$$2\pi$$', '$$3\pi$$'},...
    'FontName','times',...
	'FontSize',14);
axis([-1 11 -1 11]);
axis('square');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%%    Non-uniform oscillators
%       from KuramotoNonSync.m,  figure 4
%
figure(13)
x0=[0;
    2/3*pi;
    4/3*pi];
omega=[0.4;
       0.5;
       0.6];
Tend=50;
dt=0.1;
T=0:dt:Tend;
sim('KuramotoDiffModel3TS');
%   3 oscillator outputs
%   in vector Theta
%
subplot(6,1,[1 2])
plot(T, sin(Theta),'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
%latextitle('Non-uniform Kuramoto oscillators');
latextext(2, 1.4, '$$k=0.3$$');
yticks([-1 0 1]);
hold off;

% Synchronisation error
subplot(6,1,[3 4]);
plot(T, e,'b');
xleer;
latexylabel('$$e_i\;$$');
rotateY;
yticks([-2 0 2]);
axis([0 Tend -2 2]);

%  Average model
subplot(6,1,[5 6])
r=[];
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
yticks([0 1]);
xticks([0 10 20 30 40 50]);
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1 1],'k:');
hold off


%%  Nonlinear couplings
%   Kuramoto oscillators with internal model
%   from KuramotoIMP.m,  figure 6
%
figure(14);
N=10;
k=1.5;       % high coupling coefficient
Tend=3.5;
dt=0.01;      % simulation time step
T=0:dt:Tend;
%    Model
ke=1;
kr=1;
k0=4;
br1=kr*k0;
ke1=ke*k0;
A01=[0 1; 0 0];
B01=[ke1; kr*br1];
C01=[1 0];
D01=0;
A0=A01;      % state-space model of all oscillators
B0=B01;      % in model  KuramotoIMPModel.mdl
C0=C01;      % KuramotoIMPModel.mdl uses also  k, N
D0=D01;
theta0=[    
    3.7378
    1.6475
    3.7878
    4.4687
    1.3933
    0.7378
    1.8641
    2.0029
    2.6651
    3.1910];
omega0=4+1*rand(N,1);  % cf. omegas=1;
xr0=[
    0.9631
    0.5468
    0.5211
    0.2316
    0.4889
    0.6241
    0.6791
    0.3955
    0.3674
    0.9880];
x0=zeros(N*n,1);
for k1=1:N
    x0((k1-1)*n+1)=theta0(k1);
    x0((k1-1)*n+2)=omega0(k1)+xr0(k1);
end
x0=[3.7378
    5.0008
    1.6475
    5.4320
    3.7878
    5.4344
    4.4687
    5.0278
    1.3933
    4.5876
    0.7378
    4.8860
    1.8641
    5.0145
    2.0029
    5.0752
    2.6651
    4.5040
    3.1910
    5.7092];
sim('KuramotoIMPModel');
subplot(6,1,[1 2])
plot(T, sin(Theta),'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
%title('LL sync: Kuramoto osc. with complete communication');
yticks([-1 0 1]);
latextext(0.2, 1.4, '$$k=1.5$$');
hold off;

%  Average model
subplot(6,1,[5 6])
r=[];
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1 1],'k:');
xticks([0 1 2 3 4]);
yticks([0 1]);
hold off
%  alternative representation
subplot(6,1,[3 4]);
plot(T, Theta,'b');
hold on
%text(0.5, 100, '{\it k}=1.5');
%xlabel('{\it t}');
xleer;
latexylabel('$$\theta_i$$');   % in rad
rotateY;
%latextitle('LL sync: Kuramoto osc. with complete communication');
yticks([0 10 20]);
axis([0 Tend 0 max(max(Theta))]);
hold off


%%  The same system with weaker coupling strength
%   weak coupling
%
figure(15)
k=0.2;
Tend=8;
dt=0.01;      % simulation time step
T=0:dt:Tend;
sim('KuramotoIMPModel');
subplot(6,1,[1 2])
plot(T, sin(Theta),'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$y_i\;$$');
rotateY;
%title('LL sync: Kuramoto osc. with complete communication');
yticks([-1 0 1]);
latextext(0.2, 1.4, '$$k=0.2$$');
hold off;

%  Average model
subplot(6,1,[5 6])
r=[];
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1 1],'k:');
xticks([0 2 4 6 8 10]);
yticks([0 1]);
hold off
%  alternative representation
subplot(6,1,[3 4]);
plot(T, Theta,'b');
hold on
%text(0.5, 100, '{\it k}=1.5');
%xlabel('{\it t}');
xleer;
latexylabel('$$\theta_i$$');  % in rad
rotateY;
%latextitle('LL sync: Kuramoto osc. with complete communication');
yticks([0 20 40 60]);
axis([0 Tend 0 max(max(Theta))]);
hold off


%%  Leader-follower sync with star communication topology
%    from KuramotoIMP, figures 12, 13
figure(16)
N=10;
K=zeros(N,N);
for k1=1:N-1
    K(k1+1,1)=1;   % communication from leader to each follower
end
Ni=sum(K');
K=K-diag(Ni);
for k1=1:N
    if Ni(k1)==0
       Ni(k1)=1;
    end
end
k=0.2;
Tend=45;
dt=0.01;      % simulation time step
T=0:dt:Tend;
x0=[1.5708
    1.0000
    1.6475
    4.9284
    3.7878
    5.2866
    4.4687
    5.0268
    1.3933
    4.6758
    0.7378
    5.1139
    1.8641
    5.1247
    2.0029
    5.0418
    2.6651
    5.0768
    3.1910
    5.7427];
sim('KuramotoIMPKij');
subplot(6,1,[1 2])
plot(T, sin(Theta),'b')
axis([0, Tend, -1.1, 1.1]);
hold on
plot([0 Tend],[0 0],'k:');
xleer;
latexylabel('$$\theta_i$$');
rotateY;
%title('LL sync: Kuramoto osc. with complete communication');
yticks([-1 0 1]);
latextext(0.2, 1.4, '$$k=0.2$$');
hold off;
%  Average model
subplot(6,1,[5 6])
r=[];
for k1=1:length(Theta)
    r(k1)=abs(sum(exp(i*Theta(k1,:)))/N);
end
plot(T, r,'b')
latexxlabel('$$t$$');
latexylabel('$$r$$');
rotateY;
axis([0, Tend, -0.1, 1.1]);
hold on 
plot([0 Tend],[0 0],'k:', [0 Tend], [1 1],'k:');
xticks([0 10 20 30 40 60]);
yticks([0 1]);
hold off
%  alternative representation
subplot(6,1,[3 4]);
plot(T, Theta,'b');
hold on
xleer;
latexylabel('$$\theta_i$$ in rad');
rotateY;
%latextitle('LL sync: Kuramoto osc. with complete communication');
yticks([0  150 ]);
axis([0 Tend 0 max(max(Theta))]);
hold off


%% Representation of single outputs
figure(17)
subplot(8,1,[1 2])
plot(T, sin(Theta(:,1)),'b--', T, sin(Theta(:,2)), 'b');
hold on
plot([0 Tend],[0 0],'k:');
axis([0 Tend -1.1 1.1]);
xleer;
latexylabel('$$y_1\;$$');
rotateY;
latextitle('$${\tt SyncKuramotoDemo:}$$ Kuramoto osc. with star communication');
yticks([-1 0 1]);
hold off
subplot(8,1,[3 4])
plot(T, sin(Theta(:,1)),'b--', T, sin(Theta(:,3)), 'b');
hold on
plot([0 Tend],[0 0],'k:');
axis([0 Tend -1.1 1.1]);
xleer;
latexylabel('$$y_2\;$$');
rotateY;
yticks([-1 0 1]);
hold off
subplot(8,1,[5 6])
plot(T, sin(Theta(:,1)),'b--', T, sin(Theta(:,6)), 'b');
hold on
plot([0 Tend],[0 0],'k:');
axis([0 Tend -1.1 1.1]);
xleer;
latexylabel('$$y_5\;$$');
rotateY;
yticks([-1 0 1]);
hold off
subplot(8,1,[7 8])
plot(T, sin(Theta(:,1)),'b--', T, sin(Theta(:,8)), 'b');
hold on
plot([0 Tend],[0 0],'k:');
axis([0 Tend -1.1 1.1]);
latexxlabel('$$t$$');
latexylabel('$$y_7\;$$');
rotateY;
yticks([-1 0 1]);
xticks([0 10 20 30 40 50]);
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   figure(18) after figure(2)
%   figure(19) after figure(3)

%%
%   EPS-figure
%
epsfigc('SyncKuramotoDemo');