% Este programa simula un inversor PWM monofasico con pulso centrado en el
% periodo de conmutación, operando en lazo abierto y alimentando una carga
% no lineal como rectificador a diodos con filtro capacitivo y carga
% resistiva.

clear
close all

fr=50;					% FRECUENCIA DESEJADA DA TENSAO DE SAIDA EM Hz
w=round(2*pi*fr);		% FRECUENCIA ANGULAR DESEJADA DA TENSAO DE SAIDA EM rad/seg
mf=25;          		% NUMERO DE PULSOS PWM POR PERIODO
fs=mf*fr;               % FREQUENCIA DE FREQUENCIA DE COMUTAÇÃO
Ts=1/fs;				% PERIODO DE COMUTAÇÃO)
fa=1*fs;                % FREQUENCIA DE AMOSTRAGEM
Ta=1/fa;                % PERIODO DE AMOSTRAGEM e atualizaçao da malha de precompensaçao da planta.
ms=round(fa/fr);        % NUMERO DE AMOSTRAS POR PERIODO DA TENSAO DE SAIDA
nc=5;				    % NUMERO DE PERIODOS DE SIMULAÇAO
na=nc*ms;               % NUMERO DE ToTAL DE AMOSTRAS POR SIMULAÇAO
p=2048;                 % NUMERO DE PONTOS POR AMOSTRA (PASSO DE SIMULAÇAO) deve ser um numenro par
np=na*p;                % NUMERO TOTAL DE PONTOS DE SIMULAÇAO
Vbase=311;              % VALOR BASE DE TENSAO PARA NORMALIZAÇAO
Ibase=50;               % VALOR BASE DE CORRENTE PARA NORMALIZAÇAO
Vcc=400;

Lf=1500e-6;
Cf=240e-6;
Rc=10;

Cdc=1500e-6;
Rdc=10;
Rs=0.5;

Carga=-1/Rc/Cf;
 
A=[ 0    -1/Lf
   1/Cf  Carga ];
     
B=[1/Lf;0];

F=[0;-1/Cf];

C=[0 1];

[G,H]=c2d(A,B,Ta/p);
[G,W]=c2d(A,F,Ta/p);

%***************************** Variables initialization
kc=1:na*p;                    % Tamanho das variaveis continuas com periodo de amostragem Ta/p
ka=1:na;                      % Tamanho das variaveis discretas com periodo de amostragem Ta

x=zeros(2,length(kc));
upwm=zeros(length(kc),1);
PWM1=zeros(length(kc),1);
PWM2=zeros(length(kc),1);
io=zeros(length(kc),1);
ic=zeros(length(kc),1);
iL=zeros(length(kc),1);
t=zeros(length(kc),1);
portadora=zeros(length(kc),1);
Vdc=zeros(length(kc),1);
Vdc(1:100,1)=280*ones(100,1);
Vo=zeros(length(kc),1);
vref=zeros(length(ka),1);
uca=zeros(length(ka),1);
td=zeros(length(ka),1);

Timer1=0;
flag=1;
ComparReg1=p/2;
ComparReg2=p/2;
sinal_port=1;

km=0;

% Inicia Barra de Progressão
perc=0.01;
passo=0.01;
wb = waitbar(0,'Aguarde, Simulando...');

for k=1:np
    t(k)=k*Ta/p;
    Timer1=Timer1 + sinal_port;

    if flag==1
        flag=0;
        km=km+1;
        td(km)=km*Ta;
        vref(km)=(Vbase/Vbase)*sin(2*pi*fr*td(km));
        uca(km)= vref(km)*p/4;
    end
    
    % Event Manager
    portadora(k)=Timer1;
    
    if Timer1>=p/2
        sinal_port=-1;
        flag=1;        
        % Atualização dos comparadores
        ComparReg1= ( uca(km)+p/4)*(Vbase/Vcc);
        ComparReg2= (-uca(km)+p/4)*(Vbase/Vcc);
    end
    
    if Timer1<0
        sinal_port=1;
    end
    
    if ComparReg1>Timer1
        PWM1(k)=1;
    else
        PWM1(k)=0;
    end
    
    if ComparReg2>Timer1
        PWM2(k)=1;
    else
        PWM2(k)=0;
    end
    
    upwm(k)=(PWM1(k)-PWM2(k))*Vcc;
    %****************************************************
    %                          Plant
    %*****************************************************

    % Retificador não controlado com filtro capacitivo 
    if abs(x(2,k))>=Vdc(k)
        io(k)=sign(x(2,k))*(abs(x(2,k))-Vdc(k))/Rs;
        Vdc(k+1)= Vdc(k) + Ta/p/Cdc*(( abs(x(2,k)) - Vdc(k) )/Rs - Vdc(k)/Rdc );
    else
        io(k)=0;
        Vdc(k+1)= Vdc(k) + Ta/p/Cdc*(-Vdc(k)/Rdc );
    end
    
    Vo(k)=x(2,k);
    io(k)=Vo(k)/Rc;  % Sacar comentario para simular con carga resistiva
    
    iL(k)=x(1,k);
    
    ic(k)=iL(k) - io(k);
    
    x(:,k+1)= G*x(:,k) + H*upwm(k) + W*io(k);
       
    
    % Atualiza a barra de progressão
    if fix((round(np))*perc)==k
        waitbar(perc)
        perc=perc+passo;
    end    
end

close(wb)

figure
plot(t,PWM1*1,t,PWM2*1-1.4,t,upwm/Vcc-3.2)
axis([0 0.02 -4.5 1.5])

figure
plot(t,upwm)
 
figure
plot(t(1:k),x(2,1:k),'r')

% figure
% plot(t(1:k),iL,'r')
% 
% figure
% plot(t,ic)

% figure
% plot(t,io)
 
figure
plot(t(1:k),x(2,1:k),'r',t,io*2,'k'); grid

% figure
% plot(t(1:k),x(2,1:k),'r',t,io*2,'k',t(1:k),Vdc(1:k),'m',td(1:km),vref(1:km)*Vbase);
% grid;

figure
plot(t(1:k),x(2,1:k),'r',t,io*2,'k',td(1:km),vref(1:km)*Vbase)

figure
plot(t(1:k),Vdc(1:k))

THD_tension_salida_OSAP