%**********************************************************
%
% Synchronous Generator
%
%
% ECE 520  Final Project, Part II
%
% April 30, 2008
%
%**********************************************************

clc;
clear;

%**********************************************************
% I. Machine Parameters

% System Frequency
w_B = 2*pi()*60;    % Electical base angular frequency
w_e = 1.0;          % Electrical angular frequency, per unit

% Parameters set by Per Unit system
L_mdu = 1.0;        % Unsaturated direct-axis mutual inductance, per unit
X_mdu = w_e*L_mdu;  % Mutual direct-axis inductance, per unit
R_fd = 1.0;         % Field resistance, per unit
beta = 2334.9;      % Ration of S_sB to S_fB


% Direct-axis
L_du = 1.13;         % Direct-axis self inductance, per unit
L_dpp = 0.18;        % Direct-axis subtransient inductance, per unit
L_L = L_du - L_mdu;  % Direct-axis leakage inductance, per unit
Tau_d0p = 6.80;      % Direct-axis open-circuit transient
                     %     time constant
Tau_d0pp = 0.04;     % Direct-axis open-circuit subtransient
                     %     time constant

S_e1 = 0.033;       % Saturation at E_qp = 1.0
X = 7.57;           % Saturation power

 % Quadrature-axis
L_qu = 0.74;        % Quadrature-axis self inductance, per unit
X_qu = w_e*L_qu;    % Quadrature-axis self reactance, per unit
L_qpp = 0.13;       % Quadrature-axis subtransient inductance, per unit
Tau_q0pp = 0.06;    % Quadrature-axis open ciruit subtransien
                    %     time constant

% Field Parameters
X_ffdu = Tau_d0p*R_fd*w_B;    % Field self reactance, per unit
L_ffdu = X_ffdu/w_e;          % Field self inductance, per unit
L_Lfd = L_ffdu - beta*L_mdu;  % Field winding leakage inductance


% Direct Axis Transient
L_dp = L_du - L_mdu^2/(L_ffdu/beta);  % Direct-axis transient inductance
L_Lkd = (L_dpp-L_L)*L_mdu*L_Lfd/(L_mdu*L_Lfd-(L_dpp-L_L)*L_ffdu);

L_d0pp = L_Lkd+L_mdu*L_Lfd/L_ffdu;

% Turbine and Govenor
LoadReference = 0;
H = 3.4;
D = 1;
T_w = 1;
T_g = 0.2;
T_R = 5;
R_P = 0.05;
R_T = 0.38;

% Fault Parameters
t_fault =10;
t_clear = 0.20;


% Excitor Parameters
V_dq_cmd = 1.0341;

% System Parameters
L_sys = 0.1;
X_sys = w_e*L_sys;
V_sys = 1.0;

%**********************************************************
% II. Specified Initial Conditions
               
V_a = V_dq_cmd;                 % "a" phase terminal voltage
mV_a = abs(V_a);                % Magnitude of the terminal voltage
P = 0.77;                       % Real power



%---------------------------------------------------------------------
% Use Gauss-Seidel to solve for initial conditions
a=0;
Err = 0.1;
while abs(Err) > 1.0e-014
    I_a = (V_a-V_sys)/(j*X_sys);
    S_temp = V_a*conj(I_a);
    Q = imag(S_temp);
    S = complex(P,Q);
    I_a = conj(S/V_a);
    Va_old = V_a;
    V_a = V_sys + j*X_sys*I_a;
    V_a = mV_a*V_a/abs(V_a);
    Err = angle(Va_old)-angle(V_a);
    a = a+1;
    if a>100
        Err = 0;
        'Gauss-Seidel Did NOT solve Network Initialization'
    end
mV_a = abs(V_a);    % Initial stator voltage
end

T_m0 = P/w_e;

%**********************************************************
% III. Calculated Initial Conditions phasor quantities
%      for values independent of saturation

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

alpha = angle(V_a);

E_x = V_a + j*X_qu*I_a;

gamma0 = angle(E_x) - angle(V_sys);

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

theta_r0 = angle(E_x) - pi/2;

theta_EI = angle(E_x) - angle(I_a);

I_ad = abs(I_a)*sin(theta_EI)*exp(j*theta_r0);

I_aq = abs(I_a)*cos(theta_EI)*exp(j*angle(E_x));

psi_a = V_a/(j*w_e);

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

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


%**********************************************************
% IV. Calculated Initial Conditions dq quantities
%      for values independent of saturation

V_dq = V_a*exp(-j*theta_r0);

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

I_d0 = real(I_dq);

I_q0 = imag(I_dq);

psi_dq0 = psi_a*exp(-j*theta_r0);

psi_d0 = real(psi_dq0);

psi_q0 = imag(psi_dq0);


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

E_qp0 = E_aqp*exp(-j*theta_r0);

mE_qp0 = imag(E_qp0);

S_e = S_e1*mE_qp0^X;

L_md = L_mdu/(1+S_e);

L_d = L_md + L_L;

L_kkd = L_Lkd + L_md;

X_d = L_d;

L_ffd = beta*L_md + L_Lfd;

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

E = E_a*exp(-j*theta_r0);

I_fde = E_qp0/(j*L_md);

I_fd0 = I_fde + (L_du - L_dp)*I_d0;

V_fd0 = R_fd*real(I_fd0);

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

psi_fd0 = -beta*L_md*I_d0 + L_ffd*I_fd0;

I_kd0 =0;

psi_qpp0 = psi_q0 + L_qpp*I_q0;

%psi_kd0 = L_kkd*I_kd0 + L_md*(I_fd0+I_kd0-I_d0);
psi_kd0 = mE_qp0 - (L_dp - L_L)*I_d0;

%**********************************************************
% IV. Draw Phasor Diagram
%figure(1);
%compass(1.00+j*0,'r');
%hold on;
%compass(V_a,'b');
%compass(I_a,'g');
%compass(E_a,'m');
%hold off;



%**********************************************************
% V. Draw Complex Vector Diagram
%figure(2);
%compass(V_dq,'b');
%hold on;
%compass(I_dq,'g');
%compass(E_qp,'m');
%compass(I_fd,'k');
%hold off;


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

fprintf(OPF,'\n\n ECE 520 Midterm Examination, Problem 3 \n');

fprintf(OPF,'\n\n Parameters \n');
fprintf(OPF,'L_sys = %5.4f  \n\n', L_sys);

fprintf(OPF,'L_du = %5.4f  \n', L_du);
fprintf(OPF,'L_qu = %5.4f  \n', L_qu);
fprintf(OPF,'L_mdu = %5.4f  \n', L_mdu);
fprintf(OPF,'L_dp = %5.4f  \n', L_dp);
fprintf(OPF,'Tau_d0p = %5.4f  \n', Tau_d0p);
fprintf(OPF,'Tau_d0pp = %5.4f  \n', Tau_d0pp);
fprintf(OPF,'Tau_q0pp = %5.4f  \n', Tau_q0pp);
fprintf(OPF,'S_e1 = %5.4f  \n', S_e1);
fprintf(OPF,'X = %5.4f  \n', X);

fprintf(OPF,'\n\n Input Initial Conditions \n');
fprintf(OPF,'V_sys = %5.3f /_ %5.3f  \n\n',abs(V_sys),angle(V_sys)*180/pi);
fprintf(OPF,'P = %5.4f  \n', P);
fprintf(OPF,'|V_a| = %5.3f  \n\n',abs(V_a));

fprintf(OPF,'\n\n Calculated Initial Conditions \n');
fprintf(OPF,'S = %5.4f /_ %5.4f   \n\n', abs(S), angle(S)*180/pi);
fprintf(OPF,'T_m0 = %5.4f  \n', T_m0);
fprintf(OPF,'Q = %5.4f  \n', Q);
fprintf(OPF,'V_a = %5.3f /_ %5.3f  \n\n',abs(V_a),angle(V_a)*180/pi);

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,'gamma0 = %5.3f degrees \n', gamma0*180/pi);
fprintf(OPF,'theta_r0 = %5.3f degrees \n', theta_r0*180/pi);

fprintf(OPF,'\n\n Phase Quantities \n');
fprintf(OPF,'S_e = %5.4f  \n', S_e);
fprintf(OPF,'I_a = %5.3f /_ %5.3f  \n\n',abs(I_a),angle(I_a)*180/pi);
fprintf(OPF,'V_a = %5.3f /_ %5.3f  \n\n',abs(V_a),angle(V_a)*180/pi);
fprintf(OPF,'E_x = %5.3f /_ %5.3f  \n\n',abs(E_x),angle(E_x)*180/pi);
fprintf(OPF,'I_ad = %5.3f /_ %5.3f degrees \n\n',abs(I_ad),angle(I_ad)*180/pi);
fprintf(OPF,'I_aq = %5.3f /_ %5.3f degrees \n\n',abs(I_aq),angle(I_aq)*180/pi);
fprintf(OPF,'I_fd0 = %5.3f /_ %5.3f degrees \n\n',abs(I_fd0), angle(I_fd0)*180/pi);
fprintf(OPF,'E_a = %5.3f /_ %5.3f degrees \n\n',abs(E_a),angle(E_a)*180/pi);
fprintf(OPF,'psi_a = %5.3f /_ %5.3f degrees \n\n',abs(psi_a),angle(psi_a)*180/pi);
fprintf(OPF,'psi_fd0 = %5.3f /_ %5.3f degrees \n\n',abs(psi_fd0),...
             angle(psi_fd0)*180/pi); 
fprintf(OPF,'psi_dqarm0 = %5.3f /_ %5.3f degrees \n\n\n\n',abs(psi_dqarm0),...
             angle(psi_dqarm0)*180/pi);


fprintf(OPF,'\n\n dq Quantities \n');
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_dq0 = %5.3f /_ %5.3f degrees \n',...
             abs(psi_dq0),angle(psi_dq0)*180/pi);
fprintf(OPF,'psi_dqarm0 = %5.3f /_ %5.3f degrees \n',abs(psi_dqarm0),...
             angle(psi_dqarm0)*180/pi);
fprintf(OPF,'E = %5.3f /_ %5.3f degrees \n',abs(E),angle(E)*180/pi);
fprintf(OPF,'V_fd0 = %5.3f \n',V_fd0);
fprintf(OPF,'I_fd0 = %5.3f \n',I_fd0);
fprintf(OPF,'psi_fd0 = %5.3f \n',psi_fd0);
fprintf(OPF,'psi_kd0 = %5.3f \n',psi_kd0);
fprintf(OPF,'E_qp0 = %5.3f /_ %5.3f degrees \n',...
             abs(E_qp0),angle(E_qp0)*180/pi);
fprintf(OPF,'I_d0 = %5.3f \n',I_d0);
fprintf(OPF,'I_q0 = %5.3f \n',I_q0);
fprintf(OPF,'psi_d0 = %5.5f \n',psi_d0);
fprintf(OPF,'psi_q0 = %5.5f \n',psi_q0);
fprintf(OPF,'v_d = %5.5f \n',real(V_dq));
fprintf(OPF,'v_q = %5.5f \n',imag(V_dq));

fprintf(OPF,'\n\n Initial Conditions (repeated) \n');
fprintf(OPF,'gamma0 = %5.3f degrees \n', gamma0*180/pi);
fprintf(OPF,'V_dq_cmd = %5.3f \n\n',abs(V_dq_cmd));

fprintf(OPF,'T_m0 = %5.3f \n\n',T_m0);

fprintf(OPF,'V_fd0 = %5.3f \n\n',abs(V_fd0));

fprintf(OPF,'mE_qp0 = %5.3f \n',mE_qp0);
fprintf(OPF,'psi_qpp0 = %5.3f \n',psi_qpp0);
fprintf(OPF,'psi_kd0 = %5.3f \n',psi_kd0);

   
fclose(OPF);