% 6.691 Lecture Example: Simulation of a balanced fault
% We will use a global statement to make variables available to the
% derivative script
global ymd ymq rf rkd rkq omz ra vf H
% parameters of the machine
omz = 2*pi*60;                                % base frequency
H=4;                                          % machine inertia constant
xal = .1;                                     % armature leakage inductance
xad = 2.1;                                    % d-axis mutual inductance
xaq = 1.9;                                    % q-axis mutual inductance
xkdl = .183;                                  % d-axis damper leakage
xkql = .128;                                  % q-axis damper leakage
xfl = .42;                                    % field leakage
rkd = .00706;                                 % d-axis damper resistance
rkq = .079;                                   % q-qxis damper resistance
rf = .00134;                                  % field winding resistance
ra = .0058;                                   % armature winding resistance
xd = xad + xal;                               % assemble full reactances
xq = xaq + xal;
xkd = xad + xkdl;
xkq = xad + xkql;
xf = xad + xfl;
xmd = [xd xad xad; xad xkd xad; xad xad xf];  % d-axis reactance matrix
ymd = inv(xmd);                               % invert this for derivatives
xmq = [xq xaq; xaq xkq];                      % q-axis reactance matrix
ymq = inv(xmq);                               % and its inverse
vf = rf/xad;                                  % required field voltage

% initial state vector is [psid psikd psif psiq psikq omz theta]
% we will generate absolute machine angle as part of the simulation
psid0 = 1;
psikd0 = 1;
psif0 = 1+xfl/xad;
psiq0 = 0;
psikq0 = 0;
theta0 = 0;
state0 = [psid0 psikd0 psif0 psiq0 psikq0 omz theta0];
t0 = 0;                                       % simulation start time
tf = 0.5;                                     % simulation finish time
dt = (tf-t0)/2048;                            % to force a fine-grain output
time = t0:dt:tf;
[t, states] = ode23('sf', time, state0);      % this is the actual simuulation         

psid = states(:, 1);                          % these are the d- axis fluxes
psikd = states(:, 2);
psif = states(:, 3);
psiq = states(:, 4);                          % and the q- axis fluxes
psikq = states(:, 5);
om = states(:, 6);                            % speed
theta = states(:, 7);                         % machine angle

% now we recover the two interesting currents
id = psid .* ymd(1, 1) + psikd .* ymd(1, 2) + psif .* ymd(1, 3);
iq = psiq .* ymq(1, 1) + psikq .* ymq(1, 2);

Te = psid .* iq - psiq .* id;                  % this is generated torque

a = 2*pi/3;                                    % shorthand
ia = id .* cos (theta) - iq .* sin (theta);
ib = id .* cos (theta - a) - iq .* sin (theta - a);
ic = id .* cos (theta + a) - iq .* sin (theta + a);

figure(1)
clf
plot(t, ia)
title('Symmetrical short circuit from open circuit');
ylabel('ia, per unit');
xlabel('Time, s');

figure(2)
clf
plot(t, Te)
     title('Symmetrical short circuit from open circuit')
     ylabel('Per-Unit Torque')
     xlabel('Time, s')