function output = ppn6out(t,y,flag)

% PPN6OUT is the output function for PPN6.

% Copywright (c) John C. Polking, Rice University
% Last Modified: January 21, 2001

ppdisp = findobj(get(0,'child'),'flat','name','PPLANE6 Display'); %gca; %
dud = get(ppdisp,'user');
ppdispa = dud.axes;
ud = get(ppdispa,'user');
stop = ud.stop;
gstop = ud.gstop;
col = dud.color.temp;

if (nargin < 3)
if stop
output = 1;
ppn6out(1,y,'done');
else
L = length(t);
if gstop
% Update stop. Are we out of the compute window?
yl = y(:,L);
if any([yl;-yl] - ud.cwind < 0);
stop = 1;
end

% If the derivative function is small we assume there is a sink.
if (ud.i > 10)
yy = [ud.y,y];
dyy = yy(:,1:L) - yy(:,2:(L+1));
DY = ud.DY(:,ones(1,L));
dyy = dyy./DY;
MMf = min(sqrt(sum(dyy.^2)));
if (MMf<=ud.sinkeps*abs(t(1) - t(L)))
stop = 2;
end
else
ud.i = ud.i + 1;
end

% We look for a maximum in the randomly chosen direction. If
% we find one, we compare it with previously found maxima. If
% our new one is close to an old one we stop. If not, we
% record the position.

rr = ud.R*y;

v = [ud.rr,rr];
ud.rr = v(:,[L+1,L+2]);
[m,ii] = max(v(1,:));
%ii = ii(1);
if( 1< min(ii) & max(ii)<L+2 )
kk=0;
turn = ud.turn;
perpeps = ud.perpeps;
paraeps = ud.paraeps;
tk = ud.tk;
while ( (kk<tk) & (~stop) )
kk = kk+1;
if ((abs(v(1,ii)-turn(1,kk))<perpeps) &...
(abs(v(2,ii)-turn(2,kk))<paraeps) )
stop = 3;
end
end
ud.tk = tk + 1;
if tk >= size(turn,2)
ud.turn = [turn,zeros(2,10)];
end
ud.turn(:,tk+1) = v(:,ii);
end
end
output = 0;
ud.stop = stop;
yold = ud.y;
ud.y = y(:,L);
set(ppdispa,'user',ud);

% Finally we plot the newest line segment.

set(ud.line,'Xdata',[yold(1),y(1,:)],'Ydata',[yold(2),y(2,:)]);
drawnow
end

else
switch(flag)
case 'init' % ppn6out(tspan,y0,'init')
ud.y = y(:);
ud.i = 1;

% Set the the line handle.
figure(ppdisp);
ud.line = plot([ud.y(1),ud.y(1)],[ud.y(2),ud.y(2)],...
'color',col,'erase','none');

if gstop
% Chose a random direction & initialize the search for orbit maxima n
% that direction.

theta = rand*2*pi;
ud.R = [cos(theta), sin(theta); -sin(theta),cos(theta)];
qq = ud.R*y;
ud.rr = [qq,qq];
z = ud.DY(1) + sqrt(-1)*ud.DY(2);
w = exp(i*theta);
r = abs(z);
% phi1 = theta + angle(z);
% phi2 = theta - angle(z);
% a1 = r*[cos(phi1),sin(phi1)];
% a2 = r*[cos(phi2),sin(phi2)];
% a1 = abs(a1);
% a2 = abs(a2);
a1 = w*z;
a2 = w*(z');
% a = max(a1(1),a2(1));
% b = max(a1(2),a2(2));
a = max(abs(real([a1,a2])));
b = max(abs(imag([a1,a2])));
ud.perpeps = a/1000; % We want to be real
% close in this direction.
ud.paraeps = b/100; % We do not have to be
% too careful in this direction.
ud.sinkeps = 0.005/dud.settings.refine;
ud.tk = 0;
ud.turn = zeros(2,10);
end
output = 0;
ud.stop = 0;
set(ppdispa,'UserData',ud);

case 'done' % ppn6(t,y,'done');
nstr = get(dud.notice,'string');
if ~isempty(y)
set(ud.line,'Xdata',[ud.y(1),y(1,:)],'Ydata',[ud.y(2),y(2,:)]);
end
switch stop
case 1
nstr{5} = [nstr{5}, ' left the computation window.'];
case 2
ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').'];
nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr];
case 3
nstr{5} = [nstr{5}, ' --> a nearly closed orbit.'];
case 4
nstr{5} = [nstr{5}, ' was stopped by the user.'];
end
set(dud.notice,'string',nstr);
drawnow
output = 1;
end
end