close all 
clc

%Boost converter parameters
Vin=253;
Vout=395;
D=1-Vin/Vout;
Dp=1-D;
%RLC parameters
Rfull=19.5;
L=110e-6;
C=5e-3;
%ESR parameters
rC=40e-3;
rL=40e-3;

%Bode plot parameters
s=tf('s');
opts=bodeoptions;
opts.FreqUnits='Hz'; 
opts.Grid='on';
freq_range={1e2, 1e7};

Gd0=Vin/Dp^2;
wrC=1/(rC*C);
wRHP=(Dp^2)*(Rfull-rL)/L; 
wres=1/(sqrt(L*C))*sqrt((rL+(Dp^2)*Rfull)/Rfull);
Q=wres/(rL/L+1/(C*(Rfull+rC)));
Gdv_full=Gd0*(1+s/wrC)*(1-s/wRHP)/(1+s/(wres*Q)+(s/wres)^2);

%VMC
fc=2e3;
pm=60;
%Calculating poles and zeros
[gain_fc, phase_fc]=bode(Gdv_full/5*1/79, 2*pi*fc);
pboost=-90+pm-(phase_fc-360);
syms b k
[sol_b, sol_k]=...
    solve(b==tand((pm-phase_fc)/2-rad2deg(pi/4)),...
    sqrt(k)==b+sqrt(b^2+1),...
    b, k);
wz=double(2*pi*fc/sqrt(sol_k));
wp=double(2*pi*fc*sqrt(sol_k));
%Calculating RC values
k_cap=double(2*pi*fc/(gain_fc*sol_k));
syms C1 C2 C3 R2 R3
R1=1e3;
[sol_C1, sol_C2, sol_C3, sol_R2, sol_R3]=...
    solve(C1==wz/(wp*R1*k_cap),...
    C2==1/(R1*k_cap)-C1,...
    C3==1/R1*(1/wz-1/wp),...
    R2==(C1+C2)/(C1*C2*wp),...
    R3==1/(C3*wz)-R1,...
    C1, C2, C3, R2, R3);
%Type III compensator
C=k_cap*(s/wz+1)^2/(s*(s/wp+1)^2);
%Feedback 
H=1/79;
%Loop function
G_loop=C*Gdv_full/5*H;

bodeplot(Gdv_full/5, G_loop, freq_range, opts);
[gain_fc, phase_fc]=bode(G_loop, 2*pi*fc);
margin(G_loop);

disp(['C1 = ',num2str(double(sol_C1))]);
disp(['C2 = ',num2str(double(sol_C2))]);
disp(['C3 = ',num2str(double(sol_C3))]);
disp(['R1 = ',num2str(R1)]);
disp(['R2 = ',num2str(double(sol_R2))]);
disp(['R3 = ',num2str(double(sol_R3))]);