home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / nichols.m < prev    next >
Text File  |  1999-04-29  |  4KB  |  117 lines
  1. # Copyright (C) 1998 A. Scottedward Hodel 
  2. #
  3. # This file is part of Octave. 
  4. #
  5. # Octave is free software; you can redistribute it and/or modify it 
  6. # under the terms of the GNU General Public License as published by the 
  7. # Free Software Foundation; either version 2, or (at your option) any 
  8. # later version. 
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT 
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License 
  16. # along with Octave; see the file COPYING.  If not, write to the Free 
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
  18.  
  19. function [mag,phase,w] = nichols(sys,w,outputs,inputs)
  20. # [mag,phase,w] = nichols(sys[,w,outputs,inputs])
  21. # Produce Nichols plot of a system
  22. #
  23. # Compute the frequency response of a system.
  24. # inputs:
  25. #   sys: system data structure (must be either purely continuous or discrete;
  26. #     see is_digit)
  27. #   w: frequency values for evaluation.
  28. #      if sys is continuous, then nichols evaluates G(jw)
  29. #      if sys is discrete, then nichols evaluates G(exp(jwT)), where T=sys.tsam
  30. #         (the system sampling time)
  31. #      default: the default frequency range is selected as follows: (These
  32. #        steps are NOT performed if w is specified)
  33. #          (1) via routine bodquist, isolate all poles and zeros away from
  34. #              w=0 (jw=0 or exp(jwT)=1) and select the frequency
  35. #             range based on the breakpoint locations of the frequencies.
  36. #          (2) if sys is discrete time, the frequency range is limited
  37. #              to jwT in [0,2p*pi]
  38. #          (3) A "smoothing" routine is used to ensure that the plot phase does
  39. #              not change excessively from point to point and that singular
  40. #              points (e.g., crossovers from +/- 180) are accurately shown.
  41. #   outputs, inputs: the indices of the output(s) and input(s) to be used in
  42. #     the frequency response; see sysprune.
  43. # outputs:
  44. #    mag, phase: the magnitude and phase of the frequency response
  45. #       G(jw) or G(exp(jwT)) at the selected frequency values.
  46. #    w: the vector of frequency values used
  47. # If no output arguments are given, nichols plots the results to the screen.
  48. # Descriptive labels are automatically placed.  See xlabel, ylable, title,
  49. # and replot.
  50. #
  51. # Note: if the requested plot is for an MIMO system, mag is set to
  52. # ||G(jw)|| or ||G(exp(jwT))|| and phase information is not computed.
  53.  
  54.   # check number of input arguments given
  55.   if (nargin < 1 | nargin > 4)
  56.     usage("[mag,phase,w] = nichols(sys[,w,outputs,inputs])");
  57.   endif
  58.   if(nargin < 2)
  59.     w = [];
  60.   endif
  61.   if(nargin < 3)
  62.     outputs = [];
  63.   endif
  64.   if(nargin < 4)
  65.     inputs = [];
  66.   endif
  67.  
  68.   [f, w] = bodquist(sys,w,outputs,inputs,"nichols");
  69.  
  70.   [stname,inname,outname] = sysgetsg(sys);
  71.   systsam = sysgetts(sys);
  72.  
  73.   # Get the magnitude and phase of f.
  74.   mag = abs(f);
  75.   phase = arg(f)*180.0/pi;
  76.  
  77.   if (nargout < 1),
  78.     # Plot the information
  79.     if(gnuplot_has_multiplt)
  80.       oneplot();
  81.     endif
  82.     gset autoscale;
  83.     if(gnuplot_has_multiplt)
  84.       gset nokey;
  85.     endif
  86.     clearplot();
  87.     grid("on");
  88.     gset data style lines;
  89.     if(is_digit(sys))
  90.       tistr = "(exp(jwT)) ";
  91.     else
  92.       tistr = "(jw)";
  93.     endif
  94.     xlabel("Phase (deg)");
  95.     if(is_siso(sys))
  96.       title(["Nichols plot of |[Y/U]",tistr,"|, u=", ...
  97.     sysgetsg(sys,"in",1,1), ", y=",sysgetsg(sys,"out",1,1)]);
  98.     else
  99.       title([ "||Y(", tistr, ")/U(", tistr, ")||"]);
  100.       printf("MIMO plot from\n%s\nto\n%s\n",outlist(inname,"    "), ...
  101.         outlist(outname,"    "));
  102.     endif
  103.     if(max(mag) > 0)
  104.       ylabel("Gain in dB");
  105.       md = 20*log10(mag);
  106.     else
  107.       ylabel("Gain |Y/U|")
  108.       md = mag;
  109.     endif
  110.  
  111.     axvec = ax2dlim([vec(phase),vec(md)]);
  112.     axis(axvec);
  113.     plot(phase,md);
  114.     mag = phase = w = [];
  115.   endif
  116. endfunction
  117.