home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / FREQINT2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.7 KB  |  72 lines
  1. function [w]=freqint2(a,b,c,d,npts)
  2. %FREQINT2 Auto-ranging algorithm for Nyquist and Nichols frequency
  3. %         response.
  4. %    
  5. %    W=FREQINT2(A,B,C,D,Npts)
  6. %    W=FREQINT2(NUM,DEN,Npts)
  7.  
  8. %    Andy Grace  7-6-90
  9. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  10.  
  11. % Generate more points where graph is changing rapidly.
  12. % Calculate points based on eigenvalues and transmission zeros. 
  13.  
  14. [na,ma] = size(a);
  15.  
  16. if (nargin==3)&(na==1),            % Transfer function form
  17.   npts=c;
  18.   ep=sort([roots(b);roots(a)]);
  19.   n=length(a);
  20.   if n==length(b)
  21.     d=1;
  22.   else
  23.     d=0;
  24.   end
  25. else                % State space form
  26.   if (nargin==3), npts=c; [a,b,c,d] = tf2ss(a,b); end
  27.   z=tzero(a,b,c,d);
  28.   z=z(abs(z)<1.e6);        % Ignore zeros greater than 1.e6
  29.   ep=sort([eig(a);z]);
  30. end
  31.  
  32. [ny,nu] = size(d);  if ny*nu==0, w=[]; return, end
  33.  
  34. if isempty(ep), ep=-1000; end
  35.  
  36. ez=[ep(find(imag(ep)>=0))];
  37. integ = real(ez)<eps; % Any integrators or poles on imaginary axis
  38. ez = ez - 1e-5*(integ) + sqrt(-1)*eps*integ;
  39. [dum,ind]=sort(-abs(real(ez)));
  40. z=[];
  41. npts2=25-6*length(ind);
  42.  
  43. for i=ind'
  44.   npts2=npts2+6;
  45.   npts3=max([7+6*(d(1)~=0),npts2]);
  46.   arez=abs(real(ez(i))); iez=imag(ez(i));
  47.   r1=max([iez-8*arez,1/10*arez]);
  48.   r2=iez+8*arez;
  49.   indr=find(z<=r2&z>=r1);
  50.   npts=npts3*(1+(iez>0))+length(indr);
  51.   if 1.5*iez>arez    
  52.     f1=iez+exp(-1.5:6/npts:4)*arez-exp(-1.5)*arez;
  53.     f2=2*iez-f1; 
  54.     f=[f1,f2(find(f2>r1&f2~=iez))];
  55.   else
  56.     f=logspace(log10(r1),log10(r2),npts);
  57.   end
  58.   z=[z,f];
  59.   z(indr)=[];
  60. end
  61.  
  62. z=sort(z);
  63. mz=z(length(z));
  64. f=logspace(log10(mz),log10(6*mz),10);
  65. z=[linspace(0, 0.95*min(z),10),z,f(2:10),20*mz];
  66. if any(abs(real(ep))<eps) 
  67.     z = z(2:length(z));
  68. end
  69. w=z(:);
  70.  
  71.  
  72.