[cacsd] fix demo Graphics => Technical plots => bode
[scilab.git] / scilab / modules / cacsd / macros / bode.sci
index 27fc1c2..a9125dc 100644 (file)
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C)  1985-2010 - INRIA - Serge Steer
-// This file must be used under the terms of the CeCILL.
-// This source file is licensed as described in the file COPYING, which
-// you should have received as part of this distribution.  The terms
-// are also available at
-// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
-
-function []=bode(varargin)
-  rhs=size(varargin)
-  if type(varargin($))==10 then
-    comments=varargin($),rhs=rhs-1;
-  else
-    comments=[];
-  end
-  fname="bode";//for error messages
-  fmax=[]
-  discr=%f //for shannon limit
-  if or(typeof(varargin(1))==["state-space" "rational"]) then
-    //sys,fmin,fmax [,pas] or sys,frq
-    refdim=1 //for error message
-    discr=varargin(1).dt<>'c';
-    if rhs==1 then //sys
-      [frq,repf]=repfreq(varargin(1),1d-3,1d3)
-    elseif rhs==2 then //sys,frq
-      if size(varargin(2),2)<2 then
-        error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
-                       fname,2,1))
-      end
-      [frq,repf]=repfreq(varargin(1:rhs))
-    elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
-      [frq,repf]=repfreq(varargin(1:rhs))
+// Copyright (C) 1985 - 2016 - INRIA - Serge Steer
+// Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2018 - Samuel GOUGEON
+//
+// This file is hereby licensed under the terms of the GNU GPL v2.0,
+// pursuant to article 5.3.4 of the CeCILL v.2.1.
+// This file was originally licensed under the terms of the CeCILL v2.1,
+// and continues to be available under such terms.
+// For more information, see the COPYING file which you should have received
+// along with this program.
+
+function bode(varargin)
+    rhs = size(varargin)
+
+    if rhs == 0 then
+        s = poly(0, "s");
+        h1 = syslin("c", (s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01));
+        num = 22801+4406.18*s+382.37*s^2+21.02*s^3+s^4;
+        den = 22952.25+4117.77*s+490.63*s^2+33.06*s^3+s^4;
+        h2 = syslin("c", num/den);
+
+        bode([h1; h2], 0.01, 100, ["h1"; "h2"]);
+        return;
+    end
+
+    rad = %f;
+    if type(varargin($)) == 10 then
+        if varargin($) == "rad" then
+            rad = %t;
+            rhs = rhs-1;
+            if type(varargin($-1)) == 10 then
+                comments = varargin($-1);
+                rhs = rhs-1;
+            else
+                comments = [];
+            end
+        else
+            comments = varargin($);
+            rhs = rhs-1;
+        end
     else
-      error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
+        comments = [];
     end
-    [phi,d]=phasemag(repf)
-    if rhs>=3 then fmax=varargin(3),end
-  elseif  type(varargin(1))==1 then
-    //frq,db,phi [,comments] or frq, repf [,comments]
-    refdim=2
-    select rhs
-    case 2 then //frq,repf
-      frq=varargin(1);
-      if size(frq,2)<2 then
-        error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
-                       fname,1,1))
-      end
-      if size(frq,2)<>size(varargin(2),2) then
-        error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
-                       fname,1,2))
-      end
-      [phi,d]=phasemag(varargin(2))
-    case 3 then  //frq,db,phi
-      [frq,d,phi]=varargin(1:rhs)
-      if size(frq,2)<>size(d,2) then
-        error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
-                       fname,1,2))
-      end
-      if size(frq,2)<>size(phi,2) then
-        error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
-                       fname,1,3))
-      end
+    fname = "bode"; // For error messages
+    fmax = [];
+    discr = %f; // For Shannon limit
+    if or(typeof(varargin(1)) == ["state-space" "rational" "zpk"]) then
+        // sys, fmin, fmax [,pas] or sys, frq
+        refdim = 1; // for error message
+        discr = varargin(1).dt<>"c";
+        if rhs == 1 then // sys
+            [frq, repf] = repfreq(varargin(1), 1d-3, 1d3);
+        elseif rhs == 2 then // sys, frq
+            if size(varargin(2), 2) < 2 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"), ..
+                fname, 2, 1))
+            end
+            [frq, repf] = repfreq(varargin(1:rhs));
+        elseif or(rhs == (3:4)) then // sys, fmin, fmax [,pas]
+            [frq, repf] = repfreq(varargin(1:rhs));
+        else
+            error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 1, 5))
+        end
+        [phi, d] = phasemag(repf);
+        if rhs >= 3 then fmax = varargin(3); end
+    elseif type(varargin(1)) == 1 then
+        // frq, db, phi [,comments] or frq, repf [,comments]
+        refdim = 2;
+        select rhs
+        case 2 then // frq,repf
+            frq = varargin(1);
+            if size(frq, 2) < 2 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"), ..
+                fname, 1, 1))
+            end
+            if size(frq, 2) <> size(varargin(2), 2) then
+                error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"), ..
+                fname, 1, 2))
+            end
+            [phi, d] = phasemag(varargin(2));
+        case 3 then  // frq, db, phi
+            [frq, d, phi] = varargin(1:rhs);
+            if size(frq, 2) <> size(d, 2) then
+                error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"), ..
+                fname, 1, 2))
+            end
+            if size(frq, 2) <> size(phi, 2) then
+                error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"), ..
+                fname, 1, 3))
+            end
+        else
+            error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 2, 4))
+        end
     else
-      error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
+        ierr=execstr("%"+overloadname(varargin(1))+"_bode(varargin(:))","errcatch")
+        if ierr<>0 then
+            error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
+        end
+        return
     end
-  else
-    error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
-  end;
-  frq=frq';d=d',phi=phi'
-  [n,mn]=size(d)
-
-  if comments==[] then
-    hx=0.48
-  else
-    if size(comments,"*")<>mn then
-      error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),...
-                     fname,refdim,rhs+1))
+    frq = frq';
+    d = d';
+    phi = phi';
+    [n, mn] = size(d);
+
+    if comments <> [] & size(comments, "*") <> mn then
+        msg = _("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n")
+        error(mprintf(msg, fname, refdim, rhs+1))
     end
-    hx=0.43
-  end;
-
-  drawlater()
-  sciCurAxes=gca();
-  axes=sciCurAxes;
-  wrect=axes.axes_bounds;
-
-
-  //magnitude
-  axes.axes_bounds=[wrect(1)+0,wrect(2)+0,wrect(3)*1.0,wrect(4)*hx*0.95]
-  axes.data_bounds = [min(frq),min(d);max(frq),max(d)];
-  axes.log_flags = "lnn" ;
-  axes.grid=color("lightgrey")*ones(1,3);
-  axes.axes_visible="on";
-  axes.clip_state = "clipgrf";
-  if size(d,2)>1&size(frq,2)==1 then
-    xpolys(frq(:,ones(1,mn)),d,1:mn)
-  else
-    xpolys(frq,d,1:mn)
-  end
-  //set datatips info
-  e=gce();
-
-  for i=1:size(e.children,"*")
-    datatipInitStruct(e.children(i),"formatfunction","formatBodeMagTip")
-  end
-
-  if discr&fmax<>[]&max(frq)<fmax then
-    xpoly(max(frq)*[1;1],axes.y_ticks.locations([1 $]));e=gce();
-    e.foreground=5;
-  end
-  xtitle("",_("Frequency (Hz)"),_("Magnitude (dB)"));
-
-  //phase
-
-  axes=newaxes();
-  axes.axes_bounds=[wrect(1)+0,wrect(2)+wrect(4)*hx,wrect(3)*1.0,wrect(4)*hx*0.95];
-  axes.data_bounds = [mini(frq),mini(phi);maxi(frq),maxi(phi)];
-  axes.log_flags = "lnn" ;
-  axes.grid=color("lightgrey")*ones(1,3);
-  axes.axes_visible="on";
-  axes.clip_state = "clipgrf";
-  if size(phi,2)>1&size(frq,2)==1 then
-    xpolys(frq(:,ones(1,mn)),phi,1:mn)
-  else
-    xpolys(frq,phi,1:mn)
-  end
-  ephi=gce()
-  //set datatips info
-  for i=1:size(ephi.children,"*")
-    datatipInitStruct(ephi.children(i),"formatfunction","formatBodePhaseTip")
-  end
-
-  if discr&fmax<>[]&max(frq)<fmax then
-    xpoly(max(frq)*[1;1],axes.y_ticks.locations([1 $]));e=gce();
-    e.foreground=5;
-  end
-  xtitle("",_("Frequency (Hz)"),_("Phase (degree)"));
-  // create legend
-  if comments<>[] then
-    c=captions(ephi.children,comments,"lower_caption")
-    c.background=get(gcf(),"background")
-  end
-  drawnow()
-  // return to the previous scale
-  set( "current_axes", sciCurAxes ) ;
 
-endfunction
-function str=formatBodeMagTip(curve,pt,index)
-//this function is called by the datatips mechanism to format the tip
-//string for the magnitude bode curves
-  str=msprintf("%.4g"+_("Hz")+"\n%.4g"+_("dB"), pt(1),pt(2))
-endfunction
-function str=formatBodePhaseTip(curve,pt,index)
-//this function is called by the datatip mechanism to format the tip
-//string for the bode phase curves
-  str=msprintf("%.4g"+_("Hz")+"\n %.4g"+_("dB"), pt(1),pt(2))
+    fig = gcf();
+    immediate_drawing = fig.immediate_drawing;
+    fig.immediate_drawing = "off";
+
+    sciCurAxes = gca();
+    axes = sciCurAxes;
+    wrect = axes.axes_bounds;
+
+    // Magnitude plot
+    // --------------
+    axes.axes_bounds = [wrect(1), wrect(2), wrect(3), wrect(4)*0.5];
+    axes.data_bounds = [min(frq), min(d); max(frq), max(d)];
+    axes.log_flags = "lnn";
+    axes.grid = color("lightgrey")*ones(1, 3);
+    axes.grid_style = [8 8];
+    axes.sub_ticks(1) = 8;
+    axes.axes_visible = "on";
+    axes.clip_state = "clipgrf";
+    if size(d, 2) > 1 & size(frq, 2) == 1 then
+        xpolys(frq(:, ones(1, mn)), d, 1:mn);
+    else
+        xpolys(frq, d, 1:mn);
+    end
+    // Set datatips info
+    e = gce();
+
+    for i=1:size(e.children, "*")
+        e.children(i).display_function = "formatBodeMagTip"
+    end
+
+    if discr & fmax <> [] & max(frq) < fmax then
+        xpoly(max(frq)*[1; 1], axes.y_ticks.locations([1 $]));
+        e = gce();
+        e.foreground = 5;
+    end
+    xtitle("", _("Frequency (Hz)"), _("Magnitude (dB)"));
+    axesM = axes;
+
+    // Phase plot
+    // ----------
+    axes = newaxes();
+    axes.axes_bounds = [wrect(1), wrect(2)+wrect(4)*0.5, wrect(3), wrect(4)*0.5];
+    axes.data_bounds = [min(frq), min(phi); max(frq), max(phi)];
+    axes.log_flags = "lnn";
+    axes.grid = color("lightgrey")*ones(1, 3);
+    axes.grid_style = [8 8];
+    axes.sub_ticks(1) = 8;
+    axes.axes_visible = "on";
+    axes.clip_state = "clipgrf";
+    if size(phi, 2) > 1 & size(frq, 2) == 1 then
+        xpolys(frq(:, ones(1, mn)), phi, 1:mn);
+    else
+        xpolys(frq, phi, 1:mn);
+    end
+    ephi = gce();
+    // Set datatips info
+    for i=1:size(ephi.children, "*")
+        ephi.children(i).display_function = "formatBodePhaseTip";
+    end
+
+    if discr & fmax <> [] & max(frq) < fmax then
+        xpoly(max(frq)*[1; 1], axes.y_ticks.locations([1 $]));
+        e = gce();
+        e.foreground = 5;
+    end
+    xtitle("", _("Frequency (Hz)"), _("Phase (degree)"));
+
+    // Create legend
+    // -------------
+    if comments <> [] then
+        c = captions(ephi.children, comments, "lower_caption");
+        c.background = get(gcf(), "background");
+    end
+
+    // Rendering
+    fig.immediate_drawing = immediate_drawing;
+
+    // Post-tuning the vertical distribution of Bode plots
+    // ---------------------------------------------------
+    // (only after rendering, otherwise the legend position is [0 0] (bug))
+    if comments <> [] then
+        LgH = (1-c.position(2))*wrect(4)*0.5; // Legend height
+        // Ah = Magnitude axes height.
+        // Ah + Ah+LgH*Ah/0.5 = wrect(4)  => Ah (1 + 1 + 2*LgH) = wrect(4)
+        Ah = wrect(4)/(2+2*LgH);
+        axes.axes_bounds([2 4]) = [wrect(2)+Ah, wrect(4)-Ah];
+        axesM.axes_bounds(4) = Ah*1.05;  // Magnitude plot
+    end
+
+    // Returns to the initial axes
+    sca(sciCurAxes);
+
+    // rad/s mode
+    // ----------
+    if rad == %t then
+        // This function modifies the Bode diagrams for a rad/s display instead of Hz.
+        // h is a hanlde of a figure containing Bode diagrams.
+        // Ref: http://forge.scilab.org/index.php/p/cpge/source/tree/HEAD/macros/bode_Hz2rad_2.sci
+        labels = [_("Phase (degree)"); _("Magnitude (dB)")];
+        pos_h = [9, 5];
+        for k=1:2
+            for i=1:size(fig.children(k).children, 1)
+                if fig.children(k).children(i).type == "Compound"
+                    for j=1:size(fig.children(k).children(i).children, 1)
+                        fig.children(k).children(i).children(j).data(:, 1) = fig.children(k).children(i).children(j).data(:, 1)*2*%pi;
+                    end
+
+                    // h.children(k).title.text = h.children(k).y_label.text;
+                    xmin1 = fig.children(k).data_bounds(1)*2*%pi;
+                    xmax1 = fig.children(k).data_bounds(2)*2*%pi;
+                    ymin1 = fig.children(k).data_bounds(3);
+                    ymax1 = fig.children(k).data_bounds(4);
+
+                    rect = [xmin1, ymin1, xmax1, ymax1];
+                    nb_dec = log(xmax1/xmin1)/log(10);
+                    fig.children(k).x_label.text = _("Frequency (rad/s)");
+                    fig.children(k).x_location = "bottom";
+                    fig.children(k).y_label.text = labels(k);
+                    replot(rect, fig.children(k));
+                end
+            end
+        end
+    end
 endfunction