%**********************************************************
%
% Synchronous Generator
%
%
% ECE 520  Solution to Assignment 5
%
% April 12, 2008
%
%**********************************************************

clc;
clear;

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

w_e = 1.0;          % Electrical angular frequency
r_s = 0.00;         % Stator Resistance

% Direct-axis
L_du = 1.17;        % Direct-axis self inductance
X_du = w_e*L_du;    % Direct-axis self reactance 
L_dp = 0.360;       % Direct-axis transient inductance
T_d0p = 6.80;       % Direct-axis transient open-circuit
                    % 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 
X_qu = w_e*L_qu;    % Quadrature-axis self reactance 


% Stator-field parameters
L_mdu = 1.00;       % Mutual direct-axis inductance 
X_mdu = w_e*L_mdu;  % Mutual direct-axis inductance
beta = 2334.9;      % Ration of S_sB to S_fB

L_Lfd = 230;       % Field winding leakage inductance
R_fd = 1.0;        % Field resistance
%**********************************************************
% II. Specified Initial Conditions

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

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


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

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

V_a = V_at + r_s*I_a; 

alpha = angle(V_a);

E_x = V_a + j*X_qu*I_a;

gamma = angle(E_x) - angle(V_at);

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

mV_dq = abs(V_dq);

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

I_d0 = real(I_dq);

I_q0 = imag(I_dq);

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

psi_d = real(psi_dq);

psi_q = imag(psi_dq);


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

% E_qp = psi_d + L_dp*I_d

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

mE_qp0 = abs(E_qp);

S_e = S_e1*abs(E_qp)^X;

L_md = L_mdu/(1+S_e);

L_L = L_du - L_mdu;

L_d = L_md + L_L;

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_qp/(j*L_md);

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

V_fd = R_fd*abs(I_fd);


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)











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


%**********************************************************
% 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 Assignment 5 Initial Condition File \n');

fprintf(OPF,'\n\n Parameters \n');
fprintf(OPF,'r_s = %5.4f  \n', r_s);
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,'T_d0p = %5.4f  \n', T_d0p);
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,'S = %5.4f /_ %5.4f   \n\n', abs(S), angle(S)*180/pi);
fprintf(OPF,'P = %5.4f  \n', P);
fprintf(OPF,'Q = %5.4f  \n', Q);
fprintf(OPF,'V_at = %5.3f /_ %5.3f  \n\n',abs(V_at),angle(V_at)*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,'gamma = %5.3f degrees \n', gamma*180/pi);
fprintf(OPF,'theta_r0 = %5.3f degrees \n', theta_r0*180/pi);

fprintf(OPF,'\n\n Phase Quantities \n');
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_fd = %5.3f /_ %5.3f degrees \n\n',abs(I_fd), angle(I_fd)*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_fd = %5.3f /_ %5.3f degrees \n\n',abs(psi_fd),...
             angle(psi_fd)*180/pi); 
fprintf(OPF,'psi_dqarm = %5.3f /_ %5.3f degrees \n\n\n\n',abs(psi_dqarm),...
             angle(psi_dqarm)*180/pi);


fprintf(OPF,'\n\n dq Quantities \n');
fprintf(OPF,'E_x = %5.3f /_ %5.3f  \n',abs(E_x),angle(E_x)*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 = %5.3f /_ %5.3f degrees \n',abs(E),angle(E)*180/pi);
fprintf(OPF,'V_fd = %5.3f \n',V_fd);
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 Initial Conditions (repeated) \n');
fprintf(OPF,'mV_dq = %5.3f \n',abs(mV_dq));
fprintf(OPF,'V_fd = %5.3f \n',abs(V_fd));
fprintf(OPF,'I_d0 = %5.3f \n',I_d0);
fprintf(OPF,'I_q0 = %5.3f \n',I_q0);
fprintf(OPF,'I_fd = %5.3f \n',I_fd);
fprintf(OPF,'mE_qp0 = %5.3f \n',mE_qp0);
fprintf(OPF,'gamma = %5.3f degrees \n', gamma*180/pi);
   
fclose(OPF);