function [Ybus, Yf, Yt] = creat_Y(baseMVA, bus, branch)

% This function creates the admittance matrices Ybus, Yf and Yt for a given power system network
% Inputs:
baseMVA: Base MVA of the power system
% bus: Bus data matrix containing bus number, type, P, Q, voltage magnitude and angle
% branch: Branch data matrix containing from bus, to bus, resistance, reactance, and line charging data

% Outputs:
% Ybus: Admittance matrix of the entire power system
% Yf: Admittance matrix of the "from" end of all branches
% Yt: Admittance matrix of the "to" end of all branches

% Number of buses
nb = size(bus, 1);

% Number of branches
nl = size(branch, 1);

% Initialize the admittance matrices
Ybus = zeros(nb, nb);
Yf = zeros(nl, nb);
Yt = zeros(nl, nb);

% Compute the series admittance of each branch
y_series = (branch(:, 3) + 1i*branch(:, 4)).^-1;

% Compute the shunt admittance of each branch
y_shunt = 1i*branch(:, 5)/2;

% Compute the admittance matrix Yf for the "from" end of each branch
for i = 1:nl
Yf(i, branch(i, 1)) = -y_series(i) - y_shunt(i);
Yf(i, branch(i, 2)) = y_series(i);
end

% Compute the admittance matrix Yt for the "to" end of each branch
for i = 1:nl
Yt(i, branch(i, 1)) = y_series(i);
Yt(i, branch(i, 2)) = -y_series(i) - y_shunt(i);
end

% Compute the admittance matrix Ybus for the entire power system
Ybus = Yf.'*Yf + Yt.'*Yt;
for i = 1:nb
Ybus(i, i) = Ybus(i, i) + sum(Yf(:, i)) + sum(Yt(:, i));
end

% Scale the admittance matrices by the base MVA
Ybus = Ybus/baseMVA;
Yf = Yf/baseMVA;
Yt = Yt/baseMVA;

end