home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / damp < prev    next >
Encoding:
Text File  |  1995-11-15  |  2.0 KB  |  74 lines
  1. //---------------------------------------------------------------------------
  2. //
  3. // damp
  4. //
  5. // Syntax: s=damp(a)
  6. //
  7. //
  8. // This routine computes the natural frequency and damping for
  9. // continuous systems. It prints a table of the eigenvalues of
  10. // the matrix a, the associated damping factors, the associated natural
  11. // frequency (rad/s and Hz.).
  12. //
  13. // Also, the routine returns the natrural frequencies and damping
  14. // factors in a list:
  15. //
  16. //  a.omegan = Natural Frequencies (rad/s)
  17. //  a.zeta  = Damping Factors
  18. //
  19. // Copyright (C), by Jeffrey B. Layton, 1994
  20. // Version JBL 940918
  21. //---------------------------------------------------------------------------
  22.  
  23. rfile roots
  24. rfile esort
  25.  
  26. damp = function(a)
  27. {
  28.    local(dum1,omegan,eigvlsort,omegahz,zeta,i)
  29.  
  30. // Compute the eigenvalues of the matrix and sort them
  31.    if (a.nr == a.nc) {
  32.        dum1=esort(eig(a).val');
  33.        eigvlsort=dum1.s;
  34.    else if (a.nr == 1) {
  35.        dum1=esort(roots(a));
  36.        eigvlsort=dum1.s;
  37.    else if (a.nc == 1) {
  38.        eigvlsort=a;
  39.    else
  40.        error("damp: Must be a vector or a sqaure matrix.");
  41.    }}}
  42.  
  43. // Compute the natural frequency on rad/s.
  44.    omegan=abs(eigvlsort);
  45.  
  46. // Compute the damping factor
  47.    zeta=-cos(atan2(imag(eigvlsort),real(eigvlsort)));
  48.  
  49. // Compute the natural frequency (Hz).
  50.    omegahz=sqrt(eigvlsort)/(2.0*pi);
  51.  
  52.    printf("%s"," \n");
  53.    printf("%s","     Eigenvalue            Damping     Freq. (rad/s)    Freq. (Hz). \n");
  54.    printf("%s","     ----------            -------     -------------    -----------\n");
  55.    for (i in 1:omegan.nr) {
  56.         if (imag(eigvlsort[i]) < 0.0) {
  57.             printf("%-10.6f -j%-10.6f", real (eigvlsort[i]), abs(imag (eigvlsort[i])));
  58.         else
  59.             printf("%-10.6f  j%-10.6f", real (eigvlsort[i]), imag (eigvlsort[i]));
  60.         }
  61.         printf("%s","   ");
  62.         printf("%-10.6f",zeta[i]);
  63.         printf("%s","     ");
  64.         printf("%-10.6f",omegan[i]);
  65.         printf("%s","      ");
  66.         printf("%-10.6f",omegahz[i]);
  67.         printf("%s"," \n");
  68.    }
  69.    printf("%s"," \n");
  70.    
  71.    return << omegan=omegan; zeta=zeta >>
  72.  
  73. };
  74.