%**********************************************************
%
% Synchronous Generator
%
%
% ECE 520  Solution to Assignment 4
%
% March 3, 2008
%
%**********************************************************

clc;
clear;

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

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

% Direct-axis
X_d = 1.10;      % Direct-axis self reactance 
L_d = X_d/w_e;   % Direct-axis self inductance 

% Quadrature-axis
X_q = 0.75;      % Quadrature-axis self reactance
L_q = X_q/w_e;   % Quadrature-axis self inductance 

% Stator-field parameters
X_af = 1.00;     % Quadrature-axis self reactance
L_af = X_af/w_e;  % Quadrature-axis self inductance 


%**********************************************************
% II. Specified Initial Conditions

V_at = 1.05*exp(j*0*pi/180);    % "a" phase terminal voltage
S = 0.80*exp(j*30*pi/180);      % Complex power generated
P = real(S);                      % Real power
Q = imag(S);                      % Reactive power

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










%**********************************************************
% III. Calculated Initial Conditions phasor quantities

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

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


psi_a = V_a/(j*w_e);

psi_af = E_a/(j*w_e);

I_f  = psi_af/L_af;

psi_aarm = psi_a - psi_af;

psi_aarm2 = -L_d*I_ad - L_q*I_aq;

psi_af2 = psi_a - psi_aarm2;



%**********************************************************
% IV. Draw Complex Vector Diagram
compass(1.75+j*0,'r');
hold on;
compass(V_a,'b');
compass(I_a,'g');
compass(E_a,'k');
compass(I_f,'m');
hold off;



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

fprintf(OPF,'\n\n ECE 520 Assignment 4 \n');

fprintf(OPF,'\n\n Parameters \n');
fprintf(OPF,'r_s = %5.4f  \n', r_s);
fprintf(OPF,'L_d = %5.4f  \n', L_d);
fprintf(OPF,'L_q = %5.4f  \n', L_q);
fprintf(OPF,'L_sf = %5.4f  \n', L_af);

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_f = %5.3f /_ %5.3f degrees \n\n',abs(I_f), angle(I_f)*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_af = %5.3f /_ %5.3f degrees \n\n',abs(psi_af),...
             angle(psi_af)*180/pi); 
fprintf(OPF,'psi_aarm = %5.3f /_ %5.3f degrees \n\n\n\n',abs(psi_aarm),...
             angle(psi_aarm)*180/pi);
               
fprintf(OPF,'\n\n Calculated Using Currents \n');
fprintf(OPF,'psi_af2 = %5.3f /_ %5.3f degrees \n\n',abs(psi_af2),...
             angle(psi_af2)*180/pi);      
fprintf(OPF,'psi_aarm2 = %5.3f /_ %5.3f degrees \n\n',abs(psi_aarm2),...
             angle(psi_aarm2)*180/pi);
fclose(OPF);