Matlab code to generate the Phase and angle of the sinusoid?

I have been trying to generate the positive ,negative and zero sequence components of the input sinusoid .I have three sinusoidal phases 120 degree separated (Balanced) . The formula to calculate Positive , Negative and zero sequence is given below.

where a = 0.5 + i0.866 or 1<120* (complex number). Now if i generate Va , Vb , Vc inside of a matlab itself and use the formulae I get the correct result. Code below and correct result too .

% Define Parameters
t = 0:0.001:1; % Time vector
A = 1;         % Amplitude
f = 50;        % Frequency in Hz

% Generate Sinusoidal signals with 120 degrees phase shift
v_a = A * exp(1i * (2 * pi * f * t));  % V_a = A * e^(j * 0)
v_b = A * exp(1i * (2 * pi * f * t - 2 * pi / 3)); % V_b = A * e^(j * 120°)
v_c = A * exp(1i * (2 * pi * f * t + 2 * pi / 3)); % V_c = A * e^(j * -120°)

% Call the function to calculate sequence components
[va_pos, va_neg, va_zero] = fcn(v_a, v_b, v_c);


% Plotting the sequence components
figure;

% Plot for Positive Sequence Component
subplot(3,1,1);
plot(t, real(va_pos), 'r', 'DisplayName', 'Real Part'); 
hold on;

title('Positive Sequence Component');
xlim([0 0.3]);
ylim([-1 1]);
xlabel('Time (s)');
ylabel('Amplitude');
legend show;
grid on;

% Plot for Negative Sequence Component
subplot(3,1,2);
plot(t, real(va_neg), 'r', 'DisplayName', 'Real Part');

hold on;
%plot(t, imag(va_neg), 'b', 'DisplayName', 'Imaginary Part');
title('Negative Sequence Component');
xlabel('Time (s)');
ylim([-1 1]);
ylabel('Amplitude');
legend show;
grid on;

% Plot for Zero Sequence Component
subplot(3,1,3);
plot(t, real(va_zero), 'r', 'DisplayName', 'Real Part');
hold on;
%plot(t, imag(va_zero), 'b', 'DisplayName', 'Imaginary Part');
title('Zero Sequence Component');
xlabel('Time (s)');
ylim([-1 1]);
ylabel('Amplitude');
legend show;
grid on;

% Function to calculate sequence components
function [va_pos, va_neg, va_zero] = fcn(v_a, v_b, v_c)
% Define complex exponential for 120 degrees and 240 degrees
alpha = complex(-0.5, 0.866);         % e^(j*120 degrees)
alpha_square = complex(-0.5, -0.866);  % e^(j*240 degrees)

% Calculate positive sequence component
va_pos = (1/3) * (v_a + alpha * v_b + alpha_square * v_c);

% Calculate negative sequence component
va_neg = (1/3) * (v_a + alpha_square * v_b + alpha * v_c);

% Calculate zero sequence component
va_zero = (1/3) * (v_a + v_b + v_c);
end

Now I have Va , Vb , Vc (I generate them now in simulink using sine wave block)and use them as inputs to the simulink function written below (same as one written above).

function [va_pos, va_neg, va_zero] = fcn(v_a, v_b, v_c)

% Define complex exponential for 120 degrees and 240 degrees
alpha = complex(-0.5,0.866);         % e^(j*120 degrees)
alpha_square = complex(-0.5,-0.866);  % e^(j*240 degrees)

% Calculate positive sequence component
va_pos = (1/3) * (v_a + alpha * v_b + alpha_square * v_c);

% Calculate negative sequence component
va_neg = (1/3) * (v_a + alpha_square * v_b + alpha * v_c);

% Calculate zero sequence component
va_zero = (1/3) * (v_a + v_b + v_c);

end

It gives wrong result as show below .(First plot is Positive sequence , second is negative seq. and third is zero sequence . same as the above mentioned plot).

and i suspect it is because it is taking the Va , Vb and Vc in real time (in simulink model).

I tried abs() and phase() to extract the phase and magnitude of the sinusoidal wave separately(in real time function in simulink) but to no avail . I would like to know if there is a way of doing extracting a phase and magnitude of a sinosoidal wave in real like in simulink using MATLAB scripting ?

The error in me not being able to recreate the waveforms like in those subplots was that i was only adding the real part of the complex sinosuidal waveform .

In order for me to extract the waveform I used "PLL block" (to get the frequency 'f') and clock (to get the time values) .

Then the code goes as follows.

function [va_pos,vb_pos,vc_pos,va_neg,vb_neg, vc_neg,   va_zero,vb_zero,vc_zero,va_reconstructed,vb_reconstructed,vc_reconstructed] = fcn(t, f, va_in, vb_in, vc_in)

% Calculate the phase-shifted voltages
v_a = va_in * exp(1i * (2 * pi * f * t));               % V_a = A * e^(j * 0)
v_b = vb_in * exp(1i * (2 * pi * f * t - 2 * pi / 3));  % V_b = A * e^(j * 120°)
v_c = vc_in * exp(1i * (2 * pi * f * t + 2 * pi / 3));  % V_c = A * e^(j * -120°)

% Define complex exponential for 120 degrees and 240 degrees
alpha = exp(1i * 2 * pi / 3);          % e^(j * 120 degrees)
alpha_square = exp(-1i * 2 * pi / 3);  % e^(j * 240 degrees)

% Calculate positive sequence component (full complex value)
va_pos = (1/3) * (v_a + alpha * v_b + alpha_square * v_c);
vb_pos = alpha_square * va_pos;
vc_pos = alpha * va_pos;

% Calculate negative sequence component (full complex value)
va_neg = (1/3) * (v_a + alpha_square * v_b + alpha * v_c);
vb_neg = alpha * va_neg;
vc_neg = alpha_square * va_neg;

% Calculate zero sequence component (real part)
va_zero = (1/3) * (v_a + v_b + v_c);
vb_zero = va_zero;
vc_zero = va_zero;

% Reconstruct the original signal (optional)
% This can be done to check if va_in ≈ va_pos + va_neg + va_zero
va_reconstructed = va_pos + va_neg + va_zero;
vb_reconstructed = vb_pos + vb_neg + vb_zero;
vc_reconstructed = vc_pos + vc_neg + vc_zero;

end

Now after you reconstruct . For plotting the real time values you can plot using real() in time domain.

I was plotting using real and then reconstrucing using the real part as well (That was the error) .

Thank you for the help .