function dpsi = sf(t, psi);                   % this generates the derivatives

global ymd ymq rf rkd rkq omz ra vf H

psid = psi(1);                                % unpack the state variables
psikd = psi(2);
psif = psi(3);
psiq = psi(4);
psikq = psi(5);
om = psi(6);
theta = psi(7);

Id = ymd * [psid psikd psif]';               % d- axis quantities
Iq = ymq * [psiq psikq]';                    % q- axis quantities
id = Id(1);
ikd = Id(2);
i_f = Id(3);
iq = Iq(1);
ikq = Iq(2);


dpsid = om*psiq - omz*ra*id;                 % and here we compute the derivatives
dpsikd = -omz*rkd*ikd;
dpsif = omz*vf - omz*rf*i_f;
dpsiq = -om*psid - omz*ra*iq;
dpsikq = -omz*rkq*ikq;
dom = (omz/(2*H)) .* (psid*iq - psiq*id);
dtheta = om;

% output is a vector of derivatives: NOTE it is a COLUMN vector!

dpsi = [dpsid dpsikd dpsif dpsiq dpsikq dom dtheta]';