share/contrib/bode.mac

   1 /* bode.mac -- functions to draw Bode gain and phase plots
   2  *
   3  * copyright Robert Dodier, October 2005
   4  * Released under the terms of the GNU General Public License
   5  *
   6  * Examples (1 through 7 from http://www.swarthmore.edu/NatSci/echeeve1/Ref/Bode/BodeHow.html, 8 from Ron Crummett)
   7  *
   8    H1 (s) := 100 * (1 + s) / ((s + 10) * (s + 100));
   9    H2 (s) := 1 / (1 + s/omega0);
  10    H3 (s) := 1 / (1 + s/omega0)^2;
  11    H4 (s) := 1 + s/omega0;
  12    H5 (s) := 1/s;
  13    H6 (s) := 1/((s/omega0)^2 + 2 * zeta * (s/omega0) + 1);
  14    H7 (s) := (s/omega0)^2 + 2 * zeta * (s/omega0) + 1;
  15    H8 (s) := 0.5 / (0.0001 * s^3 + 0.002 * s^2 + 0.01 * s);
  16
  17    bode_gain (H1 (s), [w, 1/1000, 1000]);
  18    bode_phase (H1 (s), [w, 1/1000, 1000]);
  19    bode_gain (H2 (s), [w, 1/1000, 1000]), omega0 = 10;
  20    bode_phase (H2 (s), [w, 1/1000, 1000]), omega0 = 10;
  21    bode_gain (H3 (s), [w, 1/1000, 1000]), omega0 = 10;
  22    bode_phase (H3 (s), [w, 1/1000, 1000]), omega0 = 10;
  23    bode_gain (H4 (s), [w, 1/1000, 1000]), omega0 = 10;
  24    bode_phase (H4 (s), [w, 1/1000, 1000]), omega0 = 10;
  25    bode_gain (H5 (s), [w, 1/1000, 1000]);
  26    bode_phase (H5 (s), [w, 1/1000, 1000]);                              <-- carg causes an asksign here (sigh)
  27    bode_gain (H6 (s), [w, 1/1000, 1000]), omega0 = 10, zeta = 1/10;
  28    bode_phase (H6 (s), [w, 1/1000, 1000]), omega0 = 10, zeta = 1/10;
  29    bode_gain (H7 (s), [w, 1/1000, 1000]), omega0 = 10, zeta = 1/10;
  30    bode_phase (H7 (s), [w, 1/1000, 1000]), omega0 = 10, zeta = 1/10;
  31    bode_gain (H8 (s), [w, 1/1000, 1000]);
  32    block ([bode_phase_unwrap : false], bode_phase (H8 (s), [w, 1/1000, 1000]));
  33    block ([bode_phase_unwrap : true], bode_phase (H8 (s), [w, 1/1000, 1000]));
  34  */
  35
  36 bode_plot2d: 'plot2d;
  37
  38 bode_phase_unwrap : false;
  39
  40 bode_grid : false;
  41
  42 log10 (x) := log (x) / log (10);
  43
  44 carg_unwrapped (z) := block ([a: carg (z)], charfun (a < 0) * (2*%pi + a) + (1 - charfun (a < 0)) * a);
  45
  46 bode_gain (H, range, [plot_opts]) := block ([omega, L, s, my_preamble],
  47   omega : first (range),
  48   L : block ([listdummyvars : false], listofvars (H)),
  49   if length (L) # 1
  50     then throw (oops (msg = "bode_gain: failed to identify a unique variable", expr = H, variables = L))
  51     else s : first (L),
  52
  53   my_preamble : concat ("set nokey; set title \"Bode gain plot for ", string (H), "\""),
  54   if bode_grid then my_preamble : concat (my_preamble, "; set grid"),
  55
  56   H : subst (%i * omega, s, H),
  57
  58   buildq ([plot_opts],
  59     bode_plot2d (10 * log10 (cabs (H * conjugate (H))), range,
  60     [logx, true], [gnuplot_preamble, my_preamble], splice (plot_opts))),
  61   ev (%%));
  62
  63 bode_phase (H, range, [plot_opts]) := block ([omega, L, s, my_preamble],
  64   omega : first (range),
  65   L : block ([listdummyvars : false], listofvars (H)),
  66   if length (L) # 1
  67     then throw (oops (msg = "bode_phase: failed to identify a unique variable", expr = H, variables = L))
  68     else s : first (L),
  69
  70   my_preamble : concat ("set nokey; set title \"Bode phase plot for ", string (H), "\""),
  71   if bode_grid then my_preamble : concat (my_preamble, "; set grid"),
  72
  73   H : subst (%i * omega, s, H),
  74
  75   if bode_phase_unwrap
  76     then buildq ([plot_opts],
  77       bode_plot2d (180/%pi * carg_unwrapped (H), range, [logx, true], [gnuplot_preamble, my_preamble], splice (plot_opts)))
  78     else buildq ([plot_opts],
  79       bode_plot2d (180/%pi * carg (H), range, [logx, true], [gnuplot_preamble, my_preamble], splice (plot_opts))),
  80   ev (%%));