home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / ddamp < prev    next >
Encoding:
Text File  |  1995-11-15  |  3.2 KB  |  109 lines
  1. //--------------------------------------------------------------------------------
  2. //
  3. // ddamp
  4. //
  5. // Syntax: A=ddamp(a,Ts)
  6. //
  7. // This routine compures the natural frequency and damping for discrete
  8. // systems. It prints a table of the eiganvalues of the matrix a, the
  9. // z-plane magnitude, the equivalent s-plane natural frequency (rad/s),
  10. // the equivalent s-plane frequency (Hz.), and the damping factors.
  11. //
  12. // The routine can is called as,
  13. //
  14. //     A=ddamp(a,Ts)
  15. //
  16. // where a is the input variable and Ts is the sample time. The variable a
  17. // can be in several forms:
  18. //
  19. //       1) If a is square, it is assumed to be the state-space
  20. //          "a" matrix.
  21. //       2) If a is a row vector, it is assumed to be a vector of
  22. //          the polynomial coefficients from a transfer function.
  23. //       3) If a is a column vector, it is assumed to contain
  24. //          root locations.
  25. //
  26. // It returns the s-plane equivalent natural frequency (rad/s), the damping
  27. // factors, and the z-plane magnitude. If the routine is called with only the
  28. // input variable a (the sample time Ts is not input), then ddamp will not
  29. // print the table and will return values for only the magnitude.
  30. //
  31. // For a discrete system eigenvalue, lambda, the equivalent s-plane
  32. // natural frequency and damping ratio are
  33. //
  34. //           Wn = abs(log(lambda))/Ts
  35. //
  36. //            Z = -cos(angle(log(lambda)))
  37. //
  38. //
  39. // The routine returns the natrural frequencies, damping, and the z-plane
  40. // magnitudes in a list:
  41. //
  42. //  A.wn  = Natural Frequencies (rad/s)
  43. //  A.z   = Damping Factors
  44. //  A.mag = z-plane magnitudes
  45. //
  46. // Copyright (C), by Jeffrey B. Layton, 1994
  47. // Version JBL 940917
  48. //--------------------------------------------------------------------------------
  49.  
  50. rfile dsort
  51. rfile roots
  52.  
  53. ddamp = function(a,Ts)
  54. {
  55.    local(narg,r,mag,wn,wnhz,z,i,s)
  56.  
  57. // Count number of arguments
  58.    narg=0;
  59.    if (exist(a)) { narg=narg+1; }
  60.    if (exist(Ts)) { narg=narg+1; }
  61.  
  62. // Start by sorting
  63.    if (a.nr == a.nc) {
  64.        r=dsort(eig(a).val).s;
  65.    else if (a.nr == 1) {
  66.        r=dsort(roots(a)).s;
  67.    else if (a.nc == 1) {
  68.        r=a;
  69.    else
  70.        error("Must be a vector or a Square matrix.");
  71.    }}}
  72.  
  73.    mag=abs(r);
  74.  
  75.    if (narg == 2) {   // If sample time is given solve for equivalent s-plane roots
  76.        s=log(r)/Ts;
  77.        wn=abs(s);
  78.        z=-cos(atan2(imag(s),real(s)));
  79.        wnhz=sqrt(s)/(2.0*pi);
  80.  
  81.        printf("%s"," \n");
  82.    printf("%s","     Eigenvalue              Mag       Damping     Freq. (rad/s)    Freq. (Hz). \n");
  83.    printf("%s","     ----------              ---       -------     ------------     ----------- \n");
  84.        for (i in 1:wn.nr) {
  85.             if (imag(s[i]) < 0.0) {
  86.                 printf("%-10.6f -j%-10.6f", real (s[i]), abs(imag (s[i])));
  87.             else
  88.                 printf("%-10.6f  j%-10.6f", real (s[i]), imag (s[i]));
  89.             }
  90.             printf("%s","   ");
  91.             printf("%-10.6f",mag[i]);
  92.             printf("%s","   ");
  93.             printf("%-10.6f",z[i]);
  94.             printf("%s","     ");
  95.             printf("%-10.6f",wn[i]);
  96.             printf("%s","      ");
  97.             printf("%-10.6f",wnhz[i]);
  98.             printf("%s"," \n");
  99.        }
  100.        printf("%s"," \n");
  101.    else
  102.        s=[];
  103.        wn=[];
  104.        z=[];
  105.    }
  106.    
  107.    return << wn=wn; z=z; mag=mag >>
  108. };
  109.