%==========================================================================
close all
clear all
clc
format short e

%elementu parametrai:
e1m=0.54;
e2m=0.27;
j1m=18e-3;
j2m=180e-3;
phi1=225;
phi2=-270;
Ri1=1080;
Ri2=270;
Gi1=2.7e-3;
Gi2=0.9e-3;
L1=27e-3;
L2=13.5e-3;
L3=9e-3;
L4=6.75e-3;
C1=72e-9;
C2=36e-9;
C3=24e-9;
R1=22500;

w=1/sqrt(L1*C1)
f=w/(2*pi)
E1m=e1m*exp(j*phi1/180*pi)
E2m=e2m*exp(j*phi2/180*pi)
J1m=j1m*exp(j*phi1/180*pi)
J2m=j2m*exp(j*phi2/180*pi)
ZL3=j*w*L3
ZL12=j*w*(L1+L2)
ZL1=j*w*L1
ZL2=j*w*L2
ZL4=j*w*L4
ZC1=1/(j*w*C1)
ZC2=1/(j*w*C2)
ZC3=1/(j*w*C3)

%sroves ir itampos šaltiniu matrica:
B=[J1m;J2m;-J2m;0;E1m;-E2m;E2m;0;0;0];

%Kintamuju koeficientu matrica:
% iC1   iGp    iRim   iL1 iRin iR1        iL3    iGk   iC2    iL4
A=[ 1,   1,    -1,     1, -1,    1,         0,    0,    0,     0;...
    0,   0,     0,     0,  1,   -1,         1,    1,    0,     0;...
    0,   0,     1,    -1,  0,    0,         0,   -1,   -1,    -1;...
  ZC1,-1/Gi1,   0,     0,  0,    0,         0,    0,    0,     0;...
    0,   0,   Ri1,  ZL12,  0,    0,         0,    0,    0,     0;...
    0,   0,     0, -ZL12,-Ri2,   0,         0,1/Gi2,    0,     0;...
    0,   0,     0,     0, Ri2,  R1,         0,    0,    0,     0;...
    0,1/Gi1,     0,    0,   0, -R1,      -ZL3,    0,    0,     0;...
    0,   0,     0,     0,  0,    0,       ZL3,-1/Gi2,(ZC2+ZC3),0;...
    0,   0,     0,     0,  0,    0,         0,    0,(-ZC2-ZC3),     ZL4];

I=A^-1*B;

%Suskaiciuotos sroves:
%iC1   iGp    iRim   iL1 iRin iR1        iL3    iGk   iC2    iL4
iC1=I(1)
iGp=I(2)
iRim=I(3)
iL1=I(4)
iRin=I(5)
iR1=I(6)
iL3=I(7)
iGk=I(8)
iC2=I(9)
iL4=I(10)
%==========================================================================
%Galiu balansas:
%==========================================================================
SR=0.5*abs(iC1)^2*ZC1+0.5*abs(iGp)^2/Gi1+0.5*abs(iRim)^2*Ri1+...
    0.5*abs(iL1)^2*ZL12+0.5*abs(iRin)^2*Ri2+0.5*abs(iR1)^2*R1+...
    0.5*abs(iL3)^2*ZL3+0.5*abs(iGk)^2/Gi2+0.5*abs(iC2)^2*(ZC2+ZC3)+...
    0.5*abs(iL4)^2*ZL4

SS=0.5*E1m*conj(iRim)+0.5*E2m*conj(iRin)+0.5*iGp/Gi1*conj(J1m)+0.5*iGk/Gi2*conj(J2m)

%==========================================================================
%penkto konturo itampos vektoriškai
%==========================================================================
uZc1=iC1*ZC1
uGp=iGp/Gi1
uRim=iRim*Ri1
uZl1=iL1*ZL1
uZl2=iL1*ZL2
uRin=iRin*Ri2
uR1=iR1*R1
uZl3=iL3*ZL3
uGk=iGk/Gi2
uZc2=iC2*ZC2
uZc3=iC2*ZC3
uZl4=iL4*ZL4

U=[0; -uGp; -uGp+uR1; -uGp+uR1+uZl3];

figure
plot(real(U),imag(U),'-o')
grid on
axis equal
xlabel('Re(U),V')
ylabel('Im(U),V')

%==========================================================================
%penkto konturo itampos laiko srityje
%==========================================================================
T=1/f
t=0:T/100:2*T;
uGp=abs(uGp)*cos(w*t+angle(uGp));
uR1=abs(uR1)*cos(w*t+angle(uR1));
uZl3=abs(uZl3)*cos(w*t+angle(uZl3));
E1m=e1m*cos(w*t+phi1);

figure
plot(t,uGp)
hold on
plot(t,uGp)
plot(t,uR1)
plot(t,uZl3)
grid on
xlabel('t, s')
ylabel('U,V')