home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / TIMVEC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.4 KB  |  89 lines
  1. function t=timvec(a,b,c,x0,st,precision)
  2. %TIMVEC Returns a time vector for use with time response functions.
  3. %
  4. %    T=TIMVEC(a,b,c,X0,ST,PRECISION) returns a time vector T
  5. %    for use with state space systems (a,b,c) that have known 
  6. %    initial conditions X0. The time vector can be used,
  7. %    for instance, for step and impulse responses. 
  8. %    
  9. %    TIMVEC attempts to produce a vector that ensures the 
  10. %    output responses have decayed to approx. ST % of their initial 
  11. %    or peak values. 
  12. %
  13. %    More points are used for rapidly changing systems. 
  14. %    Nominal point resolution is determined by setting PRECISION
  15. %    (a good value is 30 points). For more resolution set PRECISION
  16. %    to, for example, 50. 
  17.  
  18. %    A.C.W.Grace 9-21-89
  19. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  20.  
  21.     if nargin<6, precision = 30; end
  22.     if nargin<5, st = 0.005; end
  23.     if nargin<4, x0 = ones(length(a),1); end
  24.     [r,n]=size(c);
  25. % Work out scale factors by looking at d.c. gain, 
  26. % eigenvectors and eigenvalues of the system.
  27.     [m,r]=eig(a);
  28.     r=diag(r);
  29.  
  30. % Equate the  responses to the sum of a set of 
  31. % exponentials(r) multiplied by a vector of magnitude terms(vec).
  32.     [r1,n]=size(c);
  33.     if r1>1, c1=max(abs(c)); else c1=c; end
  34. % Cater for the case when m is singular
  35.     if rcond(m)<eps 
  36.         vec = (c1*m).'.*(1e4*ones(n,1)); 
  37.     else 
  38.         vec=(c1*m).'.*(m\x0);
  39.     end
  40. % Cater for the case when poles and zeros cancel each other out:
  41.     vec=vec+1e-20*(vec==0);
  42. % Cater for the case when eigenvectors are singular:
  43. % d.c. gain (not including d matrix)
  44.     dcgain=c1*x0; % or dcgain =real(sum(vec)) or dcgain=-c(iu,:)/a*b  
  45. % If the d.c gain is small then base the scaling
  46. % factor on the estimated maximum peak in the response.
  47.     if abs(dcgain)<1e-8
  48. % Estimation of the maximum peak in the response.
  49.         estt=(atan2(imag(vec),real(vec))./(imag(r)+eps*(imag(r)==0)))';
  50. % Next line caters for unstable systems .
  51.         estt=(estt<1e12).*estt+(estt>1e12)+(estt==0);
  52.         peak=vec.'*exp(r*abs(estt));
  53.         pk=max([abs(dcgain);abs(peak');eps]);
  54.     else
  55.         pk=dcgain;
  56.     end
  57.     if ~finite(pk), pk = 1/eps; end
  58. % Work out the st% settling time for the responses.
  59. % (or st% peak value settling time for low d.c. gain systems).
  60.     lt=(st*pk./(abs(vec).*(1+(imag(r)~=0))));
  61. % The addition of (r==0) and 0.01*(real(r)==0) below, 
  62. % fixes problem for pure integrators and purely imaginary systems.
  63.     t1=real(log(lt))./(real(r)+(r==0)+0.1*imag(r).*(real(r)==0));
  64.     ts=max(t1);
  65. % Next line is just in case numerical problems are encountered
  66.     if ts<=0, ts=max(abs(t1))+eps; end
  67. % Round the maximum time to an appropriate value for plotting.
  68.     tsn=chop(ts,1,5);
  69.     if abs((ts-tsn)/ts)>0.2, 
  70.         tsn=chop(ts,1,1);
  71.         if abs((ts-tsn)/ts)>0.2
  72.             tsn=chop(ts,2,2);
  73.         end
  74.     end
  75.     ts=tsn;
  76. % Calculate the  number of points to take for the response.
  77. % The first part of this calculation is based on the speed 
  78. % and magnitude of the real eigenvalues.
  79. % The second part is concerned with the oscillatory part of the
  80. % response (imaginary eigenvalues).
  81. % The factors in the imaginary calculation are the damping ratio, the
  82. % residues, the frequency and the final time.         
  83.     r=r+(real(r)==0); % Purely imaginary roots and integrators
  84.     vec = vec + eps*(real(vec)==0); % Can occur with purely imaginary zeros 
  85.     npts=round(max(abs(precision/5*real(vec)./max(abs(real(vec)))*ts.*real(r)) ...
  86.         + precision/2*(log(1+abs(vec./(pk).*( imag(r)./real(r) ).* imag(r)/pi))*ts)));
  87.     if npts>precision*30, npts=precision*30; end
  88.     t=0:ts/npts:ts;      
  89.