// Programme P2_2
// 06-04-2010
//----------------
s=%s;
z=%z;

Te=0.0003;						// Sample time
Fmin=1;							//Hertz
Fmax=1000;						// Hertz
//--------------------------------------------------------------------------------------------------------------------------
K=0.43;
T=0.0075;
a=0.2;
//--------------------------------------------------------------------------------------------------------------------------
FT0=K*(1+T*s)/(1+a*T*s);				// Continuous system
FT0 = syslin('c', FT0);
FT10_0 = tf2ss(FT0);

//--------------------------------------------------------------------------------------------------------------------------
pz=(1-z^-1)/Te;						//  Euler Transform

FT1 = K*(1+T*pz)/(1+a*T*pz);			
FT1 = syslin(Te, FT1 );
FT10_1 = tf2ss(FT1);

//--------------------------------------------------------------------------------------------------------------------------
pz=2/Te*(1-z^-1)/(1+z^-1);				// Bilinear transform

FT2 = K*(1+T*pz)/(1+a*T*pz);
FT2 = syslin(Te, FT2 );
FT10_2 = tf2ss(FT2);

//--------------------------------------------------------------------------------------------------------------------------
b=1/(a*T);						// Invariance impulsionnelle
c=1/a;
EXP=exp(-b*Te);
N=c*(z-1) + 1 - EXP;
D=z-EXP;
FT3=K*(N/D);
FT3 = syslin(Te, FT3 );
FT10_3 = tf2ss(FT3);

//--------------------------------------------------------------------------------------------------------------------------
scf(1);clf();					// Create and clear screen 1
bode([FT10_0],Fmin,Fmax);
legends(['Continu';], [1;], with_box=%t, opt=4 );
//--------------------------------------------------------------------------------------------------------------------------
scf(2);clf();					// Create and clear screen 2
bode([FT10_0;FT10_3],Fmin,Fmax);
legends(['Continu';'Invariance'], [1;2], with_box=%t, opt=4 );

//--------------------------------------------------------------------------------------------------------------------------
scf(3);clf();					// Create and clear screen 3
bode([FT10_0;FT10_1;FT10_2;FT10_3],Fmin,Fmax);
legends(['Continu';'Euler';'Bilineaire';'Invariance';], [1;2;3;4], with_box=%t, opt=4 );
//--------------------------------------------------------------------------------------------------------------------------
frq=Fmin:0.01:Fmax;

[frq0,rep0] =repfreq(FT10_0,frq);
[db0,phi0]=dbphi(rep0);

[frq1,rep1] =repfreq(FT10_1,frq);
[db1,phi1]=dbphi(rep1);	

[frq2,rep2] =repfreq(FT10_2,frq);
[db2,phi2]=dbphi(rep2);	

[frq3,rep3] =repfreq(FT10_3,frq);
[db3,phi3]=dbphi(rep3);	

//Affichage
//------------
scf(4);clf();					// Create and clear screen 4
xgrid();
subplot(2,1,1);
plot2d('ln',frq,db0,style=1);
plot2d('ln',frq,db1,style=2);
plot2d('ln',frq,db2,style=3);
plot2d('ln',frq,db3,style=4);
legends(['Continu';'Euler';'Bilineaire';'Invariance';], [1;2;3;4], with_box=%t, opt=4 );
xtitle('', 'Fréquence Hz','Magnitude (Db)');

subplot(2,1,2);
plot2d('ln',frq,phi0,style=1);
plot2d('ln',frq,phi1,style=2);
plot2d('ln',frq,phi2,style=3);
plot2d('ln',frq,phi3,style=4);
xgrid();
xtitle('', 'Fréquence (Hz)','Phase (degré)');