home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / NYQUIST.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.2 KB  |  127 lines
  1. function [reout,im,w] = nyquist(a,b,c,d,iu,w)
  2. %NYQUIST Nyquist frequency response for continuous-time linear systems.
  3. %    NYQUIST(A,B,C,D,IU) produces a Nyquist plot from the single input
  4. %    IU to all the outputs of the system:             
  5. %               .                                    -1
  6. %               x = Ax + Bu             G(s) = C(sI-A) B + D  
  7. %               y = Cx + Du      RE(w) = real(G(jw)), IM(w) = imag(G(jw))
  8. %
  9. %    The frequency range and number of points are chosen automatically.
  10. %
  11. %    NYQUIST(NUM,DEN) produces the Nyquist plot for the polynomial 
  12. %    transfer function G(s) = NUM(s)/DEN(s) where NUM and DEN contain
  13. %    the polynomial coefficients in descending powers of s. 
  14. %
  15. %    NYQUIST(A,B,C,D,IU,W) or NYQUIST(NUM,DEN,W) uses the user-supplied
  16. %    freq. vector W which must contain the frequencies, in radians/sec,
  17. %    at which the Nyquist response is to be evaluated.  When invoked 
  18. %    with left hand arguments,
  19. %        [RE,IM,W] = NYQUIST(A,B,C,D,...)
  20. %        [RE,IM,W] = NYQUIST(NUM,DEN,...) 
  21. %    returns the frequency vector W and matrices RE and IM with as many
  22. %       columns as outputs and length(W) rows.  No plot is drawn on the 
  23. %    screen.
  24. %    See also: LOGSPACE,MARGIN,BODE, and NICHOLS.
  25.  
  26. %     J.N. Little 10-11-85
  27. %    Revised ACWG 8-15-89, CMT 7-9-90, ACWG 2-12-91, 6-21-92, 
  28. %               AFP 2-23-93
  29. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  30.  
  31. if nargin==0, eval('exresp(''nyquist'')'), return, end
  32.  
  33. % --- Determine which syntax is being used ---
  34. if (nargin==1),
  35.     error('Wrong number of input arguments.');
  36.  
  37. elseif (nargin==2),    % Transfer function form without frequency vector
  38.     num  = a; den = b; 
  39.     w = freqint2(num,den,30);
  40.     [ny,nn] = size(num); nu = 1;
  41.  
  42. elseif (nargin==3),    % Transfer function form with frequency vector
  43.     num = a; den = b;
  44.     w = c;
  45.     [ny,nn] = size(num); nu = 1;
  46.  
  47. elseif (nargin==4),    % State space system, w/o iu or frequency vector
  48.     error(abcdchk(a,b,c,d));
  49.     w = freqint2(a,b,c,d,30);
  50.     [iu,nargin,re,im]=mulresp('nyquist',a,b,c,d,w,nargout,0);
  51.     if ~iu, if nargout, reout = re; end, return, end
  52.     [ny,nu] = size(d);
  53.  
  54. elseif (nargin==5),    % State space system, with iu but w/o freq. vector
  55.     error(abcdchk(a,b,c,d));
  56.     w = freqint2(a,b,c,d,30);
  57.     [ny,nu] = size(d);
  58.  
  59. else
  60.     error(abcdchk(a,b,c,d));
  61.     [ny,nu] = size(d);
  62.  
  63. end
  64.  
  65. if nu*ny==0, im=[]; w=[]; if nargout~=0, reout=[]; end, return, end
  66.  
  67. % Compute frequency response
  68. if (nargin==2)|(nargin==3)
  69.     gt = freqresp(num,den,sqrt(-1)*w);
  70. else
  71.     gt = freqresp(a,b,c,d,iu,sqrt(-1)*w);
  72. end
  73. ret=real(gt); 
  74. imt=imag(gt);
  75.  
  76. % If no left hand arguments then plot graph.
  77. if nargout==0,
  78.    status = ishold;
  79.    plot(ret,imt,'r-',ret,-imt,'g--')
  80.    set(gca, 'YLimMode', 'auto')
  81.    limits = axis;
  82.    % Set axis hold on because next plot command may rescale
  83.    set(gca, 'YLimMode', 'auto')
  84.    set(gca, 'XLimMode', 'manual')
  85.    hold on
  86.    % Make arrows
  87.    for k=1:size(gt,2),
  88.         g = gt(:,k);
  89.         re = ret(:,k);
  90.         im = imt(:,k);
  91.     sx = limits(2) - limits(1);
  92.     [sy,sample]=max(abs(2*im));
  93.     arrow=[-1;0;-1] + 0.75*sqrt(-1)*[1;0;-1];
  94.     sample=sample+(sample==1);
  95.     reim=diag(g(sample,:));
  96.     d=diag(g(sample+1,:)-g(sample-1,:)); 
  97.     % Rotate arrow taking into account scaling factors sx and sy
  98.     d = real(d)*sy + sqrt(-1)*imag(d)*sx;
  99.     rot=d./abs(d);      % Use this when arrow is not horizontal
  100.     arrow = ones(3,1)*rot'.*arrow;
  101.     scalex = (max(real(arrow)) - min(real(arrow)))*sx/50;
  102.     scaley = (max(imag(arrow)) - min(imag(arrow)))*sy/50;
  103.     arrow = real(arrow)*scalex + sqrt(-1)*imag(arrow)*scaley;
  104.     xy =ones(3,1)*reim' + arrow;
  105.     xy2=ones(3,1)*reim' - arrow;
  106.     [m,n]=size(g); 
  107.     hold on
  108.     plot(real(xy),-imag(xy),'r-',real(xy2),imag(xy2),'g-')
  109.    end
  110.    xlabel('Real Axis'), ylabel('Imag Axis')
  111.  
  112.    limits = axis;
  113.    % Make cross at s = -1 + j0, i.e the -1 point
  114.    if limits(2) >= -1.5  & limits(1) <= -0.5 % Only plot if -1 point is not far out.
  115.     line1 = (limits(2)-limits(1))/50;
  116.     line2 = (limits(4)-limits(3))/50;
  117.     plot([-1+line1, -1-line1], [0,0], 'w-', [-1, -1], [line2, -line2], 'w-')
  118.    end
  119.  
  120.    % Axis
  121.    plot([limits(1:2);0,0]',[0,0;limits(3:4)]','w:');
  122.  
  123.    if ~status, hold off, end    % Return hold to previous status
  124.    return % Suppress output
  125. end
  126. reout = ret; 
  127.