% 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;

iao = 20*10^3*sqrt(2)/sqrt(3)/(sqrt(R^2 + Xl^2));
ibo = 20*10^3*sqrt(2)/sqrt(3)*cos(-2*pi/3)/(sqrt(R^2 + Xl^2));
ico = 20*10^3*sqrt(2)/sqrt(3)*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')

function dx = myOde(t,x)
    
    w = 377;
    S = 659*10^6;
    pf = 0.9;
    Ia = S/ 20*10^3;
    P = 0.9*S;
    R = 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,rs+R,0,0; 
    %     0,0,rs+R,0; 
    %     0,0,0,rf];
    R = diag([rs+R rs+R rs+R 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



