%     LCLCNetworkDemo.m
%     from LCLCNetwork.m  of 8.9.2017
%
%  Network of LCLC oscillators
%
%  J. Lunze
%  8.9.2017
%  Version of 1.10.2018
%  for 2nd edition: 29.3.2021
%
echo off
close all
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%
%  Neutrally stable oscillators 
% 
figure(1)
subplot(2,1,1)
%
%   Oscillators
L1=0.1;   % 100 mH
C1=0.00001;   % 10 muF
L2=L1;
C2=C1;
R1=1;   % 1 Ohm
A=[0 1/C1 0 1/C1;
   -1/L1 0 0 0;
   0 0 0 1/C2;
   -1/L2 0 -1/L2 0];
B1=[-1/C1; 0; 0; 0];
B2=[-1/C1; 0; -1/C2; 0];
C1=[-1  0 0 0];
C2=[-1 0 -1 0];
C=[C1; C2];
D=0;
n=4;
Oscillator1=ss(A, B1, C1, D);
Oscillator2=ss(A, B2, C2, D);
Oscillator=ss(A, B2, C, D);
VC01=-1;  % 1 V
IL01=0;  % 0 A
VC02=0;
IL02=0;
x01=[VC01; IL01; VC02; IL02];
Tend=0.015;
dt=0.0001;
Time=[0:dt:Tend];
Y=initial(Oscillator, x01, Time);
plot(Time, Y(:,1), 'b', Time, Y(:,2), 'k--');
latexylabel('$$V_1, V_2$$');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('$${\tt LCLCNetwork:}$$ Oscillator behaviour');
axis([0 Tend -1.4 1.4]);
yticks([-1 0 1]);
xticks([0 0.01]);

%%  Synchronisability analysis
%   with respect to the first channel (u1, y1)
figure(2)
rlocprintblue(Oscillator1);
%rlocfind(Oscillator);
kc1=0.025;
Poles1=rlocus(Oscillator1, kc1);
hold on
plot(real(Poles1), imag(Poles1), 'b^');
latextext(-2700, 1400, ['$$k=$$' num2str(kc1) '$$\;\;\Delta$$' ]);
latextitle('$${\tt LCLCNetwork:}$$ Root locus');
%
kc2=0.04;
Poles2=rlocus(Oscillator1, kc2);
plot(real(Poles2), imag(Poles2), 'b*');
latextext(-2200, 1150, [ num2str(kc2) '$$\;\;\ast$$']);
%
kc3=0.08;
Poles3=rlocus(Oscillator1, kc3);
plot(real(Poles3), imag(Poles3), 'bd');
latextext(-2200, 900, [num2str(kc3) '$$\;\;\diamond\;\;$$']);
plot([0 0], [-2000 2000],'k:');
axis([-3000 400 -1700 1700]);
yticks([-1000 0 1000]);
xticks([-2000 -1000 0]);
hold off


%%  Oscillator network
%  Bidirectional couplings by the resistor R
%  use the first channel (u1, y1)
figure(3);
subplot(4,1,[1 2]);
L=[1 -1  0;
   -1  2 -1;
   0 -1  1 ];    
N1=3;
L1=L*kc1;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L1, B1*C1);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C1);
Network=ss(Abar, Bbar, Cbar, 0);
x0=[x01; zeros((N1-1)*n, 1)];
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('$${\tt LCLCNetwork:}$$ Network behaviour');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xleer
xticksempty([0 0.01]);
latextext(0.0075, 0.8, ['$$1/R=$$' num2str(kc1)]);
%
subplot(4,1,[3 4]);
L3=L*kc3;  % hence, R1=1/kc3 = 1.25 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L3, B1*C1);
Network=ss(Abar, Bbar, Cbar, 0);
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xticks([0 0.01]);
latextext(0.0075, 0.8, ['$$1/R=$$' num2str(kc2)]);



%%  Oscillator network
%  Bidirectional couplings by the resistor R
%  use the first channel (u1, y1)
figure(4);
subplot(4,1,[1 2]);
%  Path structure
N1=10;
kc=kc1;
L=[1 -1  0 0 0 0 0 0 0 0;
   -1 2 -1 0 0 0 0 0 0 0;
   0 -1 2 -1 0 0 0 0 0 0;
   0 0 -1 2 -1 0 0 0 0 0;
   0 0 0 -1 2 -1 0 0 0 0;
   0 0 0 0 -1 2 -1 0 0 0;
   0 0 0 0 0 -1 2 -1 0 0;
   0 0 0 0 0 0 -1 2 -1 0;
   0 0 0 0 0 0 0 -1 2 -1;
   0 0 0 0 0 0 0 0 -1  1    ];    
L3=L*kc;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L3, B1*C1);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C1);
Network=ss(Abar, Bbar, Cbar, 0);
x0=[x01; zeros((N1-1)*n, 1)];
Tend=0.04;
dt=0.0001;
Time=[0:dt:Tend];
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('$${\tt LCLCNetwork:}$$ Network behaviour');
axis([0 Tend -0.5 0.5]);
yticks([-1 -0.5 0 0.5 1]);
xleer
xticksempty([0 0.02 0.04]);
latextext(0.75*Tend, 0.35, 'Path');
latextitle(['$${\tt LCLCNetwork:}$$ Network behaviour for $$1/R=$$' num2str(kc)]);
%
%  Ring structure
subplot(4,1,[3 4])
L=[2 -1  0 0 0 0 0 0 0 -1;
   -1 2 -1 0 0 0 0 0 0 0;
   0 -1 2 -1 0 0 0 0 0 0;
   0 0 -1 2 -1 0 0 0 0 0;
   0 0 0 -1 2 -1 0 0 0 0;
   0 0 0 0 -1 2 -1 0 0 0;
   0 0 0 0 0 -1 2 -1 0 0;
   0 0 0 0 0 0 -1 2 -1 0;
   0 0 0 0 0 0 0 -1 2 -1;
   -1 0 0 0 0 0 0 0 -1 2    ];  
L4=L*kc;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L4, B1*C1);
Network=ss(Abar, Bbar, Cbar, 0);
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
axis([0 Tend -0.5 0.5]);
yticks([-1 -0.5 0 0.5 1]);
xticks([0 0.02 0.04]);
latextext(0.75*Tend, 0.35, 'Ring');
hold off

%%  Synchronisability analysis
%   with respect to the second channel (u2, y2)
figure(5)
rlocprintblue(Oscillator2);
%rlocfind(Oscillator);
% kc1=0.15;
% Poles1=rlocus(Oscillator2, kc1);
% hold on
% plot(real(Poles1), imag(Poles1), 'k^');
% latextext(-90, 50, ['$$\Delta\;\;$$' num2str(kc1)]);
latextitle('$${\tt LCLCNetwork:}$$ Root locus');
%
kc2=0.015;
Poles2=rlocus(Oscillator2, kc2);
plot(real(Poles2), imag(Poles2), 'b*');
latextext(-2700, 1400, ['$$k=$$' num2str(kc2) '$$\;\;\ast$$' ]);
%
% kc3=0.8;
% Poles3=rlocus(Oscillator2, kc3);
% plot(real(Poles3), imag(Poles3), 'kd');
% latextext(-90, 30, ['$$\diamond\;\;$$' num2str(kc3)]);
plot([0 0], [-2000 2000],'k:');
axis([-3000 400 -1700 1700]);
yticks([-1000 0 1000]);
xticks([-2000 -1000 0]);
hold off


%%  Oscillator network
%  Bidirectional couplings by the resistor R
%  use the first channel (u1, y1)
figure(6);
subplot(2,1,1);
L=[1 -1  0;
   -1  2 -1;
   0 -1  1 ];    
N1=3;
L1=L*kc2;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L1, B2*C2);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C2);
Network=ss(Abar, Bbar, Cbar, 0);
x0=[x01; zeros((N1-1)*n, 1)];
Tend=0.05;
dt=0.0001;
Time=[0:dt:Tend];
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
xticks([0 0.02 0.04 0.06]);
latextitle('$${\tt LCLCNetwork:}$$ Network behaviour');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
latextext(0.035, 0.8, ['$$1/R=$$' num2str(kc2)]);


%% 
epsfigc15('LCLCNetworkDemo');

