\begin{document}

\section{Appendix: MATLAB Code}
\begin{lstlisting}[
frame=single,
numbers=left,
style=Matlab-Pyglike]
%------------------------------------------------------------------
% Code for SS Model without dampers 
%------------------------------------------------------------------
% define synch gen values
fq = 60;
vf = 377;
nrpm = fq*120/4; % speed of rotor in rev/min
n = nrpm * 0.10472; % speed of rotor in rad/s

S = 659*10^6;
vline = 20*10^3;
vphi = vline/sqrt(3);
pf = 0.9; 
w = 377;

Ia = S/sqrt(3)/vphi;
P = pf*S;
R = P/(Ia^2);
Q = S*sqrt(1-pf^2);
Xl = Q/(Ia^2);
L = Xl / w;

Vm = (20*10^3/sqrt(3))*sqrt(2);
iao = Vm/(sqrt(R^2 + Xl^2));
ibo = Vm*cos(-2*pi/3)/(sqrt(R^2 + Xl^2));
ico = Vm*cos(-4*pi/3)/(sqrt(R^2 + Xl^2));

ifo = 338 / 0.0860;
% ikdo = 0;
% ikqo = 0;

iniCon = [iao;ibo;ico;ifo]
tspan = [0 300];
[t,x] = ode45(@myOde, tspan, iniCon);


plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4));
legend('ia','ib','ic','if')

plot(t,((x(:,1).*(3*Vm*pf))./(n*Xl)))
legend('Torque')

%plot(t,((3.*x(:,1).*Vm)./(n*Xl)))

function dx = myOde(t,x)
    
    w = 377;
    S = 659*10^6;
    pf = 0.9;
    Ia = S/ 20*10^3;
    P = 0.9*S;
    Rl = P/(Ia^2);
    Q = S*sqrt(1-pf^2);
    Xl = Q/(Ia^2);
    L = Xl / w;

    % angle given time t
    
    theta = w*t;               % (2 theta)
    a = -1*(2*pi/3);      % (2 theta - 2pi/3)
    b = -1*(4*pi/3);      % (2 theta - 4pi/3)
    % define values for resistance matrix
    rs = 7.4*10^(-4);
    rf = 0.0860;
    rkd1 = 1.58*10^(-4);
    rkq1 = 1.227*10^(-4);
    
     % R = [rs+Rl, 0,     0, 0,  0,  0; 
     %      0,     rs+Rl, 0, 0,  0,  0; 
     %      0,     0,     rs+Rl, 0,  0,  0; 
     %      0,     0,     0,     rf, 0,  0;
     %      0,     0,     0,     0,  0,  0;
     %      0,     0,     0,     0,  0,  0];
    R = diag([rs+Rl rs+Rl rs+Rl rf]);
    
    % define values for inductance matrix
    pwr = 10^(-3);
    Lsa = 1.95*pwr;
    Lma = 0.80*pwr;
    Lsv = 0.05*pwr;
    Lmv = 0.05*pwr;
    
    Lafm = 26.1*pwr;
    % Lakdm1 = 0.447*pwr;
    % Lakqm1 = 0.670*pwr;
    % 
    % Lfkd1 = 5.05*pwr;
    % Lkd1f = 6.29*pwr;

    Lff = 444*pwr;
    % Lkd1kd1 = 0.1254*pwr;
    % Lkq1kq1 = 0.3762*pwr;
    
    %------------------------------
    Laa = Lsa + Lsv*cos(2*theta);
    Lab = -Lma + Lmv*cos(2*theta + a);
    Lac = -Lma + Lmv*cos(2*theta + b);
    Laf = Lafm*cos(2*theta);
    % Lakd = Lakdm1*cos(2*theta1);
    % Lakq = -Lakqm1*cos(2*theta1);

    % derivative values
    dLaa = -2*Lsv*sin(2*theta);
    dLab = -2*Lmv*sin(2*theta + a);
    dLac = -2*Lmv*sin(2*theta + b);
    dLaf = -2*Lafm*sin(2*theta);
    % dLakd = -2*Lakdm1*sin(2*theta1);
    % dLakq = 2*Lakqm1*sin(2*theta1);

    %-----------------------------------
    Lba = -Lma + Lmv*cos(2*theta + a); 
    Lbb = Lsa + Lsv*cos(2*theta + b);
    Lbc = -Lma + Lmv*cos(2*theta); 
    Lbf = Lafm*cos(2*theta + a);
    % Lbkd = Lakdm1*cos(2*theta2);
    % Lbkq = -Lakqm1*cos(2*theta2);

    % derivative
    dLba = -2*Lmv*sin(2*theta + a); 
    dLbb = -2*Lsv*sin(2*theta + b);
    dLbc = -2*Lmv*sin(2*theta); 
    dLbf = -2*Lafm*sin(2*theta + a);
    % dLbkd = -2*Lakdm1*sin(2*theta2);
    % dLbkq = 2*Lakqm1*sin(2*theta2);

    %-------------------------------
    Lca = -Lma + Lmv*cos(2*theta + b); 
    Lcb = -Lma + Lmv*cos(2*theta); 
    Lcc = Lsa + Lsv*cos(2*theta + a);
    Lcf = Lafm*cos(2*theta + b);
    % Lckd = Lakdm1*cos(2*theta3);
    % Lckq = -Lakqm1*cos(2*theta3);

    % derivative
    dLca = -2*Lmv*sin(2*theta + b); 
    dLcb = -2*Lmv*sin(2*theta); 
    dLcc = -2*Lsv*sin(2*theta + a);
    dLcf = -2*Lafm*sin(2*theta + b);
    % dLckd = -2*Lakdm1*sin(2*theta3);
    % dLckq = 2*Lakqm1*sin(2*theta3);
   
    %---------------------------
    Lfa = Lafm*cos(2*theta);
    Lfb = Lafm*cos(2*theta + a);
    Lfc = Lafm*cos(2*theta + b);

    % derivative
    dLfa = -2*Lafm*sin(2*theta);
    dLfb = -2*Lafm*sin(2*theta + a);
    dLfc = -2*Lafm*sin(2*theta + b);  

    % %-------------------------------
    % Lkda = Lakdm1*cos(2*theta1);
    % Lkdb = Lakdm1*cos(2*theta2);
    % Lkdc = Lakdm1*cos(2*theta3); 
    % 
    % % derivative
    % dLkda = -2*Lakdm1*sin(2*theta1);
    % dLkdb = -2*Lakdm1*sin(2*theta2);
    % dLkdc = -2*Lakdm1*sin(2*theta3);
    % 
    % %-------------------------------
    % Lkqa = -Lakqm1*cos(2*theta1);
    % Lkqb = -Lakqm1*cos(2*theta2);
    % Lkqc = -Lakqm1*cos(2*theta3);
    % 
    % % derivative
    % dLkqa = 2*Lakqm1*sin(2*theta1);
    % dLkqb = 2*Lakqm1*sin(2*theta2);
    % dLkqc = 2*Lakqm1*sin(2*theta3);

    % phew that was a lot
    
    % define inductance matrix and derivative inductance matrix

    L = [Laa+L, Lab, Lac, Laf;
        Lba, Lbb+L, Lbc, Lbf;
        Lca, Lcb, Lcc+L, Lcf;
        Lfa, Lfb, Lfc,   Lff];
        
    dL = [dLaa, dLab, dLac, dLaf;
         dLba, dLbb, dLbc, dLbf;
         dLca, dLcb, dLcc, dLcf;
         dLfa, dLfb, dLfc, 0];

    % V = RI + dL*I + L*dI
    % L*dI = V - RI - dL*I
    % dI = -invL(R + dL)*I + invL*V 
    
    A = -1*(R+dL)\L;
    B = 1\L;
    Vm = (20*10^3/sqrt(3))*sqrt(2);
    u = [Vm*sin(w.*t+pi/2); Vm*sin(w.*t + pi/2 - 2*pi/3); Vm*sin(w.*t + pi/2 - 4*pi/3); 338];
    %K = [Vm*cos(w.*t);Vm*cos(w.*t - (2*pi/3));Vm*cos(w.*t - (4*pi/3));338];
    %K = [0;0;0;377];
    %u = K.*sin(w*t + pi/2);

    dx = A*x + B*u;
end


\end{lstlisting}
\newpage
\begin{lstlisting}[
frame=single,
numbers=left,
style=Matlab-Pyglike]
%------------------------------------------------------------------
% Code for SS Model with dampers 
%------------------------------------------------------------------
% define synch gen values
fq = 60;
vf = 377;
nrpm = fq*120/4; % speed of rotor in rev/min
n = nrpm * 0.10472; % speed of rotor in rad/s
  
S = 659*10^6;
vline = 20*10^3;
vphi = vline/sqrt(3);
pf = 0.9; 
w = 377;

Ia = S/(sqrt(3)*vphi)
P = pf*S;
R = P/(Ia^2)
Q = S*sqrt(1-pf^2);
Xl = Q/(Ia^2)
L = Xl / w;

Vm = (20*10^(3)/sqrt(3))*sqrt(2);
iao = Vm/(sqrt(R^2 + Xl^2));
ibo = Vm*cos(-2*pi/3)/(sqrt(R^2 + Xl^2));
ico = Vm*cos(-4*pi/3)/(sqrt(R^2 + Xl^2));
ifo = 338 / 0.0860;
ikdo = 0;
ikqo = 0;

iniCon = [iao;ibo;ico;ifo;ikdo;ikqo]
tspan = [0 300];
[t,x] = ode45(@myOde, tspan, iniCon);

plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4));
legend('ia','ib','ic','if')

plot(t,x(:,5),t,x(:,6));
legend('idk','idq')

plot(t,((x(:,1).*(3*Vm*pf))./(n*Xl)))
legend('Torque')

function dx = myOde(t,x)
    
    w = 377;
    S = 659*10^6;
    vline = 20*10^3;
    vphi = vline / sqrt(3);
    pf = 0.9;
    Ia = S/(sqrt(3) * vphi);
    P = 0.9*S;
    Rl = P/(Ia^2);
    Q = S*sqrt(1-pf^2);
    Xl = Q/(Ia^2);
    L = Xl / w;

    % angle given time t
    
    theta = w*t;               % (2 theta)
    a = -1*(2*pi/3);      % (2 theta - 2pi/3)
    b = -1*(4*pi/3);      % (2 theta - 4pi/3)
    % define values for resistance matrix
    rs = 7.4*10^(-4);
    rf = 0.0860;
    rkd1 = 1.58*10^(-4);
    rkq1 = 1.227*10^(-4);
    
    % R = [rs+R,0,0,0,0,0; 
    %     0,rs+R,0,0,0,0; 
    %     0,0,rs+R,0,0,0; 
    %     0,0,0,rf,0,0; 
    %     0,0,0,0,rkd1,0; 
    %     0,0,0,0,0,rkq1];
    R = diag([rs+Rl, rs+Rl, rs+Rl, rf, rkd1, rkq1]);

    % define values for inductance matrix
    pwr = 10^(-3);
    Lsa = 1.95*pwr;
    Lma = 0.80*pwr;
    Lsv = 0.05*pwr;
    Lmv = 0.05*pwr;
    
    Lafm = 26.1*pwr;
    Lakdm1 = 0.447*pwr;
    Lakqm1 = 0.670*pwr;

    Lfkd1 = 5.05*pwr;
    Lkd1f = 6.29*pwr;

    Lff = 444*pwr;
    Lkd1kd1 = 0.1254*pwr;
    Lkq1kq1 = 0.3762*pwr;
    
    %------------------------------
    Laa = Lsa + Lsv*cos(2*theta);
    Lab = -Lma + Lmv*cos(2*theta + a);
    Lac = -Lma + Lmv*cos(2*theta + b);
    Laf = Lafm*cos(2*theta);
    Lakd = Lakdm1*cos(2*theta);
    Lakq = -Lakqm1*cos(2*theta);

    % derivative values
    dLaa = -2*Lsv*sin(2*theta);
    dLab = -2*Lmv*sin(2*theta + a);
    dLac = -2*Lmv*sin(2*theta + b);
    dLaf = -2*Lafm*sin(2*theta);
    dLakd = -2*Lakdm1*sin(2*theta);
    dLakq = 2*Lakqm1*sin(2*theta);

    %-----------------------------------
    Lba = -Lma + Lmv*cos(2*theta + a); 
    Lbb = Lsa + Lsv*cos(2*theta + b);
    Lbc = -Lma + Lmv*cos(2*theta); 
    Lbf = Lafm*cos(2*theta + a);
    Lbkd = Lakdm1*cos(2*theta + a);
    Lbkq = -Lakqm1*cos(2*theta + a);

    % derivative
    dLba = -2*Lmv*sin(2*theta + a); 
    dLbb = -2*Lsv*sin(2*theta + b);
    dLbc = -2*Lmv*sin(2*theta); 
    dLbf = -2*Lafm*sin(2*theta + a);
    dLbkd = -2*Lakdm1*sin(2*theta + a);
    dLbkq = 2*Lakqm1*sin(2*theta + a);

    %-------------------------------
    Lca = -Lma + Lmv*cos(2*theta + b); 
    Lcb = -Lma + Lmv*cos(2*theta); 
    Lcc = Lsa + Lsv*cos(2*theta + a);
    Lcf = Lafm*cos(2*theta + b);
    Lckd = Lakdm1*cos(2*theta + b);
    Lckq = -Lakqm1*cos(2*theta + b);

    % derivative
    dLca = -2*Lmv*sin(2*theta + b); 
    dLcb = -2*Lmv*sin(2*theta); 
    dLcc = -2*Lsv*sin(2*theta + a);
    dLcf = -2*Lafm*sin(2*theta + b);
    dLckd = -2*Lakdm1*sin(2*theta + b);
    dLckq = 2*Lakqm1*sin(2*theta + b);
   
    %---------------------------
    Lfa = Lafm*cos(2*theta);
    Lfb = Lafm*cos(2*theta + a);
    Lfc = Lafm*cos(2*theta + b);

    % derivative
    dLfa = -2*Lafm*sin(2*theta);
    dLfb = -2*Lafm*sin(2*theta + a);
    dLfc = -2*Lafm*sin(2*theta + b);  

    %-------------------------------
    Lkda = Lakdm1*cos(2*theta);
    Lkdb = Lakdm1*cos(2*theta + a);
    Lkdc = Lakdm1*cos(2*theta + b); 

    % derivative
    dLkda = -2*Lakdm1*sin(2*theta);
    dLkdb = -2*Lakdm1*sin(2*theta + a);
    dLkdc = -2*Lakdm1*sin(2*theta + b);

    %-------------------------------
    Lkqa = -Lakqm1*cos(2*theta);
    Lkqb = -Lakqm1*cos(2*theta + a);
    Lkqc = -Lakqm1*cos(2*theta + b);

    % derivative
    dLkqa = 2*Lakqm1*sin(2*theta);
    dLkqb = 2*Lakqm1*sin(2*theta + a);
    dLkqc = 2*Lakqm1*sin(2*theta + b);

    % phew that was a lot
    
    % define inductance matrix and derivative inductance matrix

    L = [Laa+L, Lab, Lac, Laf, Lakd, Lakq;
        Lba, Lbb+L, Lbc, Lbf, Lbkd, Lbkq;
        Lca, Lcb, Lcc+L, Lcf, Lckd, Lckq;
        Lfa, Lfb, Lfc,       Lff,   Lfkd1,   0;
        Lkda, Lkdb, Lkdc,    Lkd1f, Lkd1kd1, Lkq1kq1;
        Lkqa, Lkqb, Lkqc,    0,     0,       Lkq1kq1;];
        
    dL = [dLaa, dLab, dLac, dLaf, dLakd, dLakq;
         dLba, dLbb, dLbc, dLbf, dLbkd, dLbkq;
         dLca, dLcb, dLcc, dLcf, dLckd, dLckq;
         dLfa, dLfb, dLfc,    0,0,0;
         dLkda, dLkdb, dLkdc, 0,0,0;
         dLkqa, dLkqb, dLkqc, 0,0,0;];

    % V = RI + dL*I + L*dI
    % L*dI = V - RI - dL*I
    % dI = -invL(R + dL)*I + invL*V 
    
    A = -1*(R+dL)\L;
    B = 1\L;
    Vm = (20*10^(3)/sqrt(3))*sqrt(2);
    u = [Vm*sin(w.*t+pi/2), Vm*sin(w.*t + pi/2 - 2*pi/3), Vm*sin(w.*t + pi/2 - 4*pi/3), 338, 0, 0]';
    %u = [Vm;Vm;Vm;377;0;0];
    %u = K*sin(377.*t + pi/2);
    dx = A*x + B*u;
end

\end{lstlisting}
\end{document}