%% Solution

clearvars
format short

vVect = [toRect(1, 0); toRect(0.9629, -2.76); toRect(0.9597, -3.58); toRect(0.9742, -3.96)]
% prove toRect() works fine
% getPolar() is old and working function from
% https://medogan.com/blogs/2022/06/06/
% transmission_line_analysis/transmission_line_sending_edge_vi.html

for v = vVect'
    getPolar(v')
end
testingMatrixEdition = getPolarMatrix(vVect)

%% a)
%

% transformer
a = 0.98;
y_tr = 1/(j*0.08);
y1 = y_tr/a;
y2 = y_tr*((1 - a)/(a^2));
y3 = y_tr*((a - 1)/a);

% shunt susceptances between nodes and gnd
shunt1 = j*0.125 + j*0.1;
shunt2 = j*0.1 + y2;
shunt3 = j*0.125 + j*0.5;
shunt4 = y3;

line12 = getAntiDia(0.02, 0.06);
line13 = getAntiDia(0.02, 0.06);
line23 = getAntiDia(0.05, 0.10);
line24 = -y1;

yBus = [shunt1-(line12+line13), line12, line13, 0;
        line12, shunt2-(line12+line23+line24), line23, line24;
        line13, line23, shunt3-(line13+line23), 0;
        0, line24, 0, shunt4-line24
    ]
yBusPol = getPolarMatrix(yBus)

%% b)

iVect = yBus*vVect
sVect = (vVect).*(conj(iVect))
sVectPol = getPolarMatrix(sVect)

%

function val = toRect(magnitude, angleDrg)
    [real, imag] = pol2cart(deg2rad(angleDrg), magnitude);
    val = complex(real, imag);
end

function strToReturn = getPolar(var)
    precision = 4;
    r = abs(var); theta_deg = rad2deg(angle(var));
    if (abs(r) < 0.1 || abs(r) > 1e3) % display in exponential form
        rstr = num2str(r, "%." + num2str(precision, '%d') + "e"); 
    else
        rstr = num2str(r, "%." + num2str(precision, '%d') + "f");
    end
    if (abs(theta_deg) < 0.1) % display in exponential form
        thetaStr = num2str(theta_deg, "%." + num2str(precision, '%d') + "e");
    else
        thetaStr = num2str(theta_deg, "%." + num2str(precision, '%d') + "f");
    end
    strToReturn = rstr + "∠" + thetaStr + "°";
end

function cellToReturn = getPolarMatrix(matrix)
    mSize = size(matrix);
    cellToReturn = cell(mSize);
    for i = 1:1:mSize(1)
        for j = 1:1:mSize(2)
            cellToReturn(i, j) = cellstr(getPolar(matrix(i, j)));
        end
    end
end

function resToReturn = getAntiDia(r, x)
    z = complex(r, x);
    y = 1/z;
    resToReturn = -y;
end