Massive indent of all codes:
[scilab.git] / 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-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