How to add a fault to an asynchronous machin

can u help me to solve the problem in this code and how it run.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dy = file(t,y)

global rs rr lp lm lf j f v teta0 p n0 alfa cr rin s b0

%**********************************************************************

if t<1

cr=0;

n0=0.0;

s=0;

b0=0;

end

if t>=1&&t<20

cr=10;

n0=0 ;

s=0;

b0=0;

end

if t>=2&&t<20

n0=0 ;

s=0;

b0=1;

end

%**********************************************************************

% parametres de la MAS

rs=2.86; % résistance propre d'une phase statorique.

rr=2.756; % résistance propre d'une phase rotorique.

lp=0.3953;

lm=(1.5)*lp ;% inductance magnétisante.

lf=0.009594 ;% inductance de fuite totaliser au stator.

j=0.023976; %

f=0.0024439;%

v=380;

teta0=(0);

p=1;

alfa=(2/3)*n0;

%**********************************************************************

% les équations de tension.

va = v*sin(50*2*pi*t);

vb = v*sin((50*2*pi*t)-(2*pi/3));

vc = v*sin((50*2*pi*t)+(2*pi/3));

%**********************************************************************

% les équations de tension dans le repère de Pa.

uds=va*cos(y(6)) +(1/sqrt(3))* (vb-vc)*sin(y(6));

uqs=-va*sin(y(6)) +(1/sqrt(3))* (vb-vc)*cos(y(6));

cosl=(1/pi).*sin(pi*0.029).*cos(2*pi*140*t)+(1/2*pi)*sin(2*pi*0.029).*cos(4*pi*140*t)+(1/3*pi)*sin(3*pi*0.029).*cos(6*pi*140*t)+(1/4*pi)*sin(9*pi*0.029).*cos(8*pi*140*t)+(1/5*pi)*sin(9*pi*0.029).*cos(18*pi*140*t)+(1/6*pi)*sin(6*pi*0.048).*cos(14*pi*140*t)+(1/7*pi)*sin(7*pi*0.048).*cos(14*pi*140*t)+(1/8*pi)*sin(8*pi*0.048).*cos(16*pi*140*t);

%**********************************************************************

% le courant de la phase a.

i=(y(1)*cos(y(6)) + y(2)*cos(y(6)+(pi/2)));

%**********************************************************************

% le couple Cem

%Cem=p*(y(2).*y(3)-y(1).*y(4));

% le glissement g

g=((100*pi)-(y(5).*p))/(100*pi);

%**********************************************************************% le système d'état.

dy =[(-(rs/lf)-((rr/lf)*(1+(alfa*(sin(teta0))^2))/(1+alfa)))*(y(1))+y(5)*(y(2))+((rr/lf*lm)*((1+(alfa*(sin(teta0))^2))/(1+alfa)))*y(3)+(1/lf)*y(5)*y(4)+(1/lf)*uds;

-y(5)*y(1)-((rs/lf)+((rr/lf)*(1+(alfa*(cos(teta0))^2))/(1+alfa)))*(y(2))-(1/lf)*y(5)*y(3)+((rr/lf*lm)*((1+(alfa*(cos(teta0))^2))/(1+alfa)))*y(4)+(1/lf)*uqs;

(rr*((1+(alfa*(sin(teta0))^2))/(1+alfa)))*(y(1))-((rr/lm)*((1+(alfa*(sin(teta0))^2))/(1+alfa)))*y(3);

+(rr*((1+(alfa*(cos(teta0))^2))/(1+alfa)))*(y(2))+0*(((rr/2*lm)*((alfa*sin(2*teta0))/(1+alfa))))*y(3)-((rr/lm)*((1+(alfa*(cos(teta0))^2))/(1+alfa)))*y(4);

((p^2)/j)*(y(2))*y(3)-((p^2)/j)*y(1)*y(4)-(p/j)*(cr+s* cosl)-(f/j)*y(5);

y(5)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PARTIE A : programme à sauvegardé

clc

global rs rr lp lm lf j f v teta0 p n0 alfa cr f0 frt N fc m0

%**********************************************************************

tacq = 20;

Nb=18;

Nt=2^Nb

Te =tacq/(Nt-1);

fs=1/Te

tn = 0 :Te : tacq-Te ;

[t,y] = ode23tb('file',0:Te:tacq,[0 0 0 0 0 0]);

%**********************************************************************

%le curant de la phase a.

i=(y(:,1).*cos(y(:,6))+y(:,2).*cos(y(:,6)+(pi/2)));

% le couple

Cem=p*(y(:,2).*y(:,3)-y(:,1).*y(:,4));

% le glissement g

g=((100*pi)-(y(:,5)))/(100*pi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% résul

figure(1)

plot(t,y(:,1))

legend('Ids')

grid on

figure(2)

plot(t,y(:,2))

legend('Iqs')

grid on

figure(3)

plot(t,y(:,3))

legend('fdr')

grid on

figure(4)

plot(t,y(:,4))

axis([1.5 3 -4 4])

legend('fqr')

grid on

figure(5)

plot(t,y(:,5))

%axis([1.999 2.299 290 300])

grid on

ylabel('W (rad / s)')

xlabel('temps (s)')

% figure(6)

% plot(t,y(:,6))

% legend('teta')

% grid on

figure(7)

plot(t,i);

axis([1.2 3 -30 30]);

grid on

ylabel('Ia(A)')

xlabel('temps (s)')

figure(8)

plot(t,Cem)

% axis([1.8 6 7 11])

grid on

ylabel('Cem (N.m)')

xlabel('temps (s)')

% figure(9)

% plot(t,g)

% grid on

% ylabel('g')

% xlabel('temps (s)')

v=380*cos(2*pi*50*t);

Pins=v.*i;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% conv e r s ion +/? 2V avec Nb i t s e t +/? 1/2LSB de non l i n e a r i t e

Ucan = 4 ; Nbits = 12 ;

xn = Ucan*round( (i./Ucan )*2.^( Nbits -1) )/2.^( Nbits -1);

% c a l c u l des s p e c t r e s

Nfft = 2.^ ceil ( log2 ( length ( xn ) ) )

% ******************************

% fs=5000; % the sampling rate in samples/s

% Te=1/fs; % the sampling period

N=length(xn); % Record length of data

L=2500; % length of each segment

M=fix(N/L); % the number of segments to be averaged

K=Nt; % number of points in the frequency domain

w=hamming(L); % w[n], the Hamming window of length L

% ****************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[P,f]=psd(xn,K,fs,0); % welch periodogram with no overlap

P=((2*Te)/(L))*P; % apply the appropriate scaling units for W/Hz

P=20*log10(P); % convert to dB

figure(111) %generate a Figure window

plot(f,P);grid on %plot the estimated power spectrum P as fcn of freq f

axis([0 200 -350 30])

xlabel('frequence (Hz)') % label the frequency axis

ylabel('Amplitude (dB)') % label the amplitude axis

title(' Estimation de DSP du courant statorique ')

legend( ' la DSP du courant' )

% % NFFT = 2^nextpow2(L); % Next power of 2 from length of y

%

% Y =( fft(xn,N))/L;

% f = fs/2*linspace(0,1,N/2);

% % Plot single-sided amplitude spectrum.

% figure(12)

% plot(f,20*log10((Y(1:N/2))));grid on

% axis([0 800 -70 50])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% c a l c u l des s p e c t r e s

%Nfft = 2.^ ceil ( log2 ( length ( xn ) ) );

% f e n e t r e s r e c t a n g u l a i r e e t de Blackman

wr = ones ( size ( xn ) ) ;

wb = ( blackman ( length ( xn ) ) ) ;

wn =(hanning(length ( xn )));

wm=(hamming(length ( xn )));

xnwr = xn .*wr ;

xnwb = xn .*wb ;

xnwm= xn.*wm;

xnwn=xn.*wn;

% a j out de z e ros

xnwr = [ xnwr , zeros( 1 , Nfft-length( xnwr ) ) ] ;

xnwb = [ xnwb , zeros( 1 , Nfft-length(xnwb ) ) ] ;

xnwh = [ xnwm , zeros( 1 , Nfft-length(xnwm ) ) ] ;

% f f t

xjr = (fft ( xnwr )/ length ( xn )) ;

xjb = (fft (xnwb))/ length ( xn ) ;

xjm= (fft ( xnwm )/ length ( xn ));

xjn=(fft ( xnwn )/ length ( xn ));

% domaine s p e c t r a l

fmax = 1/Te ;

df = fmax/Nfft ;

ff = 0 : df : fmax-df ;

% i n f o s

%Nbits , tacq , Te , fmax , df

% Pac = var ( xn )

% Npoints = round( tacq /Te ) , Nfft

% % figure(1) ;

% subplot ( 3 , 1 , 1 ) ;

%

% xlabel ( ' temps [ s e c ] ' ) ;

%printdepsansptemps.eps

% s p e c t r e s

% figure(2) ; % f e n e t r e r e c t a n g u l a i r e

% subplot ( 2 , 1 , 1 ) ;

% zoom s p e c t r a l

fz1 = 0; fz2 = 200; % domaine i n t e r e s s a n t

dbmax = 80 ;

figure(13) ;

% subplot ( 2 , 1 , 1 ) ;

plot ( ff , 20*log10( abs(xjr ) ) ) ; hold on ;

axis ( [ fz1 , fz2 ,-2*dbmax ,50 ] ) ; grid ;

% title ( texte ) ;

ylabel ( 'amplitude [ dB] ' ) ;

legend( ' Fenêt re r e c t a n g u l a i r e ’' ) ;

% subplot ( 2 , 1 , 2 ) ;

figure(14)

plot ( ff , 20*log10( (xjb ) ) ) ;

axis ( [ fz1 , fz2 ,-2.5*dbmax , 50 ] ) ; grid ;

ylabel ( 'amplitude [ dB] ' ) ;

xlabel ( ' f r équenc e [Hz ] ' ) ;

legend( ' Fenêt re de Blackman ' ) ;

title (' spactre de courant statorique') ;

figure(15)

plot ( ff , 20*log10( (xjn ) ) ) ;

axis ( [ fz1 , fz2 ,-2.5*dbmax , 50 ] ) ; grid ;

ylabel ( 'amplitude [ dB] ' ) ;

xlabel ( ' f r équenc e [Hz ] ' ) ;

legend( ' Fenêt re de hanning ' ) ;

title (' spactre de courant statorique') ;

figure(16)

plot ( ff , 20*log10( (xjm) ) ) ;

axis ( [ fz1 , fz2 ,-2.5*dbmax , 50 ] ) ; grid ;

ylabel ( 'amplitude [ dB] ' ) ;

xlabel ( ' f r équenc e [Hz ] ' ) ;

legend( ' Fenêt re de hamming ' ) ;

title (' spactre de courant statorique') ;

figure(17)

plot(ff , 20*log10( (xjn ) ),ff , 20*log10( (xjb ) ))

axis ( [ fz1 , fz2 ,-2.5*dbmax , 50 ] ) ; grid ;

legend( ' Fenêt re de hanning ',' Fenêt re de Blackman ' )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%