%% Quasi-Resonant Buck 
% State plane solution to Qsw buck with arbitrary J4
%% Circuit definition
Vg = 48;
V = 5;

L = 1e-6;
C = 10e-6; %Not used due to small ripple assumption
Csw = 10e-9;



%% Solve all operating points
% Solve for all combinations of I1 and I4
I1range = linspace(0,50,100);
I4range = flip(linspace(0,50,25));

[I1, I4] = meshgrid(I1range,I4range);

%% State Plane Analysis
% Normalization
Ro = sqrt(L/Csw);
wo = 1/sqrt(L*Csw);
Vbase = Vg;
Ibase = Vg/Ro;

M = V/Vbase;

J1 = I1/Ibase;
J4 = I4/Ibase;

% Interval 1
theta1 = (J1 + J4)./(1-M);

% Interval 2
r2 = sqrt(J1.^2+(1-M)^2);
J2 = sqrt(r2.^2-(M)^2);
beta = pi-atan(J1/(1-M))-atan(J2/M);


% Interval 4
J3 = sqrt(1-2*M+J4.^2); %At ZVS Boundary
delta = pi-atan(J3/M)-atan(J4/(1-M));

% Interval 3
theta3 = (J2+J3)/M;


% Averaging
F = 2*pi./(theta1+theta3+beta+delta);
Jout = F/2/pi.*(theta1/2.*(J1-J4)+theta3/2.*(J2-J3));




%% Denormalize
tbeta = beta/wo;
tdelta = delta/wo;
t1 = theta1/wo;
t3 = theta3/wo;
I1 = J1*Ibase;
I2 = J2*Ibase;
I3 = J3*Ibase;
Iout = Jout*Ibase;
fs = F*wo/2/pi;
Ts = 1./fs;

%% Additional parameters of interest
figure(1)
contourf(I1,-I4,Iout);
axis equal;
cb = colorbar;
xlabel('I_1 [A]');
ylabel('-I_4 [A]');
cb.Label.String = 'I_{out} [A]';

figure(2)
contourf(-I4,Iout,fs*1e-6);
cb = colorbar;
ylabel('I_{out} [A]');
xlabel('-I_4 [A]');
cb.Label.String = 'f_s [MHz]';
xlim([-20, 0])
ylim([0 25])