%% EEET 3016 Power Electronics (2023)
%% Lecture 6: AC-DC single-phase conversion
% Dr. Yee Wei Law \langleyeewei.law@unisa.edu.au\rangle
% Table of Contents
% EEET 3016 Power Electronics (2023)
% Lecture 6: AC-DC single-phase conversion
%     Half-wave uncontrolled rectification
%         Example 1.1
%     Half-wave controlled rectification
%     Full-wave uncontrolled rectification
%         Example 2.1
%         Example 2.2
%     Full-wave controlled rectification
%         Example 3.1
%
%% Half-wave uncontrolled rectification

syms R L V_m omega t i(t) Z tau u
assume([R, L], 'positive');
sol1 = dsolve(diff(i,t)==-R/L*i+V_m/L*sin(omega*t), i(0)==0)
simplify(subs(sol1, R^2+omega^2*L^2, Z^2), 10)
sol2 = expand(V_m/L*int(exp(-R*(t-u)/L)*sin(omega*u), u, 0, t))
simplify(sol1-sol2, 10)

%
%% Example 1.1

clearvars; Vm = 100; omega = 377; R = 100; L = 0.1;
Z = sqrt(R^2+omega^2*L^2)
tau = L/R
theta = atan2(omega*L,R)
% syms beta;
% vpasolve(exp(-beta/(omega*tau))*sin(theta)+sin(beta-theta)==0, beta)
beta = fzero(@(beta) exp(-beta/(omega*tau))*sin(theta)+sin(beta-theta), 1)

% (a) Expression for current in the circuit:

syms vs(t) i(t) t
T = 2*pi/omega;
i = Vm/Z*(exp(-t/tau)*sin(theta)+sin(omega*t-theta));

% (b) Average current:

Iav = vpa(1/T*(int(i, t, 0, beta/omega)))

% (c) RMS current:

Irms = vpa(sqrt(1/T*(int(i^2, t, 0, beta/omega))))

% (d) Output power:

vs = Vm*sin(omega*t);
P = vpa(1/T*(int(vs*i, t, 0, beta/omega)))

% (e) Power factor:

pf = vpa(P/(100/sqrt(2)*Irms))

%
%% Half-wave controlled rectification

syms R L V_m omega i i(t) Z tau u alpha theta
assume([R, L], 'positive');
sol1 = dsolve(diff(i,t)==-R/L*i+V_m/L*sin(omega*t), i(alpha/omega)==0)
expand(sol1)
tau = L/R; theta = atan2(omega*L, R);
sol2 = V_m/Z*(sin(omega*t-theta) - sin(alpha-theta)*exp((alpha-omega*t)/(omega*tau)))
sol2 = expand(sol2)
simplify(subs(sol1-sol2, Z, sqrt(R^2+omega^2*L^2)), 10)

%
%% Full-wave uncontrolled rectification
%% Example 2.1

clearvars; V_m = 100; omega = 2*pi*60; R = 10; L = 10e-3;

% (a) Average load current:

I_0 = 2*V_m/(pi*R)

% (b) Peak-to-peak variation in the load current:

V_2 = 4*V_m/pi*1/(4-1);
I_2 = V_2/sqrt(R^2 + 2^2*omega^2*L^2)
2*I_2

% (c) Average power absorbed by the load and power factor:

I_ms = I_0^2 + 1/2*I_2^2; % Need to take square root of this later to get I_rms
n = 4;
while true
    V_n = 4*V_m/pi * 1/(n^2-1);
    I_n = V_n/sqrt(R^2 + n^2*omega^2*L^2);
    tmp = 1/2*I_n^2
    if  tmp/I_ms > 1e-4 % If tmp is more than 0.01% of I_ms so far
        I_ms = I_ms + tmp;
        n = n + 2;
    else
        break
    end
end
I_rms = sqrt(I_ms)
P_av = I_ms*R
pf = P_av/(V_m/sqrt(2) *I_rms)

% (d) Average and RMS diode currents:

[I_0/2, I_rms/sqrt(2)]

%
%% Example 2.2

clearvars; V_m = 120*sqrt(2); f = 60; omega = 2*pi*f; R = 500; C = 1e-4;
theta = atan2(omega*R*C, -1)
alpha = fzero(@(alpha) sin(alpha)-sin(theta)*exp(-(pi+alpha-theta)/(omega*R*C)), pi/2)
alpha*180/pi
deltaV_o = V_m*(1-sin(alpha))
C = 1/(2*f*R*0.01)

%
%% Full-wave controlled rectification
%% Example 3.1

clearvars; V_rms = 120; V_m = V_rms*sqrt(2); f = 60; R = 10; L = 20e-3; alpha = 60*pi/180;

% (a) DCM? Check if \beta < \pi + \alpha:

omega = 2*pi*f;
Z = sqrt(R^2+omega^2*L^2);
tau = L/R;
theta = atan2(omega*L,R);
syms t i_o(t)
i_o = vpa(V_m/Z*(sin(omega*t-theta) - exp((alpha-omega*t)/(omega*tau))*sin(alpha-theta)), 4)
5.3*exp(1.4)
beta = fzero(@(beta) sin(beta-theta) - exp((alpha-beta)/(omega*tau))*sin(alpha-theta), pi)
beta*180/pi % in degrees

% (b) Expression for load current:

i_o

% where t is between

[alpha/omega beta/omega]

% (c) Average load current, noting that the angular period is now \pi:

T = pi/omega;
I_av = 1/T*int(i_o, t, alpha/omega, beta/omega)

% (d) Power absorbed by the load:

I_rms = sqrt(1/T*int(i_o^2, t, alpha/omega, beta/omega))
P_av = vpa(I_rms^2*R)

% (e) Power factor:

pf = P_av/(V_rms*I_rms)

%
% Plot power factor for \alpha\in[10^\circ, 80^\circ]:

alphas = (10:1:80)*pi/180; pfs = zeros(size(alphas));
tic
for i = 1:numel(alphas)
	alpha = alphas(i);
	beta = fzero(@(beta) sin(beta-theta) - exp((alpha-beta)/(omega*tau))*sin(alpha-theta), pi);
	i_o = vpa(V_m/Z*(sin(omega*t-theta) - exp((alpha-omega*t)/(omega*tau))*sin(alpha-theta)), 4);
	I_rms = sqrt(1/T*int(i_o^2, t, alpha/omega, beta/omega));
	P_av = I_rms^2*R;
	pfs(i) = P_av/(V_rms*I_rms);
	if abs(alpha - 60*pi/180) < 1e-6
		pfs(i) - pf % Sanity check by comparing with pf
	end
end
toc
plot(alphas*180/pi, pfs); xlabel('\alpha (degrees)'); ylabel('Power factor'); grid on;