%**********************************************************
%
% I. Transient Synchronous Generator Model
%
% Accounts for saturation in the calculations of initial
% conditions.
%
% Includes damper windings.
%
%
% ECE 520  Final Exam
%
% April 12, 2008
%
%**********************************************************
clc;
clear;
% II. Machine Nameplate Data

% A. Nameplate
S_NP = 800*10^6;
V_NP = 24.5*10^3;
f_NP = 60;
Number_Poles = 4;

%**********************************************************
% B. Per Unit System and
%     quantities set by the choice of the per unit system
S_s = S_NP;
V_LLB = V_NP;
w_eB = 2*pi*f_NP;

V_fdB = 270;
I_fdB = 1269;

beta = 2334.9;
X_mdu = 1.0;    % Direct Axis mutual reactance 
R_fd = 1.0;     % Field Resistance



%**********************************************************
% C. Measured Machine Parameters in Per Unit
R_a = 0.006;    % Stator Resistance
H = 3.4;        % Inertia Constant


% Direct-Axis
X_du = 1.13;    % Direct-axis self reactance 
X_dpp = 0.18;   % Direct-axis subtransient reactance
S_e1 = 0.033;
X = 7.57;


% Quadrature-Axis
X_qu = 1.07;    % Quadrature-axis reactance
X_qpp = 0.13;   % Quadrature-axis subtransient reactance


% Time Constants
T_dop = 6.8;    % Direct-axis open-circuit transient time constant
T_dopp = 0.04;  % Direct-axis open-circuit subtransient time constant
T_qopp = 0.06;  % Quadrature-axis open-circuit subtransient time constant


%**********************************************************
% D. Calculated Machine Parameters in Per Unit

% Direct-Axis
L_mdu = X_mdu;
X_L = X_du - X_mdu;


% Quadrature-Axis
X_mqu = X_qu - X_L;


% Field Parameters 
X_ffdu = T_dop* R_fd*w_eB;

X_Lfd = X_ffdu - beta*X_mdu;

%**********************************************************
% E. Leakages equal Reactances in Per Unit 

% Direct-Axis
L_mdu = X_mdu;

L_du = X_du;

L_L = X_L;


% Quadrature-Axis
L_qu = X_qu;

L_mqu = X_mqu;

% Field Parameters 
L_ffdu = X_ffdu;

L_Lfd = X_Lfd;


% Direct-Axis Transient and Subtransient Inductances
L_dp =  L_du - (L_mdu^2)/(L_ffdu/beta);

L_dpp = X_dpp;

L_Lkd = (L_dpp - L_L)*L_mdu*L_Lfd/(L_mdu*L_Lfd - (L_dpp - L_L)*L_ffdu);

L_dopp = L_Lkd + L_mdu*L_Lfd/L_ffdu;

% Quadrature-Axis SubTransient
L_qpp = X_qpp;


%**********************************************************
% III. Specified Initial Conditions

V_at = 1.0*exp(j*0*pi/180);     % "a" phase system terminal voltage
P = 0.0;                        % Real power
Q = -0.4;                       % Reactive power
S = P+j*Q;                      % Complex power

I_a = conj(S/V_at);             % "a" phase current

w_r0 = 1.0;

%**********************************************************
% Calculate Initial Conditions phasor quantities
% Values Independent of Saturation

V_a = V_at + R_a*I_a;            % "a" phase voltage

alpha = angle(V_a);

phi = angle(V_a) - angle(I_a);

E_xa = V_a + j*X_qu*I_a;

gamma = angle(E_xa) - angle(V_at);         % Required I.C.

delta = angle(E_xa) - angle(V_a);

theta_r = angle(E_xa) - pi/2;

I_ad = abs(I_a)*sin(delta+phi)*exp(j*theta_r);

psi_a = -j*V_a;

psi_ad = abs(psi_a)*cos(angle(psi_a)-theta_r)*exp(j*theta_r);

E_aqp = j*(psi_ad + L_dp*I_ad);


%**********************************************************
% Calculated Initial Conditions of dq quantities
% Values Independent of Saturation
V_dq = V_a*exp(-j*theta_r);

mV_dq = abs(V_dq);               % Required I.C.

I_dq = I_a*exp(-j*theta_r);

I_d0 = real(I_dq);               % Required I.C.

I_q0 = imag(I_dq);               % Required I.C.

psi_dq = psi_a*exp(-j*theta_r);

psi_d = real(psi_dq);

psi_q = imag(psi_dq);

%**********************************************************
% Calculate Initial Conditions of quantities and
% parameters dependent on saturation

% E_qp = psi_d + L_dp*I_d

E_qp = E_aqp*exp(-j*theta_r);

E_qp0 = imag(E_qp);             % Required I.C.

S_e = S_e1*abs(E_qp0)^X;

L_md = L_mdu/(1+S_e);

L_d = L_md + L_L;

X_d = L_d;

L_ffd = beta*L_md + L_Lfd;


E_a = E_xa + j*(X_d - X_qu)*I_ad;

E_q = E_a*exp(-j*theta_r);

I_fde = E_qp/(j*L_md);

I_fd = real(I_fde + (L_du - L_dp)*I_d0);

V_fd0 = R_fd*real(I_fd);           % Required I.C.



psi_dqarm = -L_d*I_d0 - j*L_qu*I_q0;

psi_ffd = psi_dq - psi_dqarm;

psi_fd = -beta*L_md*I_d0 + L_ffd*I_fd;


%i_ffdd = E/(j*L_ad)

L_kkd = L_Lkd + L_md;
I_kd0 = 0;                                      % Required I.C.

psi_kd0 = real(L_kkd*I_kd0 + L_md*(I_fd+I_kd0-I_d0));    % Required I.C.

psi_qpp0 = (psi_q + L_qpp*I_q0);                   % Required I.C.


%**********************************************************
% Calculated of Mechanical Torque and Power
Tau_e = psi_d*I_q0 - psi_q*I_d0;

Tau_m = Tau_e;
P_m = Tau_m;                                   % Required I.C.

%P_mchech = P + R_a*(I_d^2 + I_q^2);


initial_time = -5.0;
open_breaker_time = 0.0;

%**********************************************************
% Output the initial conditions to the display
OPF = fopen('InitCond.dat','w');

fprintf(OPF,'Nameplete Quantities \n');
fprintf(OPF,'S_NP = %g  \n', S_NP);
fprintf(OPF,'V_NP = %g  \n', V_NP);
fprintf(OPF,'f_NP = %g  \n', f_NP);
fprintf(OPF,'Number_Poles = %g \n', Number_Poles);

fprintf(OPF,'\n\n Parameters \n');
fprintf(OPF,'R_a = %5.4f  \n', R_a);
fprintf(OPF,'S_e = %5.4f  \n', S_e);
fprintf(OPF,'T_dop = %5.4f secs \n', T_dop);
fprintf(OPF,'L_mdu = %5.4f \n', L_mdu);
fprintf(OPF,'L_md = %5.4f  \n', L_md);
fprintf(OPF,'L_du = %5.4f  \n', L_du);
fprintf(OPF,'L_L = %5.4f  \n', L_L);
fprintf(OPF,'L_qu = %5.4f  \n', L_qu);
fprintf(OPF,'L_mqu = %5.4f  \n', L_mqu);

fprintf(OPF,'\n\n Angles \n');
fprintf(OPF,'phi = %5.3f  \n', phi*180/pi);
fprintf(OPF,'delta = %5.3f degrees \n', delta*180/pi);
fprintf(OPF,'gamma = %5.3f degrees \n', gamma*180/pi);
fprintf(OPF,'theta_r = %5.3f degrees \n', theta_r*180/pi);

fprintf(OPF,'\n\n dq Quantities \n');
fprintf(OPF,'E_xa = %5.3f /_ %5.3f  \n',abs(E_xa),angle(E_xa)*180/pi);
fprintf(OPF,'I_ad = %5.3f /_ %5.3f degrees \n',abs(I_ad),angle(I_ad)*180/pi);
fprintf(OPF,'E_a = %5.3f /_ %5.3f degrees \n',abs(E_a),angle(E_a)*180/pi);
fprintf(OPF,'V_dq = %5.3f /_ %5.3f degrees \n',abs(V_dq),angle(V_dq)*180/pi);
fprintf(OPF,'I_dq = %5.3f /_ %5.3f degrees \n',abs(I_dq),angle(I_dq)*180/pi);
fprintf(OPF,'psi_dq = %5.3f /_ %5.3f degrees \n',abs(psi_dq),angle(psi_dq)*180/pi);
fprintf(OPF,'psi_dqarm = %5.3f /_ %5.3f degrees \n',abs(psi_dqarm),...
             angle(psi_dqarm)*180/pi);
fprintf(OPF,'E_q = %5.3f /_ %5.3f degrees \n',abs(E_q),angle(E_q)*180/pi);
fprintf(OPF,'V_fd0 = %5.3f \n',V_fd0);
fprintf(OPF,'I_fd = %5.3f \n',I_fd);
fprintf(OPF,'psi_ffd = %5.3f \n',psi_ffd);
fprintf(OPF,'psi_fd = %5.3f \n',psi_fd);
fprintf(OPF,'E_qp = %5.3f /_ %5.3f degrees \n',abs(E_qp),angle(E_qp)*180/pi);
fprintf(OPF,'I_d0 = %5.3f \n',I_d0);
fprintf(OPF,'I_q0 = %5.3f \n',I_q0);
fprintf(OPF,'psi_d = %5.5f \n',psi_d);
fprintf(OPF,'psi_q = %5.5f \n',psi_q);
fprintf(OPF,'v_d = %5.5f \n',real(V_dq));
fprintf(OPF,'v_q = %5.5f \n',imag(V_dq));

fprintf(OPF,'\n\n Mechanical Quantities \n');
fprintf(OPF,'Tau_e = %5.3f \n',Tau_e);

fprintf(OPF,'\n\n Initial Conditions (repeated) \n');
fprintf(OPF,'mV_dq = %5.3f \n',abs(mV_dq));
fprintf(OPF,'V_fd0 = %5.3f \n',abs(V_fd0));
fprintf(OPF,'I_d0 = %5.3f \n',I_d0);
fprintf(OPF,'I_q0 = %5.3f \n',I_q0);
fprintf(OPF,'I_kd0 = %5.3f \n',I_kd0);
fprintf(OPF,'E_qp0 = %5.3f \n',E_qp0);
fprintf(OPF,'psi_kd0 = %5.5f \n',psi_kd0);
fprintf(OPF,'psi_qpp0 = %5.5f \n',psi_qpp0);
fprintf(OPF,'gamma = %5.3f degrees \n', gamma*180/pi);
fprintf(OPF,'P_m = %5.3f \n',P_m);
fclose(OPF);

%**********************************************************
% Draw Complex Vector Diagram
compass(1.75+j*0,'r');
hold on;
compass(V_dq,'b');
compass(I_dq,'g');
compass(E_q,'k');
compass(I_fd,'m');
hold off;