function varargout = Xtwirl(t,y,flag)
%
%  twirl ODE forcings
%    this file modelled after "help odefile"

switch flag
case ''                      % Return dy/dt = Xf(t,y).
   varargout{1} = Xf(t,y);
case 'events'                % Return [value,isterminal,direction].
  [varargout{1:3}] = events(t,y);
otherwise
   error(['Unknown flag ''' flag '''.']);
end

% -------------------------------------------------------------

function dydt = Xf(t,y)

global r w d

dydt = zeros(size(y));

for m = 0:1:2
  a = r*(w^m)*exp(i*y(1));  za = +1 + a;

  for n = 0:1:2
    b = r*(w^n)*exp(i*y(3));  zb = -1 + b;

    v  = 1./abs(za-zb);

    dydt(2) = dydt(2) + (v.^3) .* real( i*(za-1).*conj(za-zb));
    dydt(4) = dydt(4) + (v.^3) .* real(-i*(zb+1).*conj(za-zb));
  end
end

dydt(1) = y(2);
dydt(3) = y(4);
dydt(2) = dydt(2) - d*y(2);
dydt(4) = dydt(4) - d*y(4);

% -------------------------------------------------------------

function [value,isterminal,direction] = events(t,y)
   value = y(1);   isterminal = 1;   direction = -1;