clc

% Zakladni technicke parametry transformatoru
Un1=110
Un2=22
Snt=10
dP0p=0.3
dPfep=dP0p
dPkp=1
dPcup=dPkp
Ukp=10
I0p=1
frekv=50
omega=2*pi*frekv

% Vypocty parametru nahraniho schematu
Znt=Un1^2/Snt

Rk=(dPkp/100)*(Un1^2/Snt)
Rk1=Rk/2
Zk=(Ukp/100)*(Un1^2/Snt)
Xs=sqrt(Zk^2-Rk^2)
Ls=Xs/omega
Ls1=Ls/2

Gfe=(dP0p/100)*(Snt/Un1^2)
Rfe=1/Gfe
Y0=(I0p/100)*(Snt/Un1^2)
Xh=1/sqrt(Y0^2-Gfe^2)
Lh=Xh/omega

% Kontrola pres komplexni vyjadreni impedance Zk a admitance Y0
Zk2=abs(Rk+i*Xs)
Y02=abs(1/Rfe + 1/(i*Xh))

% Vypocet chodu naprazdno
Tau0=Lh/Rk1
Z0=sqrt(Rk1^2+Xh^2)
I0m=Un1/sqrt(3)*sqrt(2)/Z0
Psi0=atan(Xh/Rk1)
% Rozdil faze oproti ryze induktivnímu proudu
OdchylkaPsi0=Psi0-pi/2
% Integracni konstanta
C0=I0m*sin(Psi0)

% Kontrola pres komplexni vyjadreni I0
I0c=Un1/sqrt(3)/(Rk1+i*Xh)
% Amplituda alternativne
I0m2=sqrt(2)*abs(I0c)
% Faze alternativne 
Psi02=angle(I0c)

t=[0:0.001:0.1]';
I0t=C0*exp(-t/Tau0) + I0m*sin(omega*t-Psi0);
plot(t,I0t);

% Vypocet chodu nakratko
TauK=Ls/Rk
Ikm=Un1*Ukp/100/sqrt(3)*sqrt(2)/Zk
PsiK=atan(Xs/Rk)
% Rozdil faze oproti ryze induktivnímu proudu
OdchylkaPsiK=PsiK-pi/2
% Integracni konstanta
Ck=Ikm*sin(PsiK)

% Kontrola pres komplexni vyjadreni Ik
Ikc=Un1*Ukp/100/sqrt(3)/(Rk+i*Xs)
% Amplituda alternativne
Ikm2=sqrt(2)*abs(Ikc)
% Faze alternativne 
PsiK2=angle(Ikc)

% Kontrola pres jmenovity proud
In=Snt/Un1/sqrt(3)
Ikm3=In*sqrt(2)
I0m3=I0p/100*In*sqrt(2)

figure;
Ikt=Ck*exp(-t/TauK) + Ikm*sin(omega*t-PsiK);
plot(t,Ikt);

% Reseni prubehu proudu numerickymi metodami
tspan = 0 : 0.0001 : 0.1;

% Reseni prubehu proudu naprazdno numerickou metodou - v1
y0 = 0;
[t,y] = ode23(@f1, tspan, y0,[]);
figure;
plot(t,y);

% Reseni prubehu proudu nakratko numerickou metodou - v1
y0 = 0;
[t,y] = ode23(@f1k, tspan, y0,[]);
figure;
plot(t,y);

% Reseni prubehu proudu nakratko numerickou metodou - v2
y0 = [0; 0];
[t,y] = ode23(@f2, tspan, y0,[]);
figure;
plot(t,y);

% Odhad rozdilu mezi proudy nakratko na pocatku prechodneho deje
DeltaIk=In*I0p/100*Ukp/100/2*sqrt(2)*1.8
