%% Calculating Sending-End Voltage & Current of a Transmission Line

clearvars
format shortEng

%

% Test the custom functions
myNum = 4 + 3*j;
getPolar(myNum)
getPolar(rotateComplexDrg(myNum, 53))
clearvars

%

%%%%%%%%%%% PARAMETERS %%%%%%%%%%%

z = 0.5 + j*0.3; % series impedance per unit length (ohm/km)
y = j*2e-6; % shunt admittance per unit length (mho/km)
L = 500; % Length (km)

S_mag = 150e6;
V_2_ll_mag = 400e3; % line-to-line RMS voltage (V)
pf = 0.91; % power factor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% Calculating the characteristic impedance and the propagation constant
gamma = sqrt(y*z)
Z_c = sqrt(z/y)

gammaPolar = getPolar(gamma)
Z_cPolar = getPolar(Z_c)

%
% b)

% Calculation using the Medium-length model
format shortEng

Z = z*L; % total series line impedance
Y = y*L; % total shunt line admittance
A = 1 + (Z*Y)/2;
B = Z;
C = (1 + (Z*Y)/4)*Y;
D = 1 + (Z*Y)/2;

A_polar = getPolar(A)
B_polar = getPolar(B)
C_polar = getPolar(C)
D_polar = getPolar(D)

phi = acos(pf); %radians
I_2_mag = S_mag/(sqrt(3)*V_2_ll_mag);
% Set v2_line_to_neutral as reference, calc I_2
V_2 = V_2_ll_mag/sqrt(3); % note that it is phase-to-neutral
I_2 = I_2_mag*cos(phi) - j*I_2_mag*sin(phi);

T = [A B; C D];
sendingMatrix = T*[V_2; I_2];

% VsendLL is sqrt(3) times larger in magnitude and 30 degree leading in phase VsenLN
VsendPol_LL = getPolar(rotateComplexDrg(sendingMatrix(1,1)*sqrt(3), 30))
IsendPol = getPolar(sendingMatrix(2,1))

% c)
%

% Calculation using the long line model
A = cosh(gamma*L);
B = Z_c*sinh(gamma*L);
C = (1/Z_c)*sinh(gamma*L);
D = cosh(gamma*L);

A_polar = getPolar(A)
B_polar = getPolar(B)
C_polar = getPolar(C)
D_polar = getPolar(D)

phi = acos(pf); %radians
I_2_mag = S_mag/(sqrt(3)*V_2_ll_mag);
% Set v2_line_to_neutral as reference, calc I_2
V_2 = V_2_ll_mag/sqrt(3); % note that it is phase-to-neutral
I_2 = I_2_mag*cos(phi) - j*I_2_mag*sin(phi);

T = [A B; C D];
sendingMatrix = T*[V_2; I_2];

% VsendLL is sqrt(3) times larger in magnitude and 30 degree leading in phase VsendLN
VsendPol_LL = getPolar(rotateComplexDrg(sendingMatrix(1,1)*sqrt(3), 30))
IsendPol = getPolar(sendingMatrix(2,1))

%

function strToReturn = getPolar(var)
    precision = 3;
    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 numToReturn = rotateComplexDrg(var, theta)
    theta = deg2rad(theta);
    rotMatrix = [cos(theta), -sin(theta); sin(theta), cos(theta)];
    myAns = rotMatrix*[real(var); imag(var)];
    numToReturn = myAns(1,1) + j*myAns(2,1);
end