home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / lsim < prev    next >
Encoding:
Text File  |  1995-11-15  |  2.9 KB  |  122 lines
  1. //---------------------------------------------------------------------------------
  2. //
  3. // lsim
  4. //
  5. // Syntax: G=lsim(A,B,C,D,U,T,X0)
  6. //
  7. // This routine plots the time response of the linear system descried by,
  8. //          .
  9. //          x = Ax + Bu
  10. //          y = Cx + Du
  11. //
  12. // to the input time history U. The routine may be called in several fashions:
  13. //
  14. // (1)  G=lsim(A,B,C,D,U,T)
  15. //      (This produces the time history with Zero initial conditions for
  16. //       a state-space model).
  17. //
  18. // (2)  G=lsim(A,B,C,D,U,T,X0)
  19. //      (This produces the time history with initial conditions for
  20. //       a state-space model).
  21. //
  22. // (3)   G=lsim(NUM,DEN,U,T)
  23. //       (This produces the time history with Zero initia conditions for
  24. //        a transfer function model, G(s)=NUM(s)/DEN(s) ).
  25. //
  26. // For all 3 cases, the matrix U must have as many columns as there
  27. // are inputs, U. Each row of U corresponds to a new time point, and
  28. // U must have length(T) rows. The input T is the time vector which
  29. // must be regularly spaced.
  30. //
  31. // Note: Two matrices are returned in a list.
  32. //
  33. //       G.x = X values in the plot.
  34. //       G.y = Y values in the plot.
  35. //
  36. // Copyright(C), by Jeffrey B. Layton, 1994
  37. // Version JBL 940906
  38. //---------------------------------------------------------------------------------
  39.  
  40. rfile abcdchk
  41. rfile tfchk
  42. rfile tf2ss
  43. rfile c2d
  44. rfile ltitr
  45.  
  46. lsim = function(a,b,c,d,u,t,x0)
  47. {
  48.    local(nargs,Dum,num,den,A,B,C,D,msg,estr,P,G,ns,n,x,y,X0)
  49.  
  50. // Count number of inputs
  51.    nargs=0;
  52.    if (exist(a)) {nargs=nargs+1;}
  53.    if (exist(b)) {nargs=nargs+1;}
  54.    if (exist(c)) {nargs=nargs+1;}
  55.    if (exist(d)) {nargs=nargs+1;}
  56.    if (exist(u)) {nargs=nargs+1;}
  57.    if (exist(t)) {nargs=nargs+1;}
  58.    if (exist(x0)) {nargs=nargs+1;}
  59.  
  60.    if (nargs < 4) {
  61.        error("LSIM: Incorrect number of input arguments.");
  62.    }
  63.  
  64. // Check for T.F. case
  65.    if (nargs == 4) {
  66. // Check if T.F. is proper
  67.        Dum=tfchk(a,b);
  68.        num=Dum.numc;
  69.        den=Dum.denc;
  70.        u=c;
  71.        t=d;
  72. // Convert T.F. to S.S.
  73.        Dum=tf2ss(num,den);
  74.        A=Dum.a;
  75.        B=Dum.b;
  76.        C=Dum.c;
  77.        D=Dum.d;
  78.    else
  79.        A=a;
  80.        B=b;
  81.        C=c;
  82.        D=d;
  83.        msg="";
  84.        msg=abcdchk(A,B,C,D);
  85.        if (msg != "") {
  86.            estr="lsim: "+msg;
  87.            error(estr);
  88.        }
  89.    }
  90.  
  91. // If X0 doesn't exist, then create a zero vector for it
  92.    if ( (nargs == 6) || (nargs == 4) ) {
  93.        X0 = zeros(A.nr,1);
  94.    else
  95.        X0=x0;
  96.    }
  97.  
  98. // Perform dimension checks
  99.    if (u.nc != D.nc) {
  100.        error("LSIM: U has to have same number of cols. as inputs");
  101.    }
  102.    if (u.nr != length(t)) {
  103.        error("LSIM: U must have same number of rows as length of T");
  104.    }
  105.  
  106. // Convert continuous to discrete using zero order hold.
  107.    Dum=c2d(A,B,t[2]-t[1]);
  108.    P=Dum.phi;
  109.    G=Dum.gamma;
  110.  
  111. // Propagate the simulation over the number of elements in t
  112.    x = ltitr(P,G,u,X0);
  113.    y = x * C.' + u * D.';
  114.  
  115. // Plot the results
  116.    xlabel("Time (secs)");
  117.    ylabel("Amplitude");
  118.    plot([t,y]);
  119.  
  120.    return << x=x; y=y >>
  121. };
  122.