%
%

clearvars
freq = 50;
vll_rms = 400;
[vtAArr, vtBArr, vtCArr, tArr] = getThreePhiSinusoid(freq, 40e-3, vll_rms);
threePhiPlot(vtAArr, vtBArr, vtCArr, tArr)
% As seen in the figure, the third harmonic term
% is the same for all the phases
[vtAArr_harm, vtBArr_harm, vtCArr_harm, tArr_harm] = calc3rdHarmonic(freq, 40e-3, vll_rms);
threePhiPlot(vtAArr_harm, vtBArr_harm, vtCArr_harm, tArr_harm)
% Therefore use ph1 harmonic for all phases
harmonicterm = vtAArr_harm;
vtA_harm = vtAArr+harmonicterm; vtB_harm = vtBArr+harmonicterm; vtC_harm = vtCArr+harmonicterm;
threePhiPlot(vtA_harm, vtB_harm, vtC_harm, tArr)

% As seen on figure, the neutral voltage is 0
neutral_wvf_normal = vtAArr + vtBArr + vtCArr;
singlePhiPlot(neutral_wvf_normal, tArr)

% Now the neutral is oscillating
neutral_wvf_harm = vtA_harm + vtB_harm + vtC_harm;
singlePhiPlot(neutral_wvf_harm, tArr)

%The line to line voltages are sinusoidal
threePhiPlot_phph(vtA_harm, vtB_harm, vtC_harm, tArr)

%

function [ph1, ph2, ph3, t] = getThreePhiSinusoid(freq, timeEnd, vllrms)
    omg = 2*pi*freq;
    t = linspace(0, timeEnd, 1e5);
    theta = omg*t;
    vll = vllrms*sqrt(2);
    vln = vll/sqrt(3);
    ph1 = vln*sin(theta);
    ph2 = vln*sin(theta - 2*pi/3);
    ph3 = vln*sin(theta + 2*pi/3);
end

function [ph1 ph2 ph3 t] = calc3rdHarmonic(freq, timeEnd, vllrms)
    omg = 2*pi*freq;
    t = linspace(0, timeEnd, 1e5);
    theta = omg*t;
    vll = vllrms*sqrt(2);
    vln = vll/sqrt(3);
    ph1 = (1/3)*vln*sin(3*theta);
    ph2 = (1/3)*vln*sin(3*(theta - 2*pi/3));
    ph3 = (1/3)*vln*sin(3*(theta + 2*pi/3));
end

function singlePhiPlot(ph1, t)
    figure()
    hold on
    plot(t*1e3, ph1)
    grid on
    grid minor
    xlabel('Time (ms)')
    ylabel('V_{phase-neutral} (V)')
    title("Voltages vs. Time Graph")
    hold off
end

function threePhiPlot(ph1, ph2, ph3, t)
    figure()
    hold on
    plot(t*1e3, ph1)
    plot(t*1e3, ph2)
    plot(t*1e3, ph3)
    legend('v_{t-A}', 'v_{t-B}', 'v_{t-C}')
    grid on
    grid minor
    xlabel('Time (ms)')
    ylabel('V_{phase-neutral} (V)')
    title("Phase to Neutral Line Voltages vs. Time Graph")
    hold off
end

function threePhiPlot_phph(ph1, ph2, ph3, t)
    figure()
    hold on
    plot(t*1e3, ph1-ph2)
    plot(t*1e3, ph2-ph3)
    plot(t*1e3, ph3-ph1)
    legend('v_{t-AB}', 'v_{t-BC}', 'v_{t-CA}')
    grid on
    grid minor
    xlabel('Time (ms)')
    ylabel('V_{phase-neutral} (V)')
    title("Line to Line Voltages vs. Time Graph")
    hold off
end