scilab/modules/cacsd/macros/bode.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C)  1985-2010 - INRIA - Serge Steer
   3 // This file must be used under the terms of the CeCILL.
   4 // This source file is licensed as described in the file COPYING, which
   5 // you should have received as part of this distribution.  The terms
   6 // are also available at
   7 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
   8
   9 function []=bode(varargin)
  10     rhs=size(varargin)
  11
  12     if rhs == 0 then
  13         s=poly(0,"s")
  14         h1=syslin("c",(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01))
  15         num=22801+4406.18*s+382.37*s^2+21.02*s^3+s^4;
  16         den=22952.25+4117.77*s+490.63*s^2+33.06*s^3+s^4
  17         h2=syslin("c",num/den);
  18
  19         bode([h1;h2],0.01,100,["h1";"h2"])
  20         return;
  21     end
  22
  23     if type(varargin($))==10 then
  24         comments=varargin($),rhs=rhs-1;
  25     else
  26         comments=[];
  27     end
  28     fname="bode";//for error messages
  29     fmax=[]
  30     discr=%f //for shannon limit
  31     if or(typeof(varargin(1))==["state-space" "rational"]) then
  32         //sys,fmin,fmax [,pas] or sys,frq
  33         refdim=1 //for error message
  34         discr=varargin(1).dt<>"c";
  35         if rhs==1 then //sys
  36             [frq,repf]=repfreq(varargin(1),1d-3,1d3)
  37         elseif rhs==2 then //sys,frq
  38             if size(varargin(2),2)<2 then
  39                 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
  40                 fname,2,1))
  41             end
  42             [frq,repf]=repfreq(varargin(1:rhs))
  43         elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
  44             [frq,repf]=repfreq(varargin(1:rhs))
  45         else
  46             error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
  47         end
  48         [phi,d]=phasemag(repf)
  49         if rhs>=3 then fmax=varargin(3),end
  50     elseif  type(varargin(1))==1 then
  51         //frq,db,phi [,comments] or frq, repf [,comments]
  52         refdim=2
  53         select rhs
  54         case 2 then //frq,repf
  55             frq=varargin(1);
  56             if size(frq,2)<2 then
  57                 error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
  58                 fname,1,1))
  59             end
  60             if size(frq,2)<>size(varargin(2),2) then
  61                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
  62                 fname,1,2))
  63             end
  64             [phi,d]=phasemag(varargin(2))
  65         case 3 then  //frq,db,phi
  66             [frq,d,phi]=varargin(1:rhs)
  67             if size(frq,2)<>size(d,2) then
  68                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
  69                 fname,1,2))
  70             end
  71             if size(frq,2)<>size(phi,2) then
  72                 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
  73                 fname,1,3))
  74             end
  75         else
  76             error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
  77         end
  78     else
  79         error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
  80     end;
  81     frq=frq';d=d',phi=phi'
  82     [n,mn]=size(d)
  83
  84     if comments==[] then
  85         hx=0.48
  86     else
  87         if size(comments,"*")<>mn then
  88             error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),...
  89             fname,refdim,rhs+1))
  90         end
  91         hx=0.43
  92     end;
  93
  94     fig=gcf();
  95     immediate_drawing=fig.immediate_drawing;
  96     fig.immediate_drawing="off";
  97
  98     sciCurAxes=gca();
  99     axes=sciCurAxes;
 100     wrect=axes.axes_bounds;
 101
 102
 103     //magnitude
 104     axes.axes_bounds=[wrect(1)+0,wrect(2)+0,wrect(3)*1.0,wrect(4)*hx*0.95]
 105     axes.data_bounds = [min(frq),min(d);max(frq),max(d)];
 106     axes.log_flags = "lnn" ;
 107     axes.grid=color("lightgrey")*ones(1,3);
 108     axes.axes_visible="on";
 109     axes.clip_state = "clipgrf";
 110     if size(d,2)>1&size(frq,2)==1 then
 111         xpolys(frq(:,ones(1,mn)),d,1:mn)
 112     else
 113         xpolys(frq,d,1:mn)
 114     end
 115     //set datatips info
 116     e=gce();
 117
 118     for i=1:size(e.children,"*")
 119         datatipInitStruct(e.children(i),"formatfunction","formatBodeMagTip")
 120     end
 121
 122     if discr&fmax<>[]&max(frq)<fmax then
 123         xpoly(max(frq)*[1;1],axes.y_ticks.locations([1 $]));e=gce();
 124         e.foreground=5;
 125     end
 126     xtitle("",_("Frequency (Hz)"),_("Magnitude (dB)"));
 127
 128     //phase
 129
 130     axes=newaxes();
 131     axes.axes_bounds=[wrect(1)+0,wrect(2)+wrect(4)*hx,wrect(3)*1.0,wrect(4)*hx*0.95];
 132     axes.data_bounds = [min(frq),min(phi);max(frq),max(phi)];
 133     axes.log_flags = "lnn" ;
 134     axes.grid=color("lightgrey")*ones(1,3);
 135     axes.axes_visible="on";
 136     axes.clip_state = "clipgrf";
 137     if size(phi,2)>1&size(frq,2)==1 then
 138         xpolys(frq(:,ones(1,mn)),phi,1:mn)
 139     else
 140         xpolys(frq,phi,1:mn)
 141     end
 142     ephi=gce()
 143     //set datatips info
 144     for i=1:size(ephi.children,"*")
 145         datatipInitStruct(ephi.children(i),"formatfunction","formatBodePhaseTip")
 146     end
 147
 148     if discr&fmax<>[]&max(frq)<fmax then
 149         xpoly(max(frq)*[1;1],axes.y_ticks.locations([1 $]));e=gce();
 150         e.foreground=5;
 151     end
 152     xtitle("",_("Frequency (Hz)"),_("Phase (degree)"));
 153     // create legend
 154     if comments<>[] then
 155         c=captions(ephi.children,comments,"lower_caption")
 156         c.background=get(gcf(),"background")
 157     end
 158     fig.immediate_drawing=immediate_drawing;
 159     // return to the previous scale
 160     set( "current_axes", sciCurAxes ) ;
 161
 162 endfunction
 163
 164 function str=formatBodeMagTip(curve,pt,index)
 165     //this function is called by the datatips mechanism to format the tip
 166     //string for the magnitude bode curves
 167     str=msprintf("%.4g"+_("Hz")+"\n%.4g"+_("dB"), pt(1),pt(2))
 168 endfunction
 169
 170 function str=formatBodePhaseTip(curve,pt,index)
 171     //this function is called by the datatip mechanism to format the tip
 172     //string for the bode phase curves
 173     str=msprintf("%.4g"+_("Hz")+"\n %.4g"+"°", pt(1),pt(2))
 174 endfunction