clear
ts=1            %time when is rotor bar breaken
a=5             %break bar - parametr
 

f=50            %Frequency Hz
P=1             %pole number
Rs=7.58         %Ohm
Rb=71.5e-5      %Ohm
Rbb=5           %Ohm
Lb=0.1e-6       %H
Re=1.5e-6       %Ohm
Ls=0,4612       %H
L1s=Ls          %H
Nr=16           %Number rotor bars
Le=1e-7         %H end ring segment
r=70e-3         %m
Ns=160          %
g=0.75e-3       %m
l=65e-3         %m lenght of rotor stack
J=0.027         %kg.m2
Us=230*sqrt(2)  %V
TL1=5
TL2=20
we=2*pi*f
wk=2*pi*f       %rad/s


Ks=sqrt(2/3)
Kr=sqrt(2/Nr)
Ksk=2/3
Krk=2/3*sqrt(Nr/3)
Ur=zeros(Nr,1)

mo=4*pi*1e-7
ar=(2*pi/Nr)
delta=ar/2
Lms=mo*l*r*Ns*Ns*pi/(4*P*P*g)
Lm=Lms*4*sin(P*ar/2)/(pi*Ns)

Rr=2*(Re+Rb*(1-cos(ar)))
Lr=(mo*l*r*ar/g+2*Le+2*Lb*(1-cos(ar)))

Inv_Beta=zeros(Nr,1)
Beta=zeros(1,Nr)
Pomocny=zeros(Nr,1)
for i=1:Nr
    Inv_Beta(i,1)=exp(-j*ar/2*(i-1))
    
    Beta(1,i)=exp(j*ar/2*(i-1))
    Pomocny(i,1)=i
end
Re_Beta=real(Beta)
Beta_k_1_1=exp(j*ar/2*(1))



% Resistance Matrix - only with fault
% k-th bar with fault
k1=2    %k1-th bar with fault
Const_dRr1=Krk*Kr*Beta(1,k1-1)*(1-Beta(1,2)) %*(Rbb-Rb)
Beta_k1_1=exp(j*ar/2*(k1-1))

k2=5    %k2-th bar with fault

Const_dRr2=Krk*Kr*Beta(1,k2-1)*(1-Beta(1,2)) %*(Rbb-Rb)
Beta_k2_1=exp(j*ar/2*(k2-1))
Beta_k2_2=exp(j*ar/2*(k2-2))


Const_Tem=Lm/2*1/Kr/Ks*sqrt(Nr/3)



%..........................
% TRansformation
%...........................
Rk=zeros(4,4)
Rk(1,1)=Rs
Rk(2,2)=Rs
Rk(3,3)=2*(Re+Rb*(1-cos(ar)))*sqrt(3/Nr)
Rk(4,4)=Rk(3,3)


Lk=zeros(4,4)
Lk(1,1)=L1s+3/2*Lms
Lk(2,2)=Lk(1,1)
Lk(1,3)=Ks/Kr*3/2*sqrt(3/Nr)*Lm
Lk(2,4)=Lk(1,3)
Lk(3,1)=Nr/2*sqrt(Nr/3)*Kr/Ks*Lm
Lk(4,2)=Lk(3,1)
Lk(3,3)=Lr
Lk(4,4)=Lk(3,3)
Inv_Lk=inv(Lk)

Lwk_12=-(L1s+3/2*Lms)
Lwk_21=-Lwk_12
Lwk_14=-3/2*Ks/Kr*sqrt(3/Nr)*Lm
Lwk_23=-Lwk_14
Lwk_32=-Nr/2*sqrt(Nr/3)*Kr/Ks*Lm
Lwk_41=-Lwk_32
Lwk_34=-Lr
Lwk_43=-Lwk_34




