% demo 8_1
% using degrees as independent variable.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Square-wave operation.
% Purely inductive load with back-emf.

% To save figure as vector image,
% go to "Print Preview" and set correct orientation (eg. landscape,
% fill page),
% then "Save As..." and pick pdf, svg etc...

f = 52; % Hz
omega = 2*pi*f; % rad/s
periods = 2;
t = linspace(0, periods/f, 1000); % generate time var for two periods at frequency f.
theta = r2d(omega*t); % calculate phase angle at time t.
fi = 120; % phase angle between phases.
L = 100e-3;

Vd = 256.5; %V

VAN = Vd*(square(d2r(theta)) + 1)/2;
VBN = Vd*(square(d2r(theta) + d2r(fi)) + 1)/2;
VCN = Vd*(square(d2r(theta) + 2*d2r(fi)) + 1)/2;
VAn = 2/3 * VAN - 1/3 * (VBN + VCN);

VAn1 = 2/pi * Vd * sin(d2r(theta));


m = 6; n = 1;
subplot(m, n, 1)
plot(theta, VAN) %plot with x axis normalized to n*pi
xticks(0:60:720);
xlabel('\theta');, ylabel('VAN'); grid on;

subplot(m, n, 2)
plot(theta, VBN) %plot with x axis normalized to n*pi
xticks(0:60:720);
xlabel('\theta'); ylabel('VBN'); grid on;

subplot(m, n, 3)
plot(theta, VCN) %plot with x axis normalized to n*pi
xticks(0:60:720);
xlabel('\theta'); ylabel('VCN'); grid on;

subplot(m, n, 4)
plot(theta, VAn, theta, VAn1, '--')
xticks(0:60:720);
ylabel('VAn'); xlabel('\theta'); grid on;

subplot(m, n, 5)
plot(theta, VAn - VAn1)
xticks(0:60:720);
ylabel('VA_{n_{ripple}}'); xlabel('\theta'); grid on; 

subplot(m, n, 6)
Z = (1/(r2d(omega).*L)) .* cumtrapz(theta,(VAn - VAn1));
Z = Z - Z(1000 / (4*periods)); % 1000/8 = 4pi/8 = pi/2 = 90deg
plot(theta, Z)
xticks(0:60:720);
ylabel('i_{A_{ripple}}'); grid on;
xlabel('\theta');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Trying to find peak current ripple using 
%i_peak = i(t_0) + 1/omega*L * integral(V_An_ripple, pi/2, pi)

fun1 = @(x) (2/3 - (2/pi)*sin(x)); 
q1 = integral(fun, pi/2, (2/3)*pi)

fun2 = @(x) (1/3 - (2/pi)*sin(x));
q2 = integral(fun, (2/3)*pi, pi)
ans = (q1 + q2) * Vd/(omega.*L)