home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / displaypol < prev    next >
Encoding:
Text File  |  1995-11-15  |  3.0 KB  |  104 lines
  1. //-----------------------------------------------------------------------
  2. // displaypoles
  3. //
  4. // Syntax:
  5. //          poles=displaypoles(A,description)
  6. //
  7. // This routine computes the poles of A and their corresponding damping
  8. // ratio and damped circular frequency to the screen. The input, description,
  9. // can be used as a heading for the output.
  10. //
  11. // The routine returns the poles of A.
  12. // function poles = displaypoles(A,description,batch)
  13. // 
  14. // Originally Written by L. D. Peterson
  15. // Modified and ported to RLaB by J. B. Layton
  16. // Version JBL 940917
  17. //-----------------------------------------------------------------------
  18.  
  19. rfile banner
  20. rfile isstr
  21.  
  22. displaypoles = function(A,description)
  23. {
  24.    local(narg,iopt,str1,poles,cmplxpoles,realpoles,...
  25.          poles_zeta,poles_freq,j,Adum,poles_tc,outstr,line_count,...
  26.          outstring1,outstring2,outstr)
  27.  
  28. // Count number of input arguments
  29.    narg=0;
  30.    if (exist(A)) {narg=narg+1;}
  31.    if (exist(description)) {narg=narg+1;}
  32.  
  33.    if (narg == 1) {
  34.        iopt=0;
  35.        str1="";
  36.    else if (narg == 2) {
  37.        if (isstr(description)) {
  38.            iopt=0;
  39.        else
  40.            iopt=description;
  41.            str1="";
  42.        }
  43.    }}
  44.  
  45.    Adum=eig(A);              // poles of the design system
  46.    poles=Adum.val;
  47.    poles=sort(poles).val';
  48.    cmplxpoles=find(abs(imag(poles))>1.e-10);
  49.    realpoles=find(abs(imag(poles))<=1.e-10);
  50.    poles_zeta=zeros(size(poles));
  51.    poles_freq=zeros(size(poles));
  52.  
  53.    if (max(size(cmplxpoles)) != 0) {
  54.        poles_zeta[cmplxpoles] = ...
  55.    sin(atan(-real(poles[cmplxpoles])./abs(imag(poles[cmplxpoles]))));
  56.        poles_freq[cmplxpoles] = ...
  57.    -real(poles[cmplxpoles])./poles_zeta[cmplxpoles]./(2*pi);
  58.    }
  59.  
  60.    if (max(size(realpoles)) != 0) {
  61.        poles_zeta[realpoles]=ones(max(size(realpoles)),1);
  62.        poles_freq[realpoles]=-real(poles[realpoles])./(2*pi);
  63.    }
  64.    for (j in 1:max(size(poles))) {
  65.        if ((poles_zeta[j] != 0.) && (poles_freq[j] != 0.)) {
  66.             poles_tc[j]=((poles_zeta[j]*poles_freq[j]*2*pi)^(-1));
  67.        else
  68.             poles_tc[j] = 0.;
  69.        }
  70.    }
  71.  
  72.    if ( (narg >= 2) && (length(str1) != 0) ) {
  73.         printf("%s"," /n");
  74.         banner(str1+" Poles");
  75.         printf("%s"," /n");
  76.    }
  77.  
  78.    outstr="     re(rad/sec)   im(rad.sec)    zeta(%)";
  79.    printf("%s",outstr);
  80.    outstr="     f(hz)    time constant (sec) \n";
  81.    printf("%s",outstr);
  82.    outstr="---------------------------------------------------";
  83.    printf("%s",outstr);
  84.    outstr="-----------------------\n";
  85.    printf("%s",outstr);
  86.    line_count = 0;
  87.    for (j in 1:max(size(poles)) ) {
  88.         line_count=line_count+1;
  89.         sprintf(outstring1,"%3.0f  %12.4e %12.4e ",...
  90.                 j,real(poles[j]),imag(poles[j]));
  91.         if (poles_tc[j] != 0) {
  92.             sprintf(outstring2,"  %6.2f   %12.4e %12.4e",...
  93.                     poles_zeta[j]*100,poles_freq[j],poles_tc[j]);
  94.         else
  95.             sprintf(outstring2,"  %6.2f   %12.4e      --inf--",...
  96.                     poles_zeta[j]*100,poles_freq[j]);
  97.         }
  98.         outstr=outstring1+outstring2+" \n";
  99.         printf("%s",outstr);
  100.    }
  101.  
  102.    return poles;
  103. };
  104.