%% Método de Tustin com Prewarping
% Professora Gabriela W Gabriel
% 5/5/2024
% Notas de Aula - EES-32
%

clear all; close all; clc;

% Considere a função G(s)=9\frac{s+2}{s+3}, com período de amostragem T=0,2~[{\rm s}].

G=tf(9*[1 2],[1 3])

% A função de transferência discretizada segundo o Método de Tustin, s=\frac{2}{T}\frac{z-1}{z+1}, é dado por
% G_d(z)=9 \frac{10\frac{z-1}{z+1}+2}{10\frac{z-1}{z+1}+3}=9\frac{12}{13}\frac{z-8/12}{z-7/13}=8,3077\frac{z-0,6667}{z-0,5385}

Gd=c2d(G,0.2,'tustin')
w=logspace(-1,2,50);
[M,F]=bode(G,w);
[Md,Fd]=bode(Gd,w);


figure(1); subplot(2,1,1);
semilogx(w,20*log10(M(:)),'k');hold on;
semilogx(w,20*log10(Md(:)),'b*');hold off;
subplot(2,1,2);
semilogx(w,F(:),'k');hold on;
semilogx(w,Fd(:),'b*');hold off;

% Considernado o ajuste de frequência, prewarping em
% \omega_c = \frac{2}{T}\tan{\left\{\omega_d\frac{T}{2}\right\}}
% O ajuste será, uilizando Tustin com prewarping na frequência de \omega_c=2,5~[{\rm rad/s}],
% G_{d,pre}(z)= 9\frac{\frac{2,5}{\tan{0,25}}\frac{z-1}{z+1}+2}{\frac{2,5}{\tan{0,25}}\frac{z-1}{z+1}+3}=9\frac{2,5+2\tan{0,25}}{2,5+3\tan{0,25}}\frac{z-\frac{2,5-2\tan{0,25}}{2,5+2\tan{0,25}}}{z-\frac{2,5-3\tan{0,25}}{2,5+3\tan{0,25}}}=8,2964\frac{z-0,6608}{z-0,5309}

wc=2.5;
Gdpre=c2d(G,0.2,'prewarp',wc)
[Mdpre,Fdpre]=bode(Gdpre,w);
figure(1);
subplot(2,1,1); hold on; semilogx(w,20*log10(Mdpre(:)),'gs'); hold off;
subplot(2,1,2); hold on; semilogx(w,Fdpre(:),'gs'); hold off;

% Porém, na frequência de \omega_c=13~[{\rm rad/s}],
% G_{d,pre}(z)= 9\frac{\frac{13}{\tan{1,3}}\frac{z-1}{z+1}+2}{\frac{13}{\tan{1,3}}\frac{z-1}{z+1}+3}=9\frac{13+2\tan{1,3}}{13+3\tan{1,3}}\frac{z-\frac{13-2\tan{1,3}}{13+2\tan{1,3}}}{z-\frac{13-3\tan{1,3}}{13+3\tan{1,3}}}=7,6382\frac{z-0,2869}{z-0,0921}

wc=13;
Gdpre=c2d(G,0.2,'prewarp',wc)
[Mdpre,Fdpre]=bode(Gdpre,w);
figure(1);
subplot(2,1,1); hold on; semilogx(w,20*log10(Mdpre(:)),'cv'); hold off;
subplot(2,1,2); hold on; semilogx(w,Fdpre(:),'cv'); hold off;

% Por outro lado, aplicando o prewarping nos polos de G(s), teremos
% G'(s)=9\frac{s+10\tan{0,2}}{s+10\tan{0,3}}
% Aplicando novamente Tustin,
% G'_{d}(z)=9 \frac{10\frac{z-1}{z+1}+\tan{0,2}}{10\frac{z-1}{z+1}+\tan{0,3}}=9\frac{1+\tan{0,2}}{1+\tan{0,3}}\frac{z-\frac{1-\tan{0,2}}{1+\tan{0,2}}}{z-\frac{1-\tan{0,3}}{1+\tan{0,3}}}=8,2671\frac{z-0,6629}{z-0,5275}


Glinha=tf(9*[1 10*tan(0.2)],[1 10*tan(0.3)])
Gdlinha=c2d(Glinha,0.2,'tustin')
% [Mdlinha,Fdlinha]=bode(Gdlinha,w);
% 
% figure(1);
% subplot(2,1,1); hold on; semilogx(w,20*log10(Mdlinha(:)),'ro'); hold off;
% subplot(2,1,2); hold on; semilogx(w,Fdlinha(:),'ro'); hold off;

%