%     LCNetworkDemo.m
%     from LCNetwork.m  of 8.9.2017
%
%  Network of LC oscillators
%
%  J. Lunze
%  7.9.2017
%  Version of 23.12.2018
%  for 2nd edition: 30.01.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
R1=1;   % 1 Ohm
A=[0 1/C1;
   -1/L1 0];
B1=[-1/C1; 0];
C1=[-1  0];
D=0;
n=2;
Oscillator=ss(A, B1, C1, D);
VC01=-1;  % 1 V
IL01=0;  % 0 A
x01=[VC01; IL01];
Tend=0.015;
dt=0.0001;
Time=[0:dt:Tend];
Y=initial(Oscillator, x01, Time);
plot(Time, Y, 'b');
latexylabel('$$V$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('$${\tt LCNetwork:}$$ Oscillator behaviour');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xticks([0 0.01]);

%%  Synchronisability analysis
%
figure(2)
rlocprintblue(Oscillator);
axis([-500 50 -275 275]);
yticks([-250 0 250]);
xticks([-500 -250 0]);
%rlocfind(Oscillator);
kc1=0.005;
Poles1=rlocus(Oscillator, kc1);
hold on
plot(real(Poles1), imag(Poles1), 'b^');
latextext(-400, 240, ['$$\Delta\;\;$$' num2str(kc1)]);
latextitle('$${\tt LCNetwork:}$$ Root locus');
%
kc2=0.04;
Poles2=rlocus(Oscillator, kc2);
plot(real(Poles2), imag(Poles2), 'b*');
latextext(-400, 200, ['$$\ast\;\;$$' num2str(kc2)]);
%
kc3=0.08;
Poles3=rlocus(Oscillator, kc3);
plot(real(Poles3), imag(Poles3), 'bd');
latextext(-400, 160, ['$$\diamond\;\;$$' num2str(kc3)]);
hold off


%%  Oscillator network
%  Bidirectional couplings by the resistor R1
figure(3);
subplot(2,1,1);
L=[1 -1  0;
   -1  2 -1;
   0 -1  1 ];    
N1=3;
L2=L*kc1;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L2, 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 LCNetwork:}$$ Network behaviour');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xleer
latextext(0.01, 0.8, ['$$1/R=$$' num2str(kc1)]);
%
subplot(2,1,2);
L3=L*kc2;  % hence, R1=1/kc2 = 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.01, 0.8, ['$$1/R=$$' num2str(kc2)]);


%%  Analysis of the coupling effects
%
figure(4)
%  rlocprint.m  of 6.3.2017
[R, K]=rlocus(Oscillator);
for i=1:n
    plot(real(R(i,:)), imag(R(i,:)),'b');
    hold on
end
plot([-100 100],[0 0],'k:', [0 0], [-300 300],'k:');
[P, Z]=pzmap(Oscillator);
for i=1:length(P)
    plot(real(P(:)), imag(P(:)), 'kx');
end
for i=1:length(Z)
    plot(real(Z(:)), imag(Z(:)), 'ko');
end
axis('square');

latexxlabel('Re');
latexylabel('Im');
rotateY;
%
%  First network
LambdaL=eig(L2);
LambdaL=sort(LambdaL);
kc=LambdaL(2);
Poles1=rlocus(Oscillator, kc);
plot(real(Poles1), imag(Poles1), 'b^');
kc=LambdaL(3);
Poles2=rlocus(Oscillator, kc);
plot(real(Poles2), imag(Poles2), 'b^');
latextext(-400, 200, ['$$\Delta\;1/R=$$' num2str(kc1)]);
%   Second network
LambdaL=eig(L3);
LambdaL=sort(LambdaL);
kc=LambdaL(2);
Poles1=rlocus(Oscillator, kc);
hold on
plot(real(Poles1), imag(Poles1), 'b*');
kc=LambdaL(3);
Poles2=rlocus(Oscillator, kc);
plot(real(Poles2), imag(Poles2), 'b*');
axis([-500 50 -275 275]);
yticks([-250 0 250]);
xticks([-500 -250 0]);
latextext(-400, 160, ['$$\ast\;\;1/R=$$' num2str(kc2)]);
latextitle('$${\tt LCNetwork:}$$ Network eigenvalues');
hold off



%%  Oscillator network
%  Bidirectional couplings by the resistor R1
figure(5);
subplot(2,1,1);
%  Path structure
N1=10;
kc=kc2;
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    ];    
L4=L*kc;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L4, 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=1;
%dt=0.005;
%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 LCNetwork:}$$ Network behaviour for $$1/R=$$' num2str(kc)]);
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xleer
latextext(0.75*Tend, 0.8, 'Path');
%
%  Ring structure
subplot(2,1,2)
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    ];  
L5=L*kc;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L5, 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.75*Tend, 0.8, 'Ring');


%%  Analysis of the coupling effects
%
figure(6)
subplot(2,2,1);
%  rlocprint.m  of 6.3.2017
[R, K]=rlocus(Oscillator);
for i=1:n
    plot(real(R(i,:)), imag(R(i,:)),'b');
    hold on
end
plot([-1000 1000],[0 0],'k:', [0 0], [-1000 1000],'k:');
% [P, Z]=pzmap(Oscillator);
% for i=1:length(P)
%     plot(real(P(:)), imag(P(:)), 'kx');
% end
% for i=1:length(Z)
%     plot(real(Z(:)), imag(Z(:)), 'ko');
% end
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
%
%  First network
LambdaL=eig(L4);
LambdaL=sort(LambdaL);
for k1=2:N1
    kc=LambdaL(k1);
    Poles=rlocus(Oscillator, kc);
    plot(real(Poles), imag(Poles), 'b^');
end
latextext(-400, 200, 'Path');
axis([-500 50 -275 275]);
yticks([-250 0 250]);
xticks([-500 -250 0]);
hold off
%
%
figure(7)
subplot(2,2,1)
%  rlocprint.m  of 6.3.2017
[R, K]=rlocus(Oscillator);
for i=1:n
    plot(real(R(i,:)), imag(R(i,:)),'b');
    hold on
end
plot([-1000 1000],[0 0],'k:', [0 0], [-1000 1000],'k:');
% [P, Z]=pzmap(Oscillator);
% for i=1:length(P)
%     plot(real(P(:)), imag(P(:)), 'kx');
% end
% for i=1:length(Z)
%     plot(real(Z(:)), imag(Z(:)), 'ko');
% end
axis('square');
latexxlabel('Re');
latexylabel('Im');
rotateY;
%  Second network
LambdaL=eig(L5);
LambdaL=sort(LambdaL);
for k1=2:N1
    kc=LambdaL(k1);
    Poles=rlocus(Oscillator, kc);
    plot(real(Poles), imag(Poles), 'b*');
end
latextext(-400, 200, 'Ring');
axis([-500 50 -275 275]);
yticks([-250 0 250]);
xticks([-500 -250 0]);
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Oscillator network
%  Bidirectional couplings by the resistor R1
%     from OscillatorNetwork.m of 21.06.2020
figure(8);
subplot(4,1,[1 2]);
L=[1 -1;
  -1  1 ];    
N1=2;
L0=0.1;   % 100 mH
C0=0.00001;   % 10 muF
Ltilde=L*kc1;  % hence, R1=1/kc = 10 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(Ltilde, 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,Time,X]=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_{{\rm C}i}$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('$${\tt LCNetworkDemo:}$$ Network behaviour');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xleer
latextext(0.01, 0.8, ['$$1/R=$$' num2str(kc1)]);
%
% Lyapunov function
V=[];
for k1=1:length(X)
    V(k1)=0.5*C0*X(k1,1)*X(k1,1)+0.5*L0*X(k1,2)*X(k1,2)+...
        0.5*C0*X(k1,3)*X(k1,3)+0.5*L0*X(k1,4)*X(k1,4);
end
subplot(4,1,[3 4])
plot(Time, V, 'b');
axis([0 Tend 0 1.1]);
latexylabel('$$v$$');
rotateY;
latexxlabel('$$t$$ in s');
yticks([0 0.000005]);
xticks([0 0.01]);
%  Final value
Vinf=0.25*C0*(x0(1)+x0(3))^2+0.25*L0*(x0(2)+x0(4))^2;
hold on
plot([0 Tend],[Vinf Vinf],'b--');
latextext(0.001,0.0000017, '$$v_\infty$$');
axis([0 Tend 0 0.000006]);
hold off



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Oscillator network
%  Bidirectional couplings by the resistor R1
%     similar to OscillatorNetwork.m of 11.7.2020, figure 5
figure(9);
%   Oscillators
L01=0.1;   % 100 mH
C01=0.00001;   % 10 muF
R1=1;   % 1 Ohm
A1=[0 1/C01;
   -1/L01 0];
B1=[-1/C01; 0];
C1=[-1  0];
D=0;
n=2;
%  Oscillator 2 with different parameters
L02=0.20;   % 0.2 = 200 mH
C02=C01;  
A2=[0 1/C02;
   -1/L02 0];
B2=[-1/C02; 0];
C2=[-1  0];
D2=0;
%  Bidirectional couplings by the resistor R1
Tend2=0.03;
dt=0.0001;
Time2=[0:dt:Tend2];
%
N1=3;
L=[1 -1  0;
   -1  2 -1;
   0 -1  1 ];
L=L*0.6*kc1;  % hence, R1=1/kc = 10 Ohm
%
Abar=[A1 zeros(n,n) zeros(n,n); 
      zeros(n,n) A2 zeros(n,n);
      zeros(n,n) zeros(n,n) A1]-...
     [B1 zeros(n,1) zeros(n,1); 
     zeros(n,1) B2 zeros(n,1);
     zeros(n,1) zeros(n,1) B1]*L*...
     [C1 zeros(1,n) zeros(1,n); 
     zeros(1,n) C2 zeros(1,n);
     zeros(1,n) zeros(1,n) C1];
Bbar=zeros(N1*n, 1);
Cbar=[C1 zeros(1,n) zeros(1,n); 
     zeros(1,n) C2 zeros(1,n);
     zeros(1,n) zeros(1,n) C1];
Network=ss(Abar, Bbar, Cbar, 0);
%x0=[x01; zeros((N1-1)*n, 1)];
x0=[1 0 0.5 0 -1 0];
[Y,Time2,X]=initial(Network, x0, Time2);
%
subplot(4,1,[1,2]);
plot(Time2, Y, 'b');
latexylabel('$$V_{{\rm C}i}$$ in V');
rotateY;
latextitle('$${\tt LCNetworkDemo:}$$ Network behaviour');
axis([0 Tend2 -1.1 1.1]);
yticks([-1 0 1]);
xleer
%
% Lyapunov function
V=[];
for k1=1:length(X)
    V(k1)=0.5*C0*X(k1,1)*X(k1,1)+0.5*L0*X(k1,2)*X(k1,2)+...
        0.5*C02*X(k1,3)*X(k1,3)+0.5*L02*X(k1,4)*X(k1,4);
end
subplot(4,1,[3,4])
plot(Time2, V, 'b');
axis([0 Tend2 0 0.000006]);
latexylabel('$$v$$');
rotateY;
latexxlabel('$$t$$ in s');
yticks([0 0.000005 0.00001]);
xticks([0 0.01 0.02 0.03 0.04 0.05 0.06]);
%


%% 
epsfigc15('LCNetworkDemo');

