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