## Copyright (C) 1996,1998 Auburn University.  All Rights Reserved
##
## This file is part of Octave. 
##
## Octave is free software; you can redistribute it and/or modify it 
## under the terms of the GNU General Public License as published by the 
## Free Software Foundation; either version 2, or (at your option) any 
## later version. 
## 
## Octave is distributed in the hope that it will be useful, but WITHOUT 
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
## for more details.
## 
## You should have received a copy of the GNU General Public License 
## along with Octave; see the file COPYING.  If not, write to the Free 
## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 

## -*- texinfo -*-
## @deftypefn {Function File } {[@var{realp}, @var{imagp}, @var{w}] =} nyquist (@var{sys}@{, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol}@})
## @deftypefnx {Function File } {} nyquist (@var{sys}@{, @var{w}, @var{out_idx}, @var{in_idx}, @var{atol}@})
## Produce Nyquist plots of a system; if no output arguments are given, Nyquist
## plot is printed to the screen.
## 
## Compute the frequency response of a system.
## @strong{Inputs} (pass as empty to get default values)
## @table @var
## @item sys
##        system data structure (must be either purely continuous or discrete;
##        see is_digital)
## @item w 
##      frequency values for evaluation.
##      if sys is continuous, then bode evaluates @math{G(jw)}
##      if sys is discrete, then bode evaluates @math{G(exp(jwT))}, where 
##         @math{@var{T}=sysgettsam(@var{sys})} (the system sampling time)
## @item default
##      the default frequency range is selected as follows: (These
##      steps are NOT performed if @var{w} is specified)
## @end table
## @enumerate 
## @item via routine bodquist, isolate all poles and zeros away from 
## @var{w}=0 (@var{jw}=0 or @math{exp(@var{jwT})=1}) and select the frequency
## range based on the breakpoint locations of the frequencies.
## @item if @var{sys} is discrete time, the frequency range is limited
## to @var{jwT} in 
## @ifinfo
## [0,2p*pi]
## @end ifinfo
## @iftex
## $[0,2p*\pi]$
## @end iftex
## @item A "smoothing" routine is used to ensure that the plot phase does
##              not change excessively from point to point and that singular
##              points (e.g., crossovers from +/- 180) are accurately shown.
## @end enumerate
##   outputs, inputs: the indices of the output(s) and input(s) to be used in
##     the frequency response; see sysprune.
## 
## @strong{Inputs} (pass as empty to get default values)
## @table @var
## @item   atol
## for interactive nyquist plots: atol is a change-in-slope tolerance 
## for the of asymptotes (default = 0; 1e-2 is a good choice).  This allows
## the user to ``zoom in'' on portions of the Nyquist plot too small to be
## seen with large asymptotes.
## @end table
## @strong{Outputs}
## @table @var
## @item    realp, imagp
## the real and imaginary parts of the frequency response
##        @math{G(jw)} or @math{G(exp(jwT))} at the selected frequency values.
## @item    w
##  the vector of frequency values used
## @end table
## 
##  If no output arguments are given, nyquist plots the results to the screen.
##  If @var{atol} != 0 and asymptotes are detected then the user is asked 
##     interactively if they wish to zoom in (remove asymptotes)
##  Descriptive labels are automatically placed.
## 
##  Note: if the requested plot is for an MIMO system, a warning message is
##  presented; the returned information is of the magnitude 
##  ||G(jw)|| or ||G(exp(jwT))|| only; phase information is not computed.
## 
## @end deftypefn
 
function [realp, imagp, w] = nyquist (sys, w, outputs, inputs, atol)
   
  ## By R. Bruce Tenison, July 13, 1994
  ## A. S. Hodel July 1995 (adaptive frequency spacing, 
  ##     remove acura parameter, etc.)
  ## Revised by John Ingram July 1996 for system format

  ## Both bode and nyquist share the same introduction, so the common
  ## parts are in a file called bodquist.m.  It contains the part that
  ## finds the number of arguments, determines whether or not the system
  ## is SISO, andd computes the frequency response.  Only the way the
  ## response is plotted is different between the two functions.

  ## check number of input arguments given
  if (nargin < 1 | nargin > 5)
    usage("[realp,imagp,w] = nyquist(sys[,w,outputs,inputs,atol])");
  endif
  if(nargin < 2)
    w = [];
  endif
  if(nargin < 3)
    outputs = [];
  endif
  if(nargin < 4)
    inputs = [];
  endif
  if(nargin < 5)
    atol = 0;
  elseif(!(is_sample(atol) | atol == 0))
    error("atol must be a nonnegative scalar.")
  endif

  ## signal to bodquist who's calling
   
  [f,w] = bodquist(sys,w,outputs,inputs,"nyquist");

  ## Get the real and imaginary part of f.
  realp = real(f);
  imagp = imag(f);

  ## No output arguments, then display plot, otherwise return data.
  if (nargout == 0)
    dnplot = 0;
    while(!dnplot)
      if(gnuplot_has_multiplot)
        oneplot();
        gset key;
      endif
      clearplot();
      grid ("on");
      gset data style lines;
  
      if(is_digital(sys))
        tstr = " G(e^{jw}) ";
      else
        tstr = " G(jw) ";
      endif
      xlabel(["Re(",tstr,")"]);
      ylabel(["Im(",tstr,")"]);
  
      [stn, inn, outn] = sysgetsignals(sys);
      if(is_siso(sys))
        title(sprintf("Nyquist plot from %s to %s, w (rad/s) in [%e, %e]", ...
	  nth(inn,1), nth(outn,1), w(1), w(length(w))) )
      endif
  
      gset nologscale xy;

      axis(axis2dlim([[vec(realp),vec(imagp)];[vec(realp),-vec(imagp)]]));
      plot(realp,imagp,"- ;+w;",realp,-imagp,"-@ ;-w;");

      ## check for interactive plots
      dnplot = 1; # assume done; will change later if atol is satisfied
      if(atol > 0 & length(f) > 2)

        ## check for asymptotes
        fmax = max(abs(f));
        fi = max(find(abs(f) == fmax));
        
        ## compute angles from point to point
        df = diff(f);
        th = atan2(real(df),imag(df))*180/pi;

        ## get angle at fmax
        if(fi == length(f)) fi = fi-1; endif
        thm = th(fi);
    
        ## now locate consecutive angles within atol of thm
        ith_same = find(abs(th - thm) < atol);
        ichk = union(fi,find(diff(ith_same) == 1));

        ## locate max, min consecutive indices in ichk
        loval = max(complement(ichk,1:fi));
        if(isempty(loval)) loval = fi;
        else               loval = loval + 1;   endif

        hival = min(complement(ichk,fi:length(th)));
        if(isempty(hival))  hival = fi+1;      endif

        keep_idx = complement(loval:hival,1:length(w));

        if(length(keep_idx))
          resp = input("Remove asymptotes and zoom in (y or n): ",1);
          if(resp(1) == "y")
            dnplot = 0;                 # plot again
            w = w(keep_idx);
            f = f(keep_idx);
            realp = real(f);
            imagp = imag(f);
          endif
        endif

     endif
    endwhile
    w = [];
    realp=[];
    imagp=[];
  endif

endfunction
