clearvars

% Summon sinusoids representing the grid and SM

freqGrid = 50;
freqSM = 51;
[phA, phB, phC, t] = getThreePhiSinusoid(freqGrid, 4, 380); % Grid
[phR, phS, phT, ~] = getThreePhiSinusoid(freqSM, 4, 380); % SM
diff_RA = phR - phA;

%

% Plot the grid and SM voltages seperately
figure()
hold on
plot(t*1e3, phA)
plot(t*1e3, phB)
plot(t*1e3, phC)
legend('v_{An}', 'v_{Bn}', 'v_{Cn}')
grid on
grid minor
xlabel('Time (ms)')
ylabel('V_{phase-neutral} (V)')
xlim([0, 4*1e3/freqGrid])
title("Grid Voltage vs. Time Graph")
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saveas(gcf, "E.d.grid.emf")
saveas(gcf, "E.d.grid.png")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure()
hold on
plot(t*1e3, phR)
plot(t*1e3, phS)
plot(t*1e3, phT)
legend('v_{Rn}', 'v_{Sn}', 'v_{Tn}')
grid on
grid minor
xlabel('Time (ms)')
ylabel('V_{phase-neutral} (V)')
xlim([0, 4*1e3/freqSM])
title("Machine Terminal Voltage vs. Time Graph")
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saveas(gcf, "E.d.sm.emf")
saveas(gcf, "E.d.sm.png")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% V_A and V_R Plot on X-T mode:
figure()
hold on
plot(t*1e3, phA)
plot(t*1e3, phR)
legend('v_{An}', 'v_{Rn}')
grid on
grid minor
xlabel('Time (ms)')
ylabel('V_{phase-neutral} (V)')
hold off
xlim([0, 4*1e3/freqSM])
title("Phase A and R Voltages vs. Time Graph")
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saveas(gcf, "E.e.emf")
saveas(gcf, "E.e.png")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% Difference plot on X-T mode:
figure()
hold on
plot(t, diff_RA)
legend('v_{R} - v_{A}')
grid on
grid minor
xlabel('Time (s)')
ylabel('v_{R} - v_{A} (V)')
hold off
xlim([0, 2])
title("Voltage Difference Between Phase A and R  vs. Time Graph")
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saveas(gcf, "E.f.emf")
saveas(gcf, "E.f.png")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% V_A and V_R plot on X-Y mode:
figure()
hold on
plot(phR, phA)
peakVal = 380*sqrt(2)/sqrt(3);
plot([-peakVal peakVal], [-peakVal peakVal] ) % Diagonal Line
grid on
grid minor
xlabel('v_{R} (V)')
ylabel('v_{A} (V)')
hold off
% xlim([0, 2])
title("x-y mode")
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saveas(gcf, "E.i,.emf")
saveas(gcf, "E.i.png")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The diagonal line corresponds to the instant when the voltages at A and R are equal. That is, the zero voltage crossings of the V_R -V_A waveform.
% The diagonal line appears periodically with period 1 second.
%

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