%**********************************************************
%
% Synchronous Generator
%
%
% ECE 520  Initial Torque Angle for a Salient Pole S/G
%
% April 29, 2008
%
%**********************************************************

clc;
clear;

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

H = 2.5;
w_eB = 2*60*pi;
T_1 = 2.1;
T_2 = 0.3;
T_m = 0.7;

%---------------------------------------------------------------------
% Use Newton-Raphson to solve for delta_0

% Initial guess
delta_0 = pi/4;
ftn = T_m - T_1*sin(delta_0) - T_2*sin(2*delta_0);

while abs(ftn) > 1.0e-014
    delta_0 = delta_0 + ftn/(T_1*cos(delta_0)+2*T_2*cos(2*delta_0));
    ftn = T_m - T_1*sin(delta_0) - T_2*sin(2*delta_0);
end

T_e = T_1*sin(delta_0) + T_2*sin(2*delta_0);


%---------------------------------------------------------------------
% Use Newton-Raphson to solve for delta_max

% Initial guess
delta_max = pi-delta_0;
ftn = T_m - T_1*sin(delta_max) - T_2*sin(2*delta_max);

while abs(ftn) > 1.0e-014
    delta_max = delta_max + ftn/(T_1*cos(delta_max)+2*T_2*cos(2*delta_max));
    ftn = T_m - T_1*sin(delta_max) - T_2*sin(2*delta_max);
end




%---------------------------------------------------------------------
% Use Newton-Raphson to solve for delta_cr

% Initial guess
delta_cr = (delta_0+delta_max)/2;
ftn = T_m*(delta_max - delta_0)...
    + T_1*(cos(delta_max) - cos(delta_cr))...
    + T_2*(cos(2*delta_max) - cos(2*delta_cr))/2.0;

while abs(ftn) > 1.0e-14
    delta_cr = delta_cr - ftn/(T_1*sin(delta_cr)+ T_2*sin(2*delta_cr));
    ftn = T_m*(delta_max - delta_0)...
        + T_1*(cos(delta_max) - cos(delta_cr))...
        + T_2*(cos(2*delta_max) - cos(2*delta_cr))/2.0;
end    

%---------------------------------------------------------------------
% Use Newton-Raphson to solve for delta_0

OPF = fopen('SalientPole.dat','w');

fprintf(OPF,'\n\n Salient Pole Data \n');
fprintf(OPF,'(SalientPole.dat) \n');

fprintf(OPF,'\n\n T_e = %5.4f radians \n', T_e);

fprintf(OPF,'\n\n delta_0 = %5.4f radians \n', delta_0);
fprintf(OPF,' delta_0 = %5.4f degrees \n\n', delta_0*180/pi);

fprintf(OPF,'\n\n delta_max = %5.4f radians \n', delta_max);
fprintf(OPF,' delta_max = %5.4f degrees \n\n', delta_max*180/pi);

fprintf(OPF,'\n\n delta_cr = %5.4f radians \n', delta_cr);
fprintf(OPF,' delta_cr = %5.4f degrees \n\n', delta_cr*180/pi);
fclose(OPF);