%% BAI 2: PHAN TICH MACH DIEN XOAY CHIEU RLC CO HO CAM

% Chuong trinh matlab tinh dong dien trong cac nhanh va cong suat cua nguon
% Su dung 3 phuong phap: Dong nhanh, Dong vong, va The dinh

clear all; clc;

%% DU LIEU DAU VAO

fprintf('=== PHAN TICH MACH DIEN XOAY CHIEU RLC - 3 PHUONG PHAP ===\n\n');

% Cac gia tri dien tro va cam khang

R1 = 200; % Ohm
R2 = 200; % Ohm
R0 = 100; % Ohm
XC1 = 100; % Ohm (cam khang tu dien C1)
L2 = 4; % H (cuon cam L2)

% Tan so goc
omega = 100; % rad/s
f = omega/(2*pi); % Hz

% Cam khang cua L2

XL2 = omega * L2; % = 400 Ohm

% Gia su L0 = 1H de tinh toan

L0 = 1; % H
XL0 = omega * L0; % = 100 Ohm

% He so ho cam M = sqrt(L2*L0)
M = sqrt(L2 * L0); % = sqrt(4*1) = 2 H
XM = omega * M; % = 200 Ohm

fprintf('Cac thong so mach:\n');

fprintf('R1 = R2 = %.0f Ohm\n', R1);
fprintf('R0 = %.0f Ohm\n', R0);
fprintf('XC1 = %.0f Ohm, XL2 = %.0f Ohm, XL0 = %.0f Ohm\n', XC1, XL2, XL0);
fprintf('Ho cam M = %.2f H, XM = %.0f Ohm\n', M, XM);
fprintf('Tan so f = %.2f Hz, omega = %.0f rad/s\n\n', f, omega);

%% NGUON DIEN
% e1 = 220*sqrt(2)*sin(100*t)
E1_max = 220*sqrt(2); % Bien do
E1_rms = 220; % Gia tri hieu dung
phi1 = 0; % Goc pha (do)

% e2 = 200*sin(100*t - 90°)
E2_max = 200; % Bien do
E2_rms = 200/sqrt(2); % Gia tri hieu dung
phi2 = -90; % Goc pha (do)

% Chuyen doi thanh dang so phuc

E1 = E1_rms * exp(1j * phi1 * pi/180); % V
E2 = E2_rms * exp(1j * phi2 * pi/180); % V

fprintf('Nguon dien:\n');
fprintf('E1 = %.2f∠%.0f° V\n', abs(E1), angle(E1)*180/pi);
fprintf('E2 = %.2f∠%.0f° V\n\n', abs(E2), angle(E2)*180/pi);

%% PHUONG PHAP 1: DONG NHANH (NODAL ANALYSIS)

fprintf('=== PHUONG PHAP 1: DONG NHANH (NODAL ANALYSIS) ===\n');
% Dinh nghia cac admittance (dao ham cua impedance)
Y1 = 1/(R1 - 1j*XC1); % Admittance nhanh 1
Y2 = 1/(R2 + R0 + 1j*(XL2 + XL0)); % Admittance nhanh 2 (khong tinh ho cam)
YM = 1/(1j*XM); % Admittance ho cam

% Ma tran admittance nut (gia su 2 nut chinh)

% Nut 1: (Y1 + YM)*V1 - YM*V2 = I1_nguon
% Nut 2: -YM*V1 + (Y2 + YM)*V2 = I2_nguon

% Nguon dong (chuyen doi tu nguon ap)

I1_nguon = E1 * Y1;
I2_nguon = E2 * Y2;

Y_nodal = [Y1 + YM, -YM; -YM, Y2 + YM];

I_nguon = [I1_nguon; I2_nguon];

V_node = Y_nodal \ I_nguon;

V1_nodal = V_node(1);
V2_nodal = V_node(2);

% Tinh dong dien cac nhanh

I1_nodal = (E1 - V1_nodal) * Y1;
I2_nodal = (E2 - V2_nodal) * Y2;
IM_nodal = (V1_nodal - V2_nodal) * YM;

fprintf('Dien ap cac nut:\n');

fprintf('V1 = %.3f∠%.1f° V\n', abs(V1_nodal), angle(V1_nodal)*180/pi);
fprintf('V2 = %.3f∠%.1f° V\n', abs(V2_nodal), angle(V2_nodal)*180/pi);
fprintf('Dong dien cac nhanh:\n');
fprintf('I1 = %.3f∠%.1f° A\n', abs(I1_nodal), angle(I1_nodal)*180/pi);
fprintf('I2 = %.3f∠%.1f° A\n', abs(I2_nodal), angle(I2_nodal)*180/pi);
fprintf('IM = %.3f∠%.1f° A\n\n', abs(IM_nodal), angle(IM_nodal)*180/pi);

%% PHUONG PHAP 2: DONG VONG (MESH ANALYSIS)

fprintf('=== PHUONG PHAP 2: DONG VONG (MESH ANALYSIS) ===\n');

% Ma tran tong tro vong

Z11_mesh = R1 - 1j*XC1; % Tong tro vong 1
Z22_mesh = R2 + R0 + 1j*(XL2 + XL0); % Tong tro vong 2
Z12_mesh = 1j*XM; % Tong tro ho cam
Z21_mesh = 1j*XM;

% Ma tran he so va vector nguon

Z_mesh = [Z11_mesh, Z12_mesh; Z21_mesh, Z22_mesh];
E_mesh = [E1; E2];

% Giai he phuong trinh: Z * I = E

I_mesh = Z_mesh \ E_mesh;
I1_mesh = I_mesh(1);
I2_mesh = I_mesh(2);

fprintf('Dong dien vong:\n');

fprintf('I1 = %.3f∠%.1f° A\n', abs(I1_mesh), angle(I1_mesh)*180/pi);
fprintf('I2 = %.3f∠%.1f° A\n\n', abs(I2_mesh), angle(I2_mesh)*180/pi);

%% PHUONG PHAP 3: THE DINH (NODAL POTENTIAL METHOD)

fprintf('=== PHUONG PHAP 3: THE DINH (NODAL POTENTIAL METHOD) ===\n');
% Chon nut goc (nut 0), xac dinh cac nut khac
% Gia su co 3 nut: nut 1, nut 2, va nut goc (0)
% Nguon E1 noi giua nut 1 va nut 0
% Nguon E2 noi giua nut 2 va nut 0

% Cac admittance
Y1 = 1/(R1 - 1j*XC1); % Y tu nut 1 ve nut 0
Y2 = 1/(R2 + R0 + 1j*(XL2 + XL0)); % Y tu nut 2 ve nut 0
YM = 1/(1j*XM); % Y ho cam giua nut 1 va nut 2

% Ma tran admittance nut (cac nut 1, 2)

% Phuong trinh nut 1: (Y1 + YM)*V1 - YM*V2 = Y1*E1
% Phuong trinh nut 2: -YM*V1 + (Y2 + YM)*V2 = Y2*E2

Y_potential = [Y1 + YM, -YM; -YM, Y2 + YM];

I_inject = [Y1*E1; Y2*E2]; % Dong dien tiem vao cac nut

% Giai he phuong trinh tim the nut

V_potential = Y_potential \ I_inject;
V1_pot = V_potential(1); % The cua nut 1
V2_pot = V_potential(2); % The cua nut 2

fprintf('The cac nut:\n');

fprintf('V1 = %.3f∠%.1f° V\n', abs(V1_pot), angle(V1_pot)*180/pi);
fprintf('V2 = %.3f∠%.1f° V\n', abs(V2_pot), angle(V2_pot)*180/pi);

% Tinh dong dien cac nhanh tu the nut

I1_pot = (E1 - V1_pot) / (R1 - 1j*XC1); % Dong nhanh 1
I2_pot = (E2 - V2_pot) / (R2 + R0 + 1j*(XL2 + XL0)); % Dong nhanh 2
IM_pot = (V1_pot - V2_pot) / (1j*XM); % Dong ho cam

fprintf('Dong dien cac nhanh:\n');

fprintf('I1 = %.3f∠%.1f° A\n', abs(I1_pot), angle(I1_pot)*180/pi);
fprintf('I2 = %.3f∠%.1f° A\n', abs(I2_pot), angle(I2_pot)*180/pi);
fprintf('IM = %.3f∠%.1f° A\n\n', abs(IM_pot), angle(IM_pot)*180/pi);

%% SO SANH KET QUA 3 PHUONG PHAP

fprintf('=== SO SANH KET QUA 3 PHUONG PHAP ===\n');
fprintf('Phuong phap | I1 (A) | I2 (A)
fprintf('-----------------|--------------|-------------
fprintf('Dong nhanh | %.3f∠%.1f° | %.3f∠%.1f°
abs(I1_nodal), angle(I1_nodal)*180/pi, abs(I2_nodal), angle(I2_nodal)*180/pi);
fprintf('Dong vong | %.3f∠%.1f° | %.3f∠%.1f°
abs(I1_mesh), angle(I1_mesh)*180/pi, abs(I2_mesh), angle(I2_mesh)*180/pi);
fprintf('The dinh | %.3f∠%.1f° | %.3f∠%.1f°

abs(I1_pot), angle(I1_pot)*180/pi, abs(I2_pot), angle(I2_pot)*180/pi);

% Kiem tra sai so

error_nodal_mesh_I1 = abs(I1_nodal - I1_mesh);
error_nodal_mesh_I2 = abs(I2_nodal - I2_mesh);
error_mesh_pot_I1 = abs(I1_mesh - I1_pot);
error_mesh_pot_I2 = abs(I2_mesh - I2_pot);

fprintf('Sai so giua cac phuong phap:\n');

fprintf('|I1_nodal - I1_mesh| = %.6f A\n', error_nodal_mesh_I1);
fprintf('|I2_nodal - I2_mesh| = %.6f A\n', error_nodal_mesh_I2);
fprintf('|I1_mesh - I1_potential| = %.6f A\n', error_mesh_pot_I1);
fprintf('|I2_mesh - I2_potential| = %.6f A\n\n', error_mesh_pot_I2);
%% TINH CONG SUAT CUA NGUON VA CAC PHAN TU
fprintf('=== CONG SUAT CUA NGUON VA CAC PHAN TU ===\n');

% Su dung ket qua tu phuong phap dong vong (chinh xac nhat)
I1 = I1_mesh;
I2 = I2_mesh;

% CONG SUAT CUA CAC NGUON

fprintf('CONG SUAT CUA CAC NGUON:\n');

% Cong suat nguon E1 (P = Re(E1 * conj(I1)))

P_nguon1 = real(E1 * conj(I1));
Q_nguon1 = imag(E1 * conj(I1));
S_nguon1 = abs(E1 * I1);

fprintf('Nguon E1:\n');
fprintf(' P1_nguon = %.2f W\n', P_nguon1);
fprintf(' Q1_nguon = %.2f VAR\n', Q_nguon1);
fprintf(' S1_nguon = %.2f VA\n', S_nguon1);

% Cong suat nguon E2

P_nguon2 = real(E2 * conj(I2));
Q_nguon2 = imag(E2 * conj(I2));
S_nguon2 = abs(E2 * I2);

fprintf('Nguon E2:\n');
fprintf(' P2_nguon = %.2f W\n', P_nguon2);
fprintf(' Q2_nguon = %.2f VAR\n', Q_nguon2);
fprintf(' S2_nguon = %.2f VA\n\n', S_nguon2);

% Tong cong suat nguon

P_nguon_total = P_nguon1 + P_nguon2;
Q_nguon_total = Q_nguon1 + Q_nguon2;
S_nguon_total = sqrt(P_nguon_total^2 + Q_nguon_total^2);

fprintf('Tong cong suat cac nguon:\n');

fprintf(' P_nguon_total = %.2f W\n', P_nguon_total);
fprintf(' Q_nguon_total = %.2f VAR\n', Q_nguon_total);
fprintf(' S_nguon_total = %.2f VA\n\n', S_nguon_total);

% CONG SUAT TIEU THU CUA CAC PHAN TU

fprintf('CONG SUAT TIEU THU CUA CAC PHAN TU:\n');

% Cong suat tac dung (P) - chi cac dien tro tieu thu
P_R1 = (abs(I1))^2 * R1; % Cong suat tai R1
P_R2 = (abs(I2))^2 * R2; % Cong suat tai R2
P_R0 = (abs(I2))^2 * R0; % Cong suat tai R0
P_load_total = P_R1 + P_R2 + P_R0; % Tong cong suat tac dung tieu thu

% Cong suat phan khang (Q)

Q_C1 = -(abs(I1))^2 * XC1; % Cong suat phan khang tu C1 (am - dung tinh)
Q_L2 = (abs(I2))^2 * XL2; % Cong suat phan khang tu L2 (duong - cam
tinh)
Q_L0 = (abs(I2))^2 * XL0; % Cong suat phan khang tu L0 (duong - cam
tinh)
Q_load_total = Q_C1 + Q_L2 + Q_L0; % Tong cong suat phan khang tieu thu

% Cong suat bieu kien cua tai

S_load_total = sqrt(P_load_total^2 + Q_load_total^2);

fprintf('Cong suat tac dung tieu thu:\n');

fprintf(' P_R1 = %.2f W\n', P_R1);
fprintf(' P_R2 = %.2f W\n', P_R2);
fprintf(' P_R0 = %.2f W\n', P_R0);
fprintf(' P_load_total = %.2f W\n\n', P_load_total);

fprintf('Cong suat phan khang tieu thu:\n');

fprintf(' Q_C1 = %.2f VAR (dung tinh)\n', Q_C1);
fprintf(' Q_L2 = %.2f VAR (cam tinh)\n', Q_L2);
fprintf(' Q_L0 = %.2f VAR (cam tinh)\n', Q_L0);
fprintf(' Q_load_total = %.2f VAR\n\n', Q_load_total);

fprintf('Cong suat bieu kien tieu thu:\n');

fprintf(' S_load_total = %.2f VA\n\n', S_load_total);

% KIEM TRA CAN BANG CONG SUAT

fprintf('KIEM TRA CAN BANG CONG SUAT:\n');
delta_P = abs(P_nguon_total - P_load_total);
delta_Q = abs(Q_nguon_total - Q_load_total);

fprintf('Can bang cong suat tac dung: %.6f W (sai so)\n', delta_P);
fprintf('Can bang cong suat phan khang: %.6f VAR (sai so)\n', delta_Q);

if delta_P < 1e-10 && delta_Q < 1e-10

fprintf('✓ Can bang cong suat CHINH XAC!\n\n');
else
fprintf('! Co sai so trong can bang cong suat\n\n');
end

% HE SO CONG SUAT CUA TOAN MACH

cos_phi_total = P_load_total / S_load_total;
fprintf('He so cong suat toan mach: cos φ = %.3f\n\n', cos_phi_total);

%% VE BIEN DO THI SO SANH

t = linspace(0, 4*pi/omega, 1000);

% Dong dien tuc thoi tu 3 phuong phap

i1_nodal_t = abs(I1_nodal) * sqrt(2) * sin(omega*t + angle(I1_nodal));
i1_mesh_t = abs(I1_mesh) * sqrt(2) * sin(omega*t + angle(I1_mesh));
i1_super_t = abs(I1_super) * sqrt(2) * sin(omega*t + angle(I1_super));

figure;
subplot(2,1,1);
plot(t*1000, i1_nodal_t, 'b-', 'LineWidth', 2); hold on;
plot(t*1000, i1_mesh_t, 'r--', 'LineWidth', 2);
plot(t*1000, i1_pot_t, 'g:', 'LineWidth', 2);
grid on;
title('So sanh dong dien I1(t) - 3 phuong phap');
xlabel('Thoi gian (ms)');
ylabel('Dong dien (A)');
legend('Dong nhanh', 'Dong vong', 'The dinh', 'Location', 'best');

i2_nodal_t = abs(I2_nodal) * sqrt(2) * sin(omega*t + angle(I2_nodal));

i2_mesh_t = abs(I2_mesh) * sqrt(2) * sin(omega*t + angle(I2_mesh));
i2_pot_t = abs(I2_pot) * sqrt(2) * sin(omega*t + angle(I2_pot));

subplot(2,1,2);
plot(t*1000, i2_nodal_t, 'b-', 'LineWidth', 2); hold on;
plot(t*1000, i2_mesh_t, 'r--', 'LineWidth', 2);
plot(t*1000, i2_pot_t, 'g:', 'LineWidth', 2);
grid on;
title('So sanh dong dien I2(t) - 3 phuong phap');
xlabel('Thoi gian (ms)');
ylabel('Dong dien (A)');
legend('Dong nhanh', 'Dong vong', 'The dinh', 'Location', 'best');

sgtitle('So sanh ket qua 3 phuong phap phan tich mach');

fprintf('=== HOAN THANH PHAN TICH 3 PHUONG PHAP ===\n');

fprintf('Tat ca 3 phuong phap cho ket qua tuong dong nhau.\n');
fprintf('Sai so nho cho thay tinh chinh xac cua cac phuong phap.\n');