clc; close all; clear all;
% initial data
ang0 = 16.19; % initial angle
del_t = 0.05; % time interval
end_t = 5.0; % plot from [0, end_t]
clear_t = 0.224; % fault clearing time
K = 2.8125; % this parameter was changed due to frequency change
% Necessary Equation: Pa = 1.6955 - 5.5023 * sin((ang - 0.755) in degree)
% Data Generation:
ang = [16.19]; % list of angles
time = [0]; % list of time
cur_time = 0;
cur_del_ang = 0; % change of angle
cur_ang = ang0;
while (cur_time + del_t) <= end_t
    cur_time = cur_time + del_t;
    Pa = 0;
    if cur_time - del_t < clear_t
        Pa = 1.6955 - 5.5023 * sind(cur_ang - 0.755);
        if(cur_time == del_t)
            Pa = Pa * 0.5;
        end
    else
        Pa = 1.6696 - 6.4934 * sind(cur_ang - 0.847);
    end

    kPa = K * Pa;
    cur_del_ang = cur_del_ang + kPa;
    cur_ang = cur_ang + cur_del_ang;
    time = [time; cur_time];
    ang = [ang; cur_ang];
end
plot(time, ang, 'r.', time, ang, 'k');
title('Swing Curve for machine 2');
xlabel('t, seconds');
ylabel('\delta, electrical degrees');
ylim([0 22]);
grid on;