## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
## 02110-1301 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 @code{is_digital})
## @item w
## frequency values for evaluation.
## If sys is continuous, then bode evaluates @math{G(@var{jw})}; 
## if sys is discrete, then bode evaluates @math{G(exp(@var{jwT}))},
## where @var{T} is the system sampling time.
## @item default
## the default frequency range is selected as follows: (These
## steps are @strong{not} performed if @var{w} is specified)
## @enumerate
## @item via routine @command{__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
## @tex
## $ [ 0,2  p \pi ] $
## @end tex
## @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
## @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
## @itemx   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 @acronym{MIMO} system, a warning message is
## presented; the returned information is of the magnitude
## @iftex
## @tex
## $ \Vert G(jw) \Vert $ or $ \Vert G( {\rm exp}(jwT) \Vert $
## @end tex
## @end iftex
## @ifinfo
## ||G(jw)|| or ||G(exp(jwT))||
## @end ifinfo
## only; phase information is not computed.
## @end deftypefn

## Author: R. Bruce Tenison <btenison@eng.auburn.edu>
## Created: July 13, 1994
## A. S. Hodel July 1995 (adaptive frequency spacing,
##     remove acura parameter, etc.)
## Revised by John Ingram July 1996 for system format

function [realp, imagp, w] = nyquist (sys, w, outputs, inputs, atol)

  ## 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, sys] = __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)
      oneplot();
      __gnuplot_set__ key;
      clearplot();
      grid ("on");
      __gnuplot_set__ 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]", ...
          inn{1}, outn{1}, w(1), w(length(w))) )
      endif

      __gnuplot_set__ 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