share/contrib/coma/coma.mac

   1 /* ======================================================================= */
   2 /* coma.mac                                                          V 2.1 */
   3 /* ----------------------------------------------------------------------- */
   4 /* COntrol engineering with MAxima            (c) Wilhelm Haager 2009-2019 */
   5 /* ----------------------------------------------------------------------- */
   6 /*                                                                         */
   7 /*  This program is free software; you can redistribute it and/or modify   */
   8 /*  it under the terms of the GNU General Public License as published by   */
   9 /*  the Free Software Foundation.                                          */
  10 /*                                                                         */
  11 /*  New:                                                                   */
  12 /*  - function bode_plot                                                   */
  13 /*  - option-processing for plot functions renewed                         */
  14 /*  - 2012-01-25: bugfix in "phase"                                        */
  15 /*  - 2012-02-27: "poles" and "zeros" enhanced                             */
  16 /*  - 2012-11-02: bugfix in "step_response": variable "t" local ()         */
  17 /*  - 2012-12-21: "ratsimp" replaced by "xthru" in many places             */
  18 /*  - 2013-01-05: "nilt" improved                                          */
  19 /*  - 2013-03-30: "pulse_response" added                                   */
  20 /*  - 2013-11-16: Rootsepsilon" set to 1.0e-7 (bug in Maxima 5.31)         */
  21 /*  - 2014-01-23: some bugfixes (nilt,...)                                 */
  22 /*  - 2014-07-14: margin adjustment in bode_plot                           */
  23 /*  - 2014-08-10: "xthru" replaced by "ratsimp" in "impedance_chain"       */
  24 /*  - 2014-09-07: "xthru" replaced by "ratsimp" in "transfer_function"     */
  25 /*  - 2014-11-07: Grid lines for Bode-plots                                */
  26 /*  - 2014-11-14: Nested lists for graphic options enabled                 */
  27 /*  - 2015-01-09: "//" - Parallel-Operator for electric resistances        */
  28 /*  - 2015-04-26: bug in "hurwitz" removed                                 */
  29 /*  - 2015-06-21: function "product_form" added                            */
  30 /*  - 2015-06-22: function "asymptotic" added, "phase" rewritten           */
  31 /*  - 2015-10-11: sum_form instead of standard_form; minor changes         */
  32 /*  - 2015-10-24  bugfixes, minor changes (product_form, plot)             */
  33 /*  - 2015-11-11  adaptations for wxMaxima 15.08 (wxplot_size)             */
  34 /*  - 2016-01-10  /_ infix operator and //_ postfix operator declared      */
  35 /*  - 2016-06-05  bug in _COMA_ppo removed (thanks to Daniel Volinski)     */
  36 /*  - 2016-07-25  sum_form enhanced; new function npartfrac                */
  37 /*  - 2016-08-05  sallrrots, spartfrac, silt added                         */
  38 /*  - 2016-09-19  bug in product_form removed (thanks to Klaus Weichinger) */
  39 /*  - 2017-02-20  new gain_crossover, new phase_crossover etc.             */
  40 /*  - 2017-05-07  bode_plot: option "logy=dB" added (P.Heidrich, S.Milani) */
  41 /*  - 2017-06-08  sum_form extended (second parameter: 1...6 )             */
  42 /*  - 2017-06-08  product_form extended (second parameter: 1...2)          */
  43 /*  - 2017-08-08  new bode_plots, phase, absval, asymptotic                */
  44 /*  - 2018-01-14  Nichols plot                                             */
  45 /*  - 2018-02-24  Weird error with cabs,carg,assume_pos avoided            */
  46 /*  - 2018-02-26  stable_area enhanced                                     */
  47 /*  - 2018-03-16  bugfix in product_form, gain_margin and time_delay       */
  48 /*  - 2018-12-03  (basic) treatment of time delays; function fullermap     */
  49 /*  - 2018-12-16  "nilt" extended, new variable "try_symbolic_ilt"         */
  50 /*  - 2019-01-03  "closed_loop" for time delayed systems                   */
  51 /*  - 2019-01-09  Digital algorithms added                                 */
  52 /*  - 2019-01-28  "normalize" and "tf" added                               */
  53 /*  - 2019-02-24  "stable_ares" depicted as transparent "region"           */
  54 /*  - 2019-03-26  "moving average" and "regan" added                       */
  55 /*  - 2019-05-08  "periodic" added                                         */
  56 /*  - 2019-05-21  "add_noise", "mult_noise" added                          */
  57 /* ======================================================================= */
  58
  59 disp("coma v.2.1 (Wilhelm Haager, 2019-05-21)")$
  60 load("draw");  /* not necessary any more but beneficial (2013-11-16) */
  61
  62 /* ======================================================================= */
  63 /* Some functions to avoid really weird errors (2018-02-24,2018-03-05)     */
  64 /* ======================================================================= */
  65 xabs(x):=lambda([u],cabs(num(u))/cabs(denom(u)))(xthru(x));
  66 xarg(x):=lambda([u],carg(num(u))-carg(denom(u)))(xthru(x));
  67
  68
  69 /* ======================================================================= */
  70 /* Default values                                                          */
  71 /* ======================================================================= */
  72
  73 fpprintprec:5;
  74 rootsepsilon:1.0e-7;
  75 coma_defaults:[grid=true, wx=true, dimensions=[700,400],
  76    user_preamble="set grid lt 3 lc '#111111',lt 3 lc '#AAAAAA';set tics font 'courier new,10'",
  77    color=[red,blue,dark-green,goldenrod,violet,gray40,dark-cyan,
  78          orange,brown,sea-green],
  79    fill_color=[red,blue,dark-green,goldenrod,violet,gray40,dark-cyan,
  80              orange,brown,sea-green]
  81    ];
  82
  83 /* ======================================================================= */
  84 /* Messages                                                                */
  85 /* ======================================================================= */
  86 _COMA_msg:[0,0,0];
  87 _COMA_msg[1]:"*** COMA error *** Coefficients must be numeric";
  88 _COMA_msg[2]:"*** COMA error *** yrange must be a list of TWO ranges";
  89 _COMA_msg[3]:"*** COMA warning *** Delay time ignored";
  90
  91 /* ======================================================================= */
  92 /* Plot Functions                                                          */
  93 /* ======================================================================= */
  94
  95 /* ----------------------------------------------------------------------- */
  96 /* plot - plots functions in one and two variables (2015-10-24)            */
  97 /* ----------------------------------------------------------------------- */
  98 plot(f,[opts]) :=
  99 block(
 100   [x,x1,x2,y,y1,y2,defs,hli,fli,h,n,t,wxx,varli,
 101     wxplot_size:assoc(dimensions,coma_defaults),ratprint:false],
 102   opts:flatten(opts),
 103   defs:[ ],
 104   [x1,x2]:assoc('xrange,opts,[0,1]),
 105   [y1,y2]:assoc('yrange,opts,[0,1]),
 106   varli:map(lambda([u],if _COMA_go(u) then [t] else listofvars(u)),flatten([f])),
 107   varli:map(lambda([u],delete(s,u)),varli),
 108   varli:map(lambda([u],if u=[] then u:[t] else u),varli), /* 10/2015 fuer Treppenfunktion */
 109   varli:map(sort,varli),
 110   n:apply(max,map(length,varli)),
 111   if n>2 then return(false),
 112   fli:flatten([f]),
 113   if n<2 then block(  /* 2d-plot */
 114     fli:map(lambda([u,w],
 115     if _COMA_go(u) then u else explicit(u,part(w,1),x1,x2)
 116     ),fli,varli))
 117   else block(  /* 3d-plot */
 118     fli:map(lambda([u,w],
 119     if _COMA_go(u) then u else explicit(u,part(w,1),x1,x2,part(w,2),y1,y2)
 120     ),fli,varli)),
 121   hli:_COMA_ppo(fli,defs,opts),
 122   if is(n>1) then block(
 123      if wxx#false then apply(wxdraw3d,hli) else apply(draw3d,hli))
 124   else block(
 125     if wxx#false then apply(wxdraw2d,hli) else apply(draw2d,hli))
 126 )$
 127
 128 /* ----------------------------------------------------------------------- */
 129 /* step_response(f,opts) - plots the step response of f                    */
 130 /* ----------------------------------------------------------------------- */
 131 step_response(f,[opts]) :=
 132 block(
 133   [wxplot_size:assoc(dimensions,coma_defaults),
 134     hli,defs:[],fli,t,tanf,tend,wxx,ratprint:false,assume_pos:false],
 135   opts:flatten(opts),
 136   if not listp(f) then f:[f],
 137   fli:flatten(map(lambda([_u], if _COMA_go(_u) then _u
 138         else if algorithmp(_u) then 'stepf(algorithm_step_response(_u))
 139         else (nilt(xthru(_u/s),s,t))),f)),
 140   fli:map(lambda([u],if _COMA_go(u) or algorithmp(u) then u
 141     /*    else explicit(if t>0 then u else 0,t,tanf,tend)),fli),  */
 142         else explicit(u,t,tanf,tend)),fli),
 143   hli:_COMA_ppo(fli,defs,opts),
 144   if option_exists('xrange,hli) then [tanf,tend]:get_option('xrange,hli)
 145   else block([tanf,tend]:[0,_COMA_npv(_COMA_srat(sublist(f,ntfp)),[2,5,10])],
 146        set_option('xrange=[tanf,tend],hli)),
 147   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 148 )$
 149 /* ----------------------------------------------------------------------- */
 150 /* step_response(f,opts) - plots the step response of f   (TESTVERSION 8.1.2019)                 */
 151 /* ----------------------------------------------------------------------- */
 152 step_response(f,[opts]) :=
 153 block(
 154   [wxplot_size:assoc(dimensions,coma_defaults),
 155     hli,defs:[],fli,t,tanf,tend,wxx,ratprint:false,assume_pos:false],
 156   opts:flatten(opts),
 157   if not listp(f) then f:[f],
 158   if option_exists('xrange,opts) then [tanf,tend]:get_option('xrange,opts)
 159   else block([tanf,tend]:[0,_COMA_npv(_COMA_srat(sublist(f,ntfp)),[2,5,10])],
 160        set_option('xrange=[tanf,tend],opts)),
 161   fli:flatten(map(lambda([_u], if _COMA_go(_u) then _u
 162     else if algorithmp(_u) then ('stepf(algorithm_step_response(_u,round(tend/ _u[2]))))
 163     else (nilt(    (_u/s),s,t))),f)),
 164   fli:map(lambda([u],if _COMA_go(u) then u
 165     /*    else explicit(if t>0 then u else 0,t,tanf,tend)),fli),  */
 166     else explicit(u,t,tanf,tend)),fli),
 167   hli:_COMA_ppo(fli,defs,opts),
 168   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 169 )$
 170 /* ----------------------------------------------------------------------- */
 171 /* response(f,opts) - plots the inverse Laplace transform of f (pulse resp)*/
 172 /* ----------------------------------------------------------------------- */
 173 response(f,[opts]) :=
 174 block(
 175   [wxplot_size:assoc(dimensions,coma_defaults),
 176     hli,defs:[],fli,t,tanf,tend,wxx,ratprint:false,assume_pos:false],
 177   opts:flatten(opts),
 178   f:flatten([f]),
 179   fli:flatten(map(lambda([u], if _COMA_go(u) then u
 180                              else nilt(xthru(u),s,t)),f)),
 181   fli:map(lambda([u],if _COMA_go(u) then u
 182         else explicit(if t>0 then u else 0,t,tanf,tend)),fli),
 183   hli:_COMA_ppo(fli,defs,opts),
 184   if option_exists('xrange,hli) then [tanf,tend]:get_option('xrange,hli)
 185   else block([tanf,tend]:[0,_COMA_npv(_COMA_srat(f),[2,5,10])],
 186        set_option('xrange=[tanf,tend],hli)),
 187   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 188 )$
 189
 190 /* ----------------------------------------------------------------------- */
 191 /* nyquist_plot(f,opts) - nyquist plot of the transfer function f          */
 192 /* ----------------------------------------------------------------------- */
 193 nyquist_plot(f,[opts]) :=
 194 block(
 195   [wxplot_size:assoc(dimensions,coma_defaults)*[1,1.25],
 196     hli,defs,t,tanf,tend,wxx,fli,ratprint:false],
 197   opts:flatten(opts),
 198   defs:[aspect_ratio=-1, nticks=500],
 199   [tanf,tend] : _COMA_bpr(sublist(flatten([f]),lambda([u],not _COMA_go(u)))),
 200   [tanf,tend] : [tanf/10,tend*10], /* ADAPT RANGE HERE */
 201   tanf:float(log(tanf)/log(10)),
 202   tend:float(log(tend)/log(10)),
 203   fli:map(lambda([u], if _COMA_go(u) then u else
 204           [float(realpart(ev(u,s=%i*t))), float(imagpart(ev(u,s=%i*t)))]
 205           ),flatten([f])),
 206   fli:map(lambda([u],if _COMA_go(u) then u else ev(u,t=10**t)),fli),
 207   fli:map(lambda([u],if _COMA_go(u) then u
 208           else parametric(part(u,1),part(u,2),t,tanf,tend)),fli),
 209   hli:ev(_COMA_ppo(fli,defs,opts)),
 210   if wxx#false then apply(wxdraw2d,hli) else apply(draw2d,hli)
 211 )$
 212
 213 /* ----------------------------------------------------------------------- */
 214 /* nichols_plot(f,opts) - Nichols plot of the Bode-diagram                 */
 215 /* ----------------------------------------------------------------------- */
 216 nichols_plot(f,[opts]) :=
 217 block(
 218   [wxplot_size:assoc(dimensions,coma_defaults),hli,defs,t,tanf,tend,fli,
 219        fli1,fli2,wxx,ratprint:false,db:false],
 220   opts:flatten(opts),
 221   if option_exists('logy,opts) and assoc('logy,opts)='dB then /* dB-scale, 2017-05-07 */
 222   (
 223     set_option('logy=false,opts),
 224     db:true
 225   ),
 226   fli:ratsimp(flatten([f])),
 227   [tanf,tend] : _COMA_bpr(sublist(flatten([f]),lambda([u],not _COMA_go(u)))),
 228   [tanf,tend] : [tanf/10,tend*10], /* ADAPT RANGE HERE */
 229   tanf:float(log(tanf)/log(10)),
 230   tend:float(log(tend)/log(10)),
 231   fli1:map(lambda([u],if _COMA_go(u) then u else
 232        float(if db then (20*log(absval(u))/log(10)) else (absval(u)))),fli),
 233   fli2:map(lambda([u], if _COMA_go(u) then u else ev(phase(u))),fli),
 234   fli1:map(lambda([u],if _COMA_go(u) then u else ev(u,omega=10**t)),fli1),
 235   fli2:map(lambda([u],if _COMA_go(u) then u else ev(u,omega=10**t)),fli2),
 236   fli:map(lambda([u,v],if _COMA_go(u) then u else
 237               parametric(u,v,t,tanf,tend)),fli2,fli1),
 238   defs:['nticks=500,'logy=true,grid=false,'user_preamble="set grid xtics mxtics ytics mytics"],
 239   hli:_COMA_ppo(fli,defs,opts),
 240   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 241 )$
 242
 243 /* ----------------------------------------------------------------------- */
 244 /* magnitude_plot(f,opts) - magnitude plot of the Bode-diagram             */
 245 /* ----------------------------------------------------------------------- */
 246 magnitude_plot(f,[opts]) :=
 247 block(
 248   [wxplot_size:assoc(dimensions,coma_defaults),hli,defs,t,tanf,tend,fli,wxx,ratprint:false,db:false],
 249     opts:flatten(opts),
 250   if option_exists('logy,opts) and assoc('logy,opts)='dB then /* dB-scale, 2017-05-07 */
 251   (
 252     set_option('logy=false,opts),
 253     db:true
 254   ),
 255   fli:ratsimp(flatten([f])),
 256   fli:map(lambda([u],if _COMA_go(u) then u else
 257          float(if db then (20*log(absval(u))/log(10)) else (absval(u)))),fli),
 258   fli:map(lambda([u],if _COMA_go(u) then u else explicit(u,omega,tanf,tend)),fli),
 259   defs:['logx=true,'logy=true,grid=false,'user_preamble="set grid xtics mxtics ytics mytics"],
 260   hli:_COMA_ppo(fli,defs,opts),
 261   if option_exists('xrange,opts) then [tanf,tend]:assoc('xrange,opts)
 262   else block([tanf,tend] : _COMA_bpr(f),
 263        set_option('xrange=[tanf,tend],hli)),
 264   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 265 )$
 266
 267 /* ----------------------------------------------------------------------- */
 268 /* phase_plot(f,opts) - phase plot of the Bode-diagram                     */
 269 /* ----------------------------------------------------------------------- */
 270 phase_plot(f,[opts]) :=
 271 block(
 272   [wxplot_size:assoc(dimensions,coma_defaults),
 273     hli,defs,t,tanf,tend,fli,ratprint:false,ph,wxx,omega],
 274   opts:flatten(opts),
 275   defs:['logx=true,grid=false,'user_preamble="set grid xtics mxtics ytics mytics",'nticks=500],
 276   fli:map(ratsimp,flatten([f])),
 277   fli:map(lambda([u], if _COMA_go(u) then u else ev(phase(u),assume_pos=true)),fli),
 278   fli:map(lambda([u], if _COMA_go(u) then u else explicit(u,omega,tanf,tend)),fli),
 279   hli:_COMA_ppo(fli,defs,opts),
 280   if option_exists('xrange,opts) then [tanf,tend]:assoc('xrange,opts)
 281   else block(
 282      [tanf,tend] : _COMA_bpr(f),
 283      set_option('xrange=[tanf,tend],hli)),
 284   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 285 )$
 286
 287 /* ----------------------------------------------------------------------- */
 288 /* bode_plot(f,opts) - Bode-diagram (magnitude_plot and phase_plot)        */
 289 /* ----------------------------------------------------------------------- */
 290 bode_plot(f,[opts]) :=
 291 block(
 292   [wxplot_size:assoc(dimensions,coma_defaults)*[1,2],
 293     fli,fli1,fli2,hli1,hli2,defs1,defs2,t,tanf,tend,ratprint:false,wxx,
 294     xdim,ydim,dims,g0,db:false],
 295   opts:flatten(opts),
 296   if option_exists('logy,opts) and assoc('logy,opts)='dB then /* dB-scale, 2017-05-07 */
 297   (
 298     set_option('logy=false,opts),
 299     db:true
 300   ),
 301   fli:ratsimp(flatten([f])),
 302   defs1:['logx=true,'logy=true,grid=false,'user_preamble=["set tmargin 0.5",
 303         "set lmargin 6","set rmargin 3","set bmargin 1.5","set grid xtics mxtics ytics mytics"]],
 304   defs2:['logx=true, 'nticks=500,grid=false,'user_preamble=["set tmargin 0.5",
 305         "set lmargin 6","set rmargin 3","set bmargin 1.5","set grid xtics mxtics ytics mytics"]],
 306   xdim:600, ydim:500, /* default values (for both diagrams together) */
 307   if listp(coma_defaults)=false then coma_defaults:['grid=true],
 308   if option_exists('dimensions,coma_defaults) then block(
 309      xdim:assoc('dimensions,coma_defaults)[1],
 310      ydim:2.0*assoc('dimensions,coma_defaults)[2]),
 311   if option_exists('dimensions,opts) then block(
 312      xdim:assoc('dimensions,opts)[1],
 313      ydim:assoc('dimensions,opts)[2]),
 314   dims:'dimensions=[xdim,ydim],
 315   /* Generating the function list */
 316   fli1:map(lambda([u],if _COMA_go(u) then u else
 317          float(if db then (20*log(absval(u))/log(10)) else (absval(u)))),fli),
 318   fli1:map(lambda([u],if _COMA_go(u) then u else explicit(u,omega,tanf,tend)),fli1),
 319   hli1:_COMA_ppo(fli1,defs1,opts),
 320   delete_option('dimensions,hli1),
 321   fli2:map(lambda([u], if _COMA_go(u) then u else ev(phase(u))),fli),
 322   fli2:map(lambda([u], if _COMA_go(u) then u else explicit(u,omega,tanf,tend)),fli2),
 323   hli2:_COMA_ppo(fli2,defs2,opts),
 324   delete_option('dimensions,hli2),
 325   /* Processing the ranges (xrange, yrange) */
 326   if option_exists('xrange,opts) then [tanf,tend]:assoc('xrange,opts)
 327   else block(
 328      [tanf,tend] : _COMA_bpr(f),
 329      set_option('xrange=[tanf,tend],hli1),
 330      set_option('xrange=[tanf,tend],hli2)),
 331   if option_exists('yrange,hli1) then block([yr],
 332      yr:assoc('yrange,opts),
 333      delete_option('yrange,hli1), delete_option('yrange,hli2),
 334      if listp(yr) and listp(yr[1]) and listp(yr[2]) then block(
 335         set_option('yrange=yr[1],hli1), set_option('yrange=yr[2],hli2))
 336      else disp(_COMA_msg[2])),
 337   /* Put the scenes together */
 338   g0:[apply(gr2d,ev(hli1)),apply(gr2d,ev(hli2)),dims],
 339   if wxx#false then apply(wxdraw,g0) else apply(draw,g0)
 340 )$
 341
 342 /* ----------------------------------------------------------------------- */
 343 /* absval - absolute value of a frequency response (treats asymptotic)     */
 344 /* ----------------------------------------------------------------------- */
 345 absval(f,[opts]) :=
 346 block([fli,ratprint:false,assume_pos:true,dt],
 347   /* --- treatment of asymptotic, 24.7.2017 --- */
 348   [fli,dt]:_COMA_strip(f),
 349   fli:map(lambda([u],
 350         if not mapatom(u) and op(u)='asymptotic then block
 351         ([assume_pos:true,inflag:true,f1,zz,nn,zl,nl,hz,hn,n,l0,l1],
 352           f1:product_form(first(args(u))),
 353           zz:num(f1),
 354           nn:denom(f1),
 355           zl:if not mapatom(zz) and op(zz)="*" then args(zz) else [zz],
 356           nl:if not mapatom(nn) and op(nn)="*" then args(nn) else [nn],
 357           hz:map(lambda([u],cl:coefficient_list(u,s),
 358           n:length(cl),
 359           [l0,l1]:[first(cl),last(cl)],
 360           if l0 # 0 then (if n<2 or omega<1/abs(l1)^(1/(n-1))
 361                           then l0 else s^(n-1)*l1)
 362                  else s^(n-1)*l1),zl),
 363           hn:map(lambda([u],cl:coefficient_list(u,s),
 364           n:length(cl),
 365           [l0,l1]:[first(cl),last(cl)],
 366           if l0 # 0 then (if n<2 or omega<1/abs(l1)^(1/(n-1))
 367                            then l0 else s^(n-1)*l1)
 368                   else s^(n-1)*l1),nl),
 369           apply("*",hz)/apply("*",hn)
 370         ) else u
 371       ),fli),
 372   /* --- END treatment of asymptotic --- */
 373   fli:float(map(xabs,subst(s=%i*omega,fli))),
 374   if length(fli) = 1 then fli[1] else fli)$
 375
 376 /* ----------------------------------------------------------------------- */
 377 /* phase(f) - phase of a frequency response (treats asymptotic)            */
 378 /* ----------------------------------------------------------------------- */
 379 phase(f):=block(
 380   [inflag:true,n,li,ass:false,res,ratprint:false,assume_pos:true,dt],
 381   [f,dt]:_COMA_strip(f),
 382   res:map(lambda([ff],
 383     if not mapatom(ff) and op(ff)='asymptotic then (ass:true, ff:product_form(first(args(ff))))
 384                                            else ff:product_form(ff),
 385     if not mapatom(ff) and op(ff)="*" then li:args(ff) else li:[ff],
 386     n:0,
 387     while n<length(li) do /* Zerlegen von Potenzen */
 388     (
 389       n:n+1,
 390       if not mapatom(li[n]) and op(li[n])="^" and abs(args(li[n])[2])>1 then
 391       block([prae,post,mid],
 392         prae:rest(li,-length(li)+n-1),
 393         post:rest(li,n),
 394         mid:makelist(args(li[n])[1]**signum(args(li[n])[2]),k,1,abs(args(li[n])[2])),
 395         li:append(prae,mid,post)
 396       )
 397     ),
 398     /* --- treatment of asymptotic 24.7.2017 --- */
 399     li:map(lambda([u],
 400         if ass then block([lu,lo,xcf:_COMA_cf(u)],
 401           if length(xcf)>0 then
 402           ( lu:limit(xarg(subst(s=%i*omega,u)),omega,0,plus),
 403             lo:limit(xarg(subst(s=%i*omega,u)),omega,inf),
 404             if omega<xcf[1] then lu else lo) else xarg(subst(s=%i*omega,u)) )
 405         else xarg(subst(s=%i*omega,u))
 406      ),li),
 407     /* --- END treatment of asymptotic --- */
 408     apply("+", li)*180/%pi),
 409     f) + log(dt)/s*omega*180/%pi,
 410   if length(res) = 1 then res[1] else res)$
 411
 412 /* ----------------------------------------------------------------------- */
 413 /* poles_and_zeros(f,opts) - plots the poles/zeros distribution            */
 414 /* ----------------------------------------------------------------------- */
 415 poles_and_zeros(f,[opts]) :=
 416 block(
 417   [wxplot_size:assoc(dimensions,coma_defaults),dt,
 418     hli,defs,zli,pli,wxx,x1,x2,y1,y2,range,ratprint:false],
 419   opts:flatten(opts),
 420   defs:[aspect_ratio=-1,points_joined=false],
 421   [f,dt]:_COMA_strip(f),     /* _COMA_(u)...Wrapper around graphic object */
 422   zli:map(lambda([u],if _COMA_go(u) then _COMA_(u) else zeros(u)),f),
 423   pli:map(lambda([u],if _COMA_go(u) then _COMA_(u) else poles(u)),f),
 424   range:_COMA_rlpr(flatten(sublist(append(zli,pli),listp)),1.2,1.2),
 425   [x1,x2]:assoc('xrange,opts,rhs(range[1])),
 426   [y1,y2]:assoc('yrange,opts,rhs(range[2])),
 427   set_option('xrange=[x1,x2],opts),
 428   set_option('yrange=[y1,y2],opts),
 429   zli:map(lambda([u],if op(u)=_COMA_ then u
 430           else map(lambda([v],[realpart(v),imagpart(v)]),u)),zli),
 431   zli:map(lambda([u],if op(u)=_COMA_ then u else [points(u)]),zli),
 432   zli:map(lambda([u],if op(u)=_COMA_ then u else cons(point_type=6,u)),zli),
 433   pli:map(lambda([u],if op(u)=_COMA_ then u
 434           else map(lambda([v],[realpart(v),imagpart(v)]),u)),pli),
 435   pli:map(lambda([u],if op(u)=_COMA_ then u else [points(u)]),pli),
 436   pli:map(lambda([u],if op(u)=_COMA_ then u else cons(point_type=2,u)),pli),
 437   pli:map(lambda([u,v],if op(u)=_COMA_ then args(u)
 438                 else append(u,v)),zli,pli),
 439      /*    else flatten([u,v])),zli,pli), */
 440   hli:_COMA_ppo(pli,defs,opts),
 441   hli:delete(points([]),flatten(hli)),
 442   if wxx#false then apply(wxdraw2d,hli) else apply(draw2d,hli)
 443 )$
 444
 445 /* ----------------------------------------------------------------------- */
 446 /* contourplot                                                             */
 447 /* ----------------------------------------------------------------------- */
 448 contourplot(f,p1,p2,[opts]):=
 449 block(
 450   [wxplot_size:assoc(dimensions,coma_defaults),
 451     x1,x2,y,y1,y2,defs,cli,hli,fli,ratprint:false,wxx],
 452   opts:flatten(opts),
 453   defs:[contours=[0],color=[red],line_width=[1]],
 454   hli:_COMA_ppo(1,defs,opts),
 455   [x1,x2]:assoc('xrange,opts,[0,1]),
 456   [y1,y2]:assoc('yrange,opts,[0,1]),
 457   cli:assoc('contours,hli),
 458   hli:map(lambda([u],f=u),cli),
 459   hli:map(lambda([u],implicit(u,p1,x1,x2,p2,y1,y2)),hli),
 460   hli:flatten(_COMA_ppo(hli,defs,opts)),
 461   delete_option('contours,hli),
 462   if wxx#false then apply(wxdraw2d,hli) else apply(draw2d,hli)
 463 )$
 464
 465 /* ----------------------------------------------------------------------- */
 466 /* root_locus(f,opts) - root locus plot of the transfer function f         */
 467 /* ----------------------------------------------------------------------- */
 468 root_locus(f,[opts]) :=
 469 block(
 470   [wxplot_size:assoc(dimensions,coma_defaults),
 471     i,range,nnticks,fac,pars,nn,pl,poli:[],gpoli:[],t,t1,t2,defs,hli,fli,h,
 472     ratprint:false,startpoints,endpoints,wxx,dt],
 473   opts:flatten(opts),
 474   defs:[aspect_ratio=-1,trange=[0.001,100],nticks=500],
 475   hli:_COMA_ppo(0,defs,opts),
 476   [t1,t2]:assoc('trange,hli),
 477   nnticks:assoc('nticks,hli),
 478   delete_option(nticks,hli),
 479   delete_option(trange,hli),
 480   fac:float((t2/t1)**(1/nnticks)),
 481   [fli,dt]:_COMA_strip(f),
 482   pars:map(lambda([u],first(delete(s,listofvars(u)))),fli), /* variable */
 483   range:_COMA_rlpr(
 484         flatten(map(lambda([u,v],poles(ev(u,ev(v)=t1))),fli,pars)),2,1.5),
 485   if not option_exists('xrange,hli) then set_option(first(range),hli),
 486   if not option_exists('yrange,hli) then set_option(last(range),hli),
 487   for i:1 step 1 thru length(fli) do block(
 488     nn:denom(ratsimp(fli[i])),
 489     poli:[zeros(ev(nn,ev(pars[i])=t1))],
 490     t:t1,
 491     while t<t2 do block(
 492       t:t*fac,
 493       poli:endcons(_COMA_rlsp(zeros(ev(nn,ev(pars[i])=t)),last(poli)),poli)),
 494     poli:map(lambda([v],map(lambda([u],[realpart(u),imagpart(u)]),v)),poli),
 495     [startpoints,endpoints]:[first(poli),last(poli)],
 496     poli:map('points,apply("[",args(transpose(apply(matrix,poli))))),
 497     poli:append(['points_joined=true,'point_type=-1],poli),
 498     poli:append(poli,['points_joined=false,'point_type=2,points(startpoints)]),
 499     poli:append(poli,['point_type=6,points(endpoints)]),
 500     gpoli:endcons(poli,gpoli)
 501   ),
 502   hli:_COMA_ppo(gpoli,hli,[]),
 503  if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 504 )$
 505 _COMA_rlpr(poles,fx,fy) :=   /* root locus plot range */
 506 block(                       /* fx,fy ... oversizing factors */
 507   [immax,remin,remax,xmin,xmax,ymax,xm,aspect_ratio:0.6],
 508   immax:apply(max,imagpart(poles)),
 509   remax:apply(max,realpart(poles)),
 510   remin:apply(min,realpart(poles)),
 511   ymax:fy*immax,
 512   xm:(remax+remin)/2,
 513   xmax:xm+fx*(remax-remin)/2,
 514   xmin:xm-fx*(remax-remin)/2,
 515   if is(xmax=xmin) then block(xmax:xmax+1,xmin:xmin-1),
 516   ymax:max(ymax,(xmax-xmin)/2*aspect_ratio),
 517   xmax:max(xmax,xm+ymax/aspect_ratio),
 518   xmin:min(xmin,xm-ymax/aspect_ratio),
 519   ['xrange=[xmin,xmax],'yrange=[-ymax,ymax]]
 520 )$
 521 _COMA_rlsp(z1,z2) :=  /* root locus sort points */
 522 block([l:z1], for i:1 step 1 thru length(z1-1) do block(
 523   for j:i+1 step 1 thru length(z1) do
 524   if cabs(z2[i]-l[j])<cabs(z2[i]-l[i]) then [l[j],l[i]]:[l[i],l[j]]),l)$
 525
 526 /* ======================================================================= */
 527 /* Stability and related functions                                         */
 528 /* ======================================================================= */
 529
 530 /* ----------------------------------------------------------------------- */
 531 /* poles(f) - poles of the transfer function f(s)                          */
 532 /* ----------------------------------------------------------------------- */
 533 poles(f):=
 534 block([res,polyfactor:false,ratprint:false,dt],
 535   res:map(lambda([ff],
 536     if ntfp(ff) then map(rhs,chop(allroots(%i*expand(float(denom(ff))))))
 537        else []
 538   ), xthru(first(_COMA_strip(float(f))))),
 539   if listp(f) then res else res[1]
 540 )$
 541
 542 /* ----------------------------------------------------------------------- */
 543 /* zeros(f) - zeros of the transfer function f(s)                          */
 544 /* ----------------------------------------------------------------------- */
 545 zeros(f):=
 546 block([res,polyfactor:false,ratprint:false,dt],
 547   res:map(lambda([ff],
 548     if ntfp(ff) then map(rhs,chop(allroots(%i*expand(float(num(ff))))))
 549        else []),xthru(first(_COMA_strip(float(f))))),
 550   if listp(f) then res else res[1]
 551 )$
 552
 553 /* ----------------------------------------------------------------------- */
 554 /* hurwitz(p) - Hurwitz-determinants of the polynomial p(s)                */
 555 /* ----------------------------------------------------------------------- */
 556 hurwitz(p):=
 557 block(
 558   [cli,n,i,k,hli,mat,res],
 559   if not listp(p) then p:[p],
 560   p:map(expand,p),
 561   res:map(lambda([pp],
 562   block(
 563     hli:[],
 564     n:hipow(pp,'s),
 565     cli:reverse(coefficient_list(pp,s)),
 566     cli:append(makelist(0,i,1,n-1),cli,makelist(0,i,1,n  )),
 567     hli:makelist(mat:apply(matrix,makelist(
 568       /*  makelist('cli[''((2*(n-j)+1)+(i-1))],j,1,k),i,1,k)) ,k,2,n), HAA, 14.2.2018 */
 569         makelist('cli[((2*(n-j)+1)+(i-1))],j,1,k),i,1,k)) ,k,2,n),
 570     hli:ev(hli,nouns), /* trick to prevent substituting variables k and n */
 571     hli:map(determinant,hli),
 572     if n>1 then hli else [])
 573   ), p),
 574   if length(res)=1 then res[1] else res
 575 )$
 576
 577 /* ----------------------------------------------------------------------- */
 578 /* stable_area(f,p1,p2) - plots the stable area with respect to p1 and p2  */
 579 /* ----------------------------------------------------------------------- */
 580 stable_area_outlined(f,p1,p2,[opts]):=
 581 block(
 582   [wxplot_size:assoc(dimensions,coma_defaults),
 583     x,x1,x2,y,y1,y2,defs,hli,fli,h,wxx,ratprint:false],
 584   opts:flatten(opts),
 585   defs:[xrange=[0,1],yrange=[0,1]],
 586   f:flatten([f]),
 587   h:map(lambda([u],if _COMA_go(u) then u else append(hurwitz(denom(u)),
 588               coefficient_list(denom(u)))),f),
 589   fli:map(lambda([u],if _COMA_go(u) then u else apply('min,u)),h),
 590   fli:map(lambda([u],if _COMA_go(u) then u else implicit(u,p1,x1,x2,p2,y1,y2)),fli),
 591   hli:_COMA_ppo(fli,defs,opts),
 592   [x1,x2]:get_option('xrange,hli),
 593   [y1,y2]:get_option('yrange,hli),
 594   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 595 )$
 596
 597 stable_area(f,p1,p2,[opts]):=   /* depiction as region 24.2.2019 */
 598 block(
 599   [wxplot_size:assoc(dimensions,coma_defaults),
 600     x,x1,x2,y,y1,y2,defs,hli,fli,h,wxx,ratprint:false],
 601   opts:flatten(opts),
 602   defs:[xrange=[0,1],yrange=[0,1],x_voxel=50,y_voxel=50,
 603        user_preamble="set style fill transparent solid 0.2 noborder"],
 604   f:flatten([f]),
 605   h:map(lambda([u],if _COMA_go(u) then u else
 606       map(lambda([v],v>0),append(hurwitz(denom(u)),coefficient_list(denom(u))))),f),
 607   fli:map(lambda([u],if _COMA_go(u) then u else apply("and",u)),h),
 608   fli:map(lambda([u],if _COMA_go(u) then u else region(u,p1,x1,x2,p2,y1,y2)),fli),
 609   hli:_COMA_ppo(fli,defs,opts),
 610   [x1,x2]:get_option('xrange,hli),
 611   [y1,y2]:get_option('yrange,hli),
 612   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 613 )$
 614
 615 unstable_area(f,p1,p2,[opts]):=   /* depiction as region 24.2.2019 */
 616 block(
 617   [wxplot_size:assoc(dimensions,coma_defaults),
 618     x,x1,x2,y,y1,y2,defs,hli,fli,h,wxx,ratprint:false],
 619   opts:flatten(opts),
 620   defs:[xrange=[0,1],yrange=[0,1],x_voxel=50,y_voxel=50,
 621        user_preamble="set style fill transparent solid 0.2 noborder"],
 622   f:flatten([f]),
 623   h:map(lambda([u],if _COMA_go(u) then u else
 624       map(lambda([v],v<0),append(hurwitz(denom(u)),coefficient_list(denom(u))))),f),
 625   fli:map(lambda([u],if _COMA_go(u) then u else apply("or",u)),h),
 626   fli:map(lambda([u],if _COMA_go(u) then u else region(u,p1,x1,x2,p2,y1,y2)),fli),
 627   hli:_COMA_ppo(fli,defs,opts),
 628   [x1,x2]:get_option('xrange,hli),
 629   [y1,y2]:get_option('yrange,hli),
 630   if wxx#false then apply(wxdraw2d,ev(hli)) else apply(draw2d,ev(hli))
 631 )$
 632
 633 /* ----------------------------------------------------------------------- */
 634 /* stablep(f) - checks whether the transfer function f(s) is stable        */
 635 /* ----------------------------------------------------------------------- */
 636 stablep(f) :=
 637 block([ratprint:false],
 638   if listp(f) then
 639     map(lambda([ff],
 640       not (member(true,map(lambda([u],is(realpart(u)>=0)),poles(ff))))
 641     ), f)
 642   else not (member(true,map(lambda([u],is(realpart(u)>=0)),poles(f))))
 643 )$
 644
 645 /* ----------------------------------------------------------------------- */
 646 /* phase_crossover(f) - calculates the phase crossover frequencies of f    */
 647 /* ----------------------------------------------------------------------- */
 648 phase_crossover(f):=
 649 block([re,im,sol,res,ratprint:false,omega,dt],
 650   [f,dt]:_COMA_strip(f),
 651   res:map(lambda([ff,dd],
 652   if dd=1 then block(     /* without time delay */
 653     re:realpart(ev(ff,s=%i*omega)),
 654     im:imagpart(ev(ff,s=%i*omega)),
 655     sol:algsys([num(im)=0],[omega]),
 656     sol:sublist(flatten(sol),lambda([u],subst(u,denom(re))#0 and
 657          is(ev(re,u)<0) and is(rhs(u)>0))),
 658     sort(sol,lambda([u,v],rhs(u)<rhs(v))))
 659   else block([],    /*  with time delay */
 660     ['omega=find_root(phase(ff*dd)+180,omega,0,1e10)])
 661   ), f,dt),
 662   if length(res)=1 then res[1] else res
 663 )$
 664
 665 /* ----------------------------------------------------------------------- */
 666 /* gain_crossover(f) - calculates the gain crossover frequencies of f      */
 667 /* ----------------------------------------------------------------------- */
 668 gain_crossover(f):=
 669 block(
 670   [roots,sol,res,polyfactor:false,ratprint:false,zz,omega,dt],
 671   [f,dt]:_COMA_strip(f),
 672   res:map(lambda([ff,dd],
 673     if ntfp(ff)
 674     then block
 675     (
 676       ff:xthru(subst(s=%i*omega,ff)),
 677       zz:num(xthru((cabs(num(ff))/cabs(denom(ff)))**2-1)),
 678       roots:if hipow(zz,omega)>0 then realroots(zz) else [],
 679       sol:float(sublist(roots,lambda([u],is(part(u,2)>0)))),
 680       if listp(sol) then sort(sol,lambda([u,v],rhs(u)<rhs(v)))
 681     ) else []
 682   ), f,dt),
 683   if length(res)=1 then res[1] else res
 684 )$
 685
 686 /* ----------------------------------------------------------------------- */
 687 /* phase_margin(f) - calculates the phase margin of f                      */
 688 /* ----------------------------------------------------------------------- */
 689 phase_margin(f):=
 690 block(
 691   [res,h,ratprint:false],
 692   if not listp(f) then f:[f],
 693   res:map(lambda([ff],
 694     h:gain_crossover(ff),
 695     180+map(lambda([v],float(ev(phase(ff),v))),h)),f),
 696   if length(res)=1 then res[1] else res
 697 )$
 698
 699 /* ----------------------------------------------------------------------- */
 700 /* gain_margin(f) - calculates the gain margin of f                        */
 701 /* ----------------------------------------------------------------------- */
 702  gain_margin(f):=
 703  block([res,h,ratprint:false,omega],
 704    if not listp(f) then f:[f],
 705    res:map(lambda([ff],
 706        h:float(phase_crossover(ff)),
 707        (-1)/realpart(map(lambda([v],ev(ff,s=%i*omega,v,eval)),h))),f),
 708    if length(res)=1 then res[1] else res
 709 )$
 710
 711 /* ----------------------------------------------------------------------- */
 712 /* stability_limit(f,par) - calculates conditions for imaginary poles      */
 713 /* ----------------------------------------------------------------------- */
 714 stability_limit(f,par,[opt]):=
 715 block(
 716   [eqs,nn,sol,res,assume_pos:true,omega,ratprint:false,dt],
 717   [f,dt]:_COMA_strip(f),
 718   if setify(dt)#{1} then disp(_COMA_msg[3]),
 719   res:map(lambda([ff],
 720   block(
 721     nn:ev(denom(ff),s=%i*omega),
 722     eqs:map(lambda([u],u(nn)=0),[realpart,imagpart]),
 723     sol:ev(solve(eqs,[par,omega]),realonly=true),
 724     sol:sublist(sol,lambda([u], assoc(omega,u)>0)))
 725   ), f),
 726   if length(res)=1 then res[1] else res
 727 )$
 728
 729 /* ----------------------------------------------------------------------- */
 730 /* damping(f) - negative realpart of the rightmost pole of f(s)            */
 731 /* ----------------------------------------------------------------------- */
 732 damping(f) :=
 733 block(
 734   [polyfactor:true,ratprint:false,expr,p,res],
 735   res:map(lambda([ff],
 736   block(
 737     p:poles(float(xthru(ff))),
 738     expr:map(realpart,p),
 739     -apply(max,expr))
 740   ), first(_COMA_strip(f))),
 741   if length(res)=1 then res[1] else res
 742 )$
 743
 744 /* ----------------------------------------------------------------------- */
 745 /* damping_ratio(f) - minimal value of all damping ratios                  */
 746 /* ----------------------------------------------------------------------- */
 747 damping_ratio(f) :=
 748 block(
 749   [polyfactor:true,ratprint:false,expr,p,res],
 750   res:map(lambda([ff],
 751   block(
 752     p:poles(float(xthru(ff))),
 753     expr:map(lambda([u],realpart(u)/sqrt(realpart(u)**2+imagpart(u)**2)),p),
 754     -apply(max,expr))
 755   ), first(_COMA_strip(f))),
 756   if length(res)=1 then res[1] else res
 757 )$
 758
 759 /* ======================================================================= */
 760 /* Transfer Function Related                                               */
 761 /* ======================================================================= */
 762
 763 /* ----------------------------------------------------------------------- */
 764 /* normalize(f,a) - applies the similarity theorm on a transfer functi     */
 765 /* ----------------------------------------------------------------------- */
 766 normalize(f,[ _a]):=
 767 block([alpha:1,k,cl,res,assume_pos:true],
 768   if length(_a)>0 then alpha: _a[1],
 769   if not listp(f) then f:[f],
 770   res:map(lambda([u],
 771     u:sum_form(u),
 772     cl:coefficient_list(denom(u),s),
 773     k:last(cl)**(1/(length(cl)-1)),
 774     subst(s=s/alpha/k,u)),f),
 775   res:fullmap(lambda([u],if numberp(u) and abs('round(u)-u) < 1.0e-10
 776     and round(u)>0 then 'round(u) else u),sum_form(res)),
 777   if length(res)=1 then res[1] else res)$
 778
 779 /* ----------------------------------------------------------------------- */
 780 /* tf(type,a,..) - generates special transferfunctions                     */
 781 /* ----------------------------------------------------------------------- */
 782
 783 tf(f,[p]):=
 784 block([k,T1,T2,D,omega_n,Td,Tn,Tv,Ti,kr,Ta,Tb,nz:1,nn:2,h,pz,pn,c,d,eps,a,b],
 785   local(c,d),
 786   if f='pt1 then
 787      (if length(p)>0 then k:p[1],if length(p)>1 then T1:p[2],return(k/(1+s*T1))),
 788   if f='pt2  and length(p)>0 and p[1]='prod then
 789      (if length(p)>1 then k:p[2],if length(p)>2 then T1:p[3], if length(p)>3 then T2:p[4],
 790      return(k/((1+s*T1)*(1+s*T2)))),
 791   if f='pt2 then
 792      (if length(p)>0 then k:p[1],if length(p)>1 then D:p[2],if length(p)>2 then omega_n:p[3],
 793      return(k/(1+2*D/omega_n*s+s**2/omega_n**2))),
 794   if f='dt1 then
 795      (if length(p)>0 then Td:p[1],if length(p)>1 then T1:p[2],
 796      return(s*Td/(1+s*T1))),
 797   if f='it1 then
 798      (if length(p)>0 then Ti:p[1],if length(p)>1 then T1:p[2],
 799      return(1/(s*Ti*(1+s*T1)))),
 800   if f='pdt1 then
 801      (if length(p)>0 then k:p[1],if length(p)>1 then T1:p[2],if length(p)>2 then T2:p[3],
 802      return(k*(1+s*T1)/(1+s*T2))),
 803   if f='pd then
 804      (if length(p)>0 then k:p[1],if length(p)>1 then T1:p[2],
 805      return(k*(1+s*T1))),
 806   if f='allpass then
 807      (if length(p)>0 then k:p[1],if length(p)>1 then T1:p[2],
 808      return(k*(1-s*T1)/(1+s*T1))),
 809   if f='pi then
 810      (if length(p)>0 then kr:p[1],if length(p)>1 then Tn:p[2],
 811      return(kr*(1+1/(s*Tn)))),
 812   if f='pid  and length(p)>0 and p[1]='prod then
 813      (if length(p)>1 then Ta:p[2],if length(p)>2 then Tb:p[3], if length(p)>3 then Ti:p[4],
 814      return((1+s*Ta)*(1+s*Tb)/(s*Ti))),
 815   if f='pid then
 816      (if length(p)>0 then kr:p[1],if length(p)>1 then Tn:p[2],if length(p)>2 then Tv:p[3],
 817      return(kr*(1+1/(s*Tn)+s*Tv))),
 818   if f='pidt1  and length(p)>0 and p[1]='prod then
 819      (if length(p)>1 then Ta:p[2],if length(p)>2 then Tb:p[3], if length(p)>3 then Ti:p[4],
 820      if length(p)>4 then T1:p[5],return((1+s*Ta)*(1+s*Tb)/(s*Ti*(1+s*T1)))),
 821   if f='pidt1 then
 822      (if length(p)>0 then kr:p[1],if length(p)>1 then Tn:p[2],if length(p)>2 then Tv:p[3],
 823      if length(p)>3 then T1:p[4], return(kr*(1+1/(s*Tn)+s*Tv/(1+s*T1)))),
 824   if f='generic then
 825      (if length(p)>0 then nz:p[1],if length(p)>1 then nn:p[2],if length(p)>2 then c:p[3],
 826      if length(p)>3 then d:p[4],return(sum(c[i]*s**i,i,0,nz)/sum(d[i]*s**i,i,0,nn))),
 827   if f='random then
 828   (
 829     nn:3,h:20,
 830     if length(p)>1 and map(numberp,rest(p,length(p)-2))=[true,true]
 831     then (h:last(p),nn:second(reverse(p)))
 832     else if length(p)>0 and numberp(last(p)) then (nn:last(p)),
 833     nz:random(nn),
 834     if member('normalized,p) then
 835     (
 836       pz:1, pn:1+s**nn,
 837       for i:1 thru nn-1 do pn: pn + (1+random(h)) * s**i,
 838       for i:1 thru nz do pz: pz + (1+random(h)) * s**i
 839     )
 840     else
 841     (
 842       pz:0, pn:0,
 843       for i:0 thru nn do pn: pn + (1+random(h)) * s**i,
 844       for i:0 thru nz do pz: pz + (1+random(h)) * s**i
 845     ),
 846     f:pz/pn,
 847     while not stablep(f) and member('stable,p) do
 848     (
 849       if member('normalized,p) then (pn:1+s**nn,for i:1 thru nn-1 do pn: pn+(1+random(h))*s**i)
 850       else (pn:0,for i:0 thru nn do pn: pn + (1+random(h))*s**i),
 851       f:pz/pn
 852     ),return(f)
 853   ),
 854   if f='butterworth then
 855   (
 856     nn:if length(p)>0 then p[1] else 3,
 857     h:if length(p)>1 then p[2] else 1,
 858     return(float(h**nn/([1,s+h,1][1+mod(nn,2)]*
 859        product(expand((s**2-2*s*h*cos((2*k+nn-1)/(2*nn)*%pi)+h**2)),k,1,nn/2))))
 860   ),
 861   if f='pseudobutt then
 862   (
 863     nn:if length(p)>0 then p[1] else 3, h:if length(p)>1 then p[2] else 1,
 864     return(float(h**nn/([1,s+h,1][1+mod(nn,2)]*
 865        product(expand((s+h*(sin(k*%pi/(nn+1))+%i*cos(k*%pi/(nn+1))))*
 866        (s+h*(sin(k*%pi/(nn+1))-%i*cos(k*%pi/(nn+1))))),k,1,nn/2))))
 867   ),
 868 if f='ptn then /* similar to butterworth, poles confined to a given angle */
 869   (
 870     nn:if length(p)>0 then p[1] else 3,
 871     h:if length(p)>1 then p[2] else 1,
 872     eps:if length(p)>2 then p[3] else 0,
 873     if evenp(nn) then return(float(h**nn/(product(expand(
 874             (s**2+2*s*h*cos((k-1/2)*%pi/(nn-1)*eps/90)+h**2)),k,1,nn/2))))
 875     else return(float(h**nn/((h+s)*product(expand(
 876             (s**2+2*s*h*cos(k/(nn-1)*%pi*eps/90)+h**2)),k,1,nn/2))))
 877   ),
 878 if f='chebychev or f='tschebytscheff then
 879   (
 880     nn:if length(p)>0 then p[1] else 3,
 881     h:if length(p)>1 then p[2] else 1,
 882     eps:if length(p)>2 then p[3] else 0.2,
 883     return(subst(s=s/h,
 884       if evenp(nn)
 885       then 1/float(product(
 886         ( b:1/(cosh(eps)**2-cos((2*k-1)*%pi/(2*nn))**2),
 887           a:2*b*sinh(eps)*cos((2*k-1)*%pi/(2*nn)),
 888           1+a*s+b*s**2),k,1,nn/2))
 889       else 1/float((1+s/sinh(eps)) * product(
 890         ( b:1/(cosh(eps)**2-cos((k-1)*%pi/nn)**2),
 891           a:2*b*sinh(eps)*cos((k-1)*%pi/nn),
 892           1+a*s+b*s**2),k,2,(nn+1)/2))))
 893   ),
 894 if f='bessel then
 895   (
 896     nn:if length(p)>0 then p[1] else 3,
 897     h:if length(p)>1 then p[2] else 1,
 898     return(normalize(1/(1+sum(product(2*(nn-i+1)/i/(2*nn-i+1),i,1,k)*(s)**k,k,1,nn)),h))
 899   ),
 900 if f='ise then
 901   (
 902     nn:if length(p)>0 then p[1] else 3,
 903     h:if length(p)>1 then p[2] else 1,
 904     return(normalize(1/sum(binomial(nn-floor((k+1)/2),floor(k/2))*s^(nn-k),k,0,nn),h))
 905   )
 906 )$
 907
 908 /* ----------------------------------------------------------------------- */
 909 /* tranftype(f) - returns the type of the transfer function f as a string  */
 910 /* ----------------------------------------------------------------------- */
 911 tranftype(f) :=
 912 block([res],
 913   res:map(lambda([ff],
 914   block(
 915     [zz,nn,z0:1,z1,n0:1,n1,ratprint:false],
 916     zz:coefficient_list(num(xthru(ff)),s),
 917     nn:coefficient_list(denom(xthru(ff)),s),
 918     z1:length(zz),
 919     n1:length(nn),
 920     while is(zz[z0]=0) do z0:z0+1,
 921     while is(nn[n0]=0) do n0:n0+1,
 922     if is([z1,n1]=[1,1]) then return("P"),
 923     if is([z1,n0]=[1,1]) then return(concat("PT",n1-1)),
 924     if is([z1,n1]=[1,2]) then return(concat("I")),
 925     if is([z1,n0]=[1,n1]) then return(concat("I",n1-1)),
 926     if is([z1,n0]=[1,2]) then return(concat("IT",n1-n0)),
 927     if is([z1]=[1]) then return(concat("I",n0-1,"T",n1-n0)),
 928     if is([z0,z1,n1]=[2,2,1]) then return("D"),
 929     if is([z0,n1]=[z1,1]) then return(concat("D",z0-1)),
 930     if is([z0,z1,n0]=[2,2,1]) then return(concat("DT",n1-1)),
 931     if is([z1,n1]=[2,1]) then return("PD"),
 932     if is([z0,n1]=[1,1]) then return(concat("PD",z1-1)),
 933     if is([z0,z1,n0]=[1,2,1]) then return(concat("PDT",n1-1)),
 934     if is([n1]=[1]) then return(concat("D",n0-1,"T",n1-n0)),
 935     if is([z1,n0,n1]=[2,2,2]) then return("PI"),
 936     if is([z1,n0,n1]=[3,2,2]) then return("PID"),
 937     if is([z1,n0]=[3,2]) then return(concat("PIDT",n1-n0)),
 938     if is([z0,z1,n1]=[1,2,1]) then return("PD"),
 939     if is([z0,n0]=[1,1]) then return(concat("PD",z1-z0,"T",n1-n0)),
 940     return("any"))
 941   ), first(_COMA_strip(f))),
 942   if length(res)=1 then res[1] else res
 943 )$
 944 tftype(f):=tranftype(f)$
 945
 946 /* ----------------------------------------------------------------------- */
 947 /* rantranf(n) - Random transfer function of order n                       */
 948 /* ----------------------------------------------------------------------- */
 949 rantranf(n) :=
 950 block(
 951   [zz:0, nn:0, kz:random(n), kn:n, i, f],
 952   for i:0 step 1 thru kn do nn: nn + (1+random(10)) * s**i,
 953   for i:0 step 1 thru kz do zz: zz + (1+random(10)) * s**i,
 954   zz/nn
 955 )$
 956
 957 /* ----------------------------------------------------------------------- */
 958 /* stable_rantranf(n) - stable random transfer function of order n         */
 959 /* ----------------------------------------------------------------------- */
 960
 961 stable_rantranf(n) :=
 962 block(
 963   [ntest:1,zz:0, nn:0, kz, kn, i, f],
 964   kn:min(n,7),
 965   kz:random(kn),
 966   for i:0 step 1 thru kn do nn: nn + (1+random(10)) * s**i,
 967   for i:0 step 1 thru kz do zz: zz + (1+random(10)) * s**i,
 968   f:(zz/nn),
 969   unless stablep(f) do block(
 970     ntest:ntest+1,
 971     nn:0,
 972     for i:0 step 1 thru kn do nn: nn + (1+random(10)) * s**i,
 973     f:(zz/nn)),
 974   return(f)
 975 )$
 976
 977 /* ----------------------------------------------------------------------- */
 978 /* gentranf(c,nz,d,nn) - general transfer function of order nn             */
 979 /* ----------------------------------------------------------------------- */
 980 gentranf(c,nz,d,nn) := sum(c[i]*s**i,i,0,nz)/sum(d[i]*s**i,i,0,nn)$
 981
 982 /* ----------------------------------------------------------------------- */
 983 /* impedance_chain(z1,z2,...n) - transfer function of an impedance chain   */
 984 /* ----------------------------------------------------------------------- */
 985 impedance_chain(z1,z2,[zi]) :=
 986 block(
 987   [_z1:[],_num:1,_sol,_vars,_n,_k,_ue,_i,ratprint:false,_eqs,
 988              _z:flatten(append([z1,z2],zi))],
 989 /* Repetition of the chain at odd number of parameters */
 990   if oddp(length(_z)) then block(
 991      _num:last(_z),
 992      _z:rest(_z,-1),
 993      for _j:1 thru _num do _z1:append(_z,_z1),
 994      _z:_z1),
 995   _n:length(_z)/2,
 996 /* Building and solving the mesh-equations */
 997   _eqs:makelist(_z[2*(_k-1)]*_i[_k-1]=_z[2*_k]*_i[_k]+_z[2*_k-1]*sum(_i[_j],_j,_k,_n),_k,2,_n),
 998   _eqs:cons(_ue=_z[1]*sum(_i[_j],_j,1,_n)+_z[2]*_i[1],_eqs),
 999   _vars:makelist(_i[_k],_k,1,_n),
1000   _sol:linsolve(_eqs,_vars), /* linsolve works better than solve (why?) */
1001 /* Calculation of the transfer function as the ratio Ua/Ue */
1002   ratsimp(subst(_sol,_z[2*_n]*_i[_n]/ _ue))
1003 )$
1004
1005 /* ----------------------------------------------------------------------- */
1006 /* nilt(f,s,t) - inverse Laplace transform of f with numerically           */
1007 /*               calculated poles                                          */
1008 /* ----------------------------------------------------------------------- */
1009 oilt(f,[st]):=  /* Old version - just for test purposes */
1010 block([polyfactor:true,ratprint:false,ft,res,dt,fl,lvar,tvar,dummy],
1011   if length(st)=2 then [lvar,tvar]:st else [lvar,tvar]:['s,'t],
1012   if not listp(f) then f:[f],
1013   res:map(lambda([ff],
1014     [fl,dt]: if op(ff)="+" then _COMA_strip(args(ff)) else _COMA_strip(ff),
1015     ft:map(lambda([v],(if unit_step_included then unit_step(t) else 1)*
1016         if try_symbolic_ilt
1017            then block([r:ilt(v,lvar,tvar)],
1018                 if member("ilt",map(string,operators(r))) /* cannot be transformed analytically */
1019                    then ev(expand(ilt(float(num(v)/allroots(float(denom(v)))),lvar,tvar)),numer)
1020                    else r)
1021            else ev(expand(ilt(float(num(v)/allroots(float(denom(v)))),lvar,tvar)),numer)
1022                        ),xthru(fl)),
1023    ft:map(lambda([u,v],ev(u,t=t+v)),ft,log(dt)/s),
1024  /*  apply("+",ft)),expandwrt(f,s)),      15.2.2019 */
1025    apply("+",ft)),f),
1026     if length(res)=1 then res[1] else res
1027 )$
1028 nilt(f,[st]):=  /* New version - 18.2019 */
1029 block([polyfactor:true,ratprint:false,ft,res,dt,fl,lvar,tvar,dummy],
1030   if length(st)=2 then [lvar,tvar]:st else [lvar,tvar]:['s,'t],
1031   if not listp(f) then f:[f],
1032   res:map(lambda([ff],
1033     [fl,dt]: _COMA_sdt(ff),
1034     ft:map(lambda([v],(if unit_step_included then unit_step(t) else 1)*
1035         if try_symbolic_ilt
1036            then block([r:ilt(v,lvar,tvar)],
1037                 if member("ilt",map(string,operators(r))) /* cannot be transformed analytically */
1038                    then ev(expand(ilt(float(num(v)/allroots(float(denom(v)))),lvar,tvar)),numer)
1039                    else r)
1040            else ev(expand(ilt(float(num(v)/allroots(float(denom(v)))),lvar,tvar)),numer)
1041                        ),xthru(fl)),
1042    ft:map(lambda([u,v],ev(u,t=t+v)),ft,log(dt)/s),
1043  /*  apply("+",ft)),expandwrt(f,s)),      15.2.2019 */
1044    apply("+",ft)),f),
1045     if length(res)=1 then res[1] else res
1046 )$
1047 try_symbolic_ilt:true$
1048 unit_step_included:true$
1049
1050 /* ----------------------------------------------------------------------- */
1051 /* closed_loop(Fo) - calculates the closed loop transfer function Fw       */
1052 /* ----------------------------------------------------------------------- */
1053 closed_loop(f,[x]):=
1054 block(
1055   [res,ratprint:false,hf,ht,n:5],
1056   if length(x)>0 then n:first(x),
1057   if not listp(f) then f:[f],
1058   res:map(lambda([ff],
1059     [hf,ht]:_COMA_strip(ff),
1060     if ht=[1] then (if not mapatom(ff/(1+ff)) then map(expand,xthru(ff/(1+ff))) else ff/(1+ff))
1061               else sum((-1)**((_k)+1)*first(hf)**(_k)*first(ht)**(_k),(_k),1,n)
1062    ), f),
1063   if length(res)=1 then res[1] else res
1064 )$
1065
1066 /* ----------------------------------------------------------------------- */
1067 /* open_loop(Fw) - calculates the open loop transfer function Fo           */
1068 /* ----------------------------------------------------------------------- */
1069 open_loop(f):=
1070 block(
1071   [res,ratprint:false],
1072   if not listp(f) then f:[f],
1073   res:map(lambda([ff],map(expand,xthru(ff/(1-ff)))), f),
1074   if length(res)=1 then res[1] else res
1075 )$
1076
1077 /* ----------------------------------------------------------------------- */
1078 /* time_delay(tt,n,k) - Pade-approximation of order n (numerator order k)  */
1079 /* ----------------------------------------------------------------------- */
1080 time_delay(tt,n,[k]):=
1081 if n=0 then 1
1082 else block([ratprint:false,numer:false,keepfloat:false,h],
1083     h:if emptyp(k) then n-1 else first(k),
1084     lambda([u],num(u)/denom(u))(first(pade(taylor(exp(-s*tt),s,0,h+n),h,n)))
1085 )$
1086
1087 /* ----------------------------------------------------------------------- */
1088 /* ntfp(f) - checks whether f has only numeric coefficients             */
1089 /* ----------------------------------------------------------------------- */
1090 ntfp(f) :=
1091 block([ratprint:false,f1,zz,nn,res],
1092   if not listp(f) then f:[f],
1093   res:map(lambda([ff],
1094      f1:xthru(ff),
1095      zz:coefficient_list(num(f1),s),
1096      nn:coefficient_list(denom(f1),s),
1097      if is(setify(map(numberp,zz))={true}) and
1098        is(setify(map(numberp,nn))={true})
1099        then true else false
1100   ), float(f)),
1101   if length(res)=1 then res[1] else res
1102 )$
1103
1104 /* ----------------------------------------------------------------------- */
1105 /* transfer_function - transfer function from state space or equations     */
1106 /* ----------------------------------------------------------------------- */
1107 transfer_function(a,[d]) :=
1108 block(
1109 [ratprint:false,A,B,C,D,sp:systemp(a)],
1110 if sp then block(
1111    if length(a)=4 then [A,B,C,D]:a,
1112    if length(a)=3 then
1113       [A,B,C,D]:endcons(zeromatrix(length(a[3]),length(a[2][1])),a)),
1114 if sp and matrixp(D) and length(D)=1 and length(D[1])=1 then D:D[1][1],
1115 if sp then return(map(expand,xthru(C.invert(s*ident(length(A))-A).B+D))),
1116 if matrixp(a) then block(
1117    if length(d)=3 then [A,B,C,D]:[a,d[1],d[2],d[3]]
1118    else [A,B,C,D]:[a,d[1],d[2],zeromatrix(length(d[2]),length(d[1][1]))],
1119    if matrixp(D) and length(D)=1 and length(D[1])=1 then D:D[1][1],
1120    map(expand,ratsimp(C.invert(s*ident(length(A))-A).B+D)))
1121 else block([m],
1122 m:coefmatrix(subst(linsolve(a,d[1]),flatten([d[3]])),flatten([d[2]])),
1123 if not listp(d[2]) and not listp(d[3]) then m[1][1] else m
1124 ))$
1125
1126 /* ----------------------------------------------------------------------- */
1127 /* sum_form - One of the 4 standard forms of a transfer function           */
1128 /*    (making one of the leading or absolute coefficients to 1)            */
1129 /* ----------------------------------------------------------------------- */
1130 sum_form(f,[nli]):=
1131 block([_zz,_nn,_zli,_nli,_sz,_kz,_kn,_n,_zli1,_nli1,_res,dt,_h],
1132   _n: if emptyp(nli) then 4 else nli[1],
1133   [f,dt]:_COMA_strip(f),
1134   _res:map(lambda([ff],
1135     _sz:if member(s,listofvars(ff)) then s
1136         elseif member(z,listofvars(ff)) then z else s,
1137     [_zz,_nn]:lambda([u],[num(u),denom(u)])(xthru(ff)),
1138     _zli:coefficient_list(_zz,_sz),
1139     _nli:coefficient_list(_nn,_sz),
1140     _zli1:delete(0,_zli),
1141     _nli1:delete(0,_nli),
1142     _kz:[last(_zli1),first(_zli1),last(_nli1),first(_nli1),last(_zli1),first(_zli1)][_n],
1143     _kn:[last(_zli1),first(_zli1),last(_nli1),first(_nli1),last(_nli1),first(_nli1)][_n],
1144     [_zli,_nli]:[_zli,_nli]/ [_kz,_kn],
1145     fullermap(lambda([u],if numberp(u) and equal(fix(u),u) then fix(u) else u),
1146       [1,1,1,1,_kz/ _kn,_kz/ _kn][_n]*apply("+",float(_zli)*makelist(_sz**i,i,0,length(_zli)-1))/
1147       apply("+",float(_nli)*makelist(_sz**i,i,0,length(_nli)-1)))),f)*dt,
1148     if length(_res)=1 then _res[1] else _res
1149 )$
1150
1151 /* ----------------------------------------------------------------------- */
1152 /* product_form - Factorized forms of a transfer function (2017-06-08)     */
1153 /* ----------------------------------------------------------------------- */
1154 product_form(f,[n]):=
1155 block([res,dt],
1156   [f,dt]:_COMA_strip(f),
1157   res:map(lambda([ff],
1158   if ntfp(ff) and not mapatom(ff)
1159     then block([inflag:true,polyfactor:true,ratprint:false,fff,zz,nn,zli,nli,zdiv,ndiv],
1160     fff:xthru(float(ff)),
1161     fff:chop(allroots(num(fff))/allroots(denom(fff))),
1162     zz:num(fff),
1163     nn:denom(fff),
1164     if not emptyp(n) and n[1]=1 then return(fullmap(lambda([u],if numberp(u)
1165        and abs(round(u)-u)<1e-10 then round(u) else u),(zz/nn))), /* high coefficient is 1 */
1166     zli:if not mapatom(zz) and op(zz)="*" then args(zz) else [zz],
1167     nli:if not mapatom(nn) and op(nn)="*" then args(nn) else [nn],
1168     zli:map(lambda([u],coefficient_list(u,s)),zli),
1169     nli:map(lambda([u],coefficient_list(u,s)),nli),
1170     zdiv:map(lambda([u],first(delete(0,u))),zli),
1171     ndiv:map(lambda([u],first(delete(0,u))),nli),
1172     fullmap(lambda([u],if numberp(u) and abs('round(u)-u) < 1.0e-10 and round(u)>0 then 'round(u) else u),
1173       ((apply("*",zdiv)*
1174         apply("*",map(lambda([u],u.makelist(s^k,k,0,length(u)-1)),zli/zdiv)))/
1175         apply("*",ndiv))/
1176         apply("*",map(lambda([u],u.makelist(s^k,k,0,length(u)-1)),nli/ndiv)))
1177         ) else ff),f)*dt,
1178              if length(res)=1 then res[1] else res)$
1179
1180 /* ----------------------------------------------------------------------- */
1181 /* npartfrac - Partial fraction expansion with numeric coefficients        */
1182 /* ----------------------------------------------------------------------- */
1183 npartfrac(G):=
1184 block([polyfactor:true,ratprint:false],
1185    map(lambda([u],sum_form(u,3)),partfrac(num(G)/allroots(denom(G)),
1186       first(listofvars(G))))
1187 )$
1188
1189 /* ======================================================================= */
1190 /* Optimization, Controller Design                                         */
1191 /* ======================================================================= */
1192
1193 /* ----------------------------------------------------------------------- */
1194 /* ise(f) - integral of squared error                                      */
1195 /* ----------------------------------------------------------------------- */
1196 ise(f):=
1197 block(
1198   [n,sol,res,ratprint:false,dn,varlist,dt],
1199   [f,dt]:_COMA_strip(f),
1200   if setify(dt)#{1} then disp(_COMA_msg[3]),
1201   res:map(lambda([ff],
1202   block(
1203     block(                  /* Trick: solve darf gleichnamige Symbole in f */
1204       [a,b,aa,bb,c,d,i],    /* nicht evaluieren! (macht es aber sonst :-)  */
1205       [c,d]:[num(ff),denom(ff)],
1206       n:hipow(d,s),
1207       dn:coeff(expand(d),s,n),
1208       aa:sum('a[i]*s**i,i,0,n-1),
1209       bb:sum('b[i]*s**i,i,0,n-1),
1210       eq:coefficient_list(expand(c*subst(s=-s,c)-aa*subst(s=-s,d)-bb*d),s),
1211       varlist:flatten([makelist('a[i],i,0,n-1),makelist('b[i],i,0,n-1)])),
1212     sol:solve(eq,varlist),
1213     at('a[n-1]/dn,first(sol)))
1214   ), ratsimp(f)),
1215   if length(res)=1 then res[1] else res
1216 )$
1217
1218 /* ----------------------------------------------------------------------- */
1219 /* gain_optimum(fs,fr) - controller design according to the gain optimum   */
1220 /* ----------------------------------------------------------------------- */
1221 gain_optimum(fs,fr,[v]):=
1222 block(
1223   [fw,omega,vars,eq,pli,li,ratprint:false,assume_pos:true,res],
1224   fw:ratsimp(closed_loop(fr*fs)),
1225   if emptyp(v) then vars:delete(s,listofvars(fr)) else vars:flatten([v]),
1226   res:map(lambda([ff],
1227   block(
1228     ff:ratsimp(ev(ff,s=%i*omega)),
1229     eq:cabs(num(ff))**2-cabs(denom(ff))**2,
1230     li:coefficient_list(expand(eq),omega),
1231     li:sublist(li,lambda([u],not atom(u))),
1232     li:sublist(li,lambda([u],is(op(u)="+"))),
1233     li:sublist(li,lambda([v],not apply("and",map(lambda([u],freeof(u,v)),vars)))),
1234     li:makelist(li[k],k,1,min(length(li),length(vars))),
1235     li:solve(li,vars),
1236     pli:sublist(li,lambda([x],
1237                 (not member(false,map(lambda([u],is(rhs(u)>0)),x))))),
1238     if emptyp(pli) then li else first(pli))
1239   ), flatten([fw])),
1240   if length(res)=1 then res[1] else res
1241 )$
1242 /* ======================================================================= */
1243 /* Digital Algorithms                                                      */
1244 /* ======================================================================= */
1245
1246 /* ----------------------------------------------------------------------- */
1247 /* stepf(f,[tt]) - step function of f with the optional sampling time tt   */
1248 /* ----------------------------------------------------------------------- */
1249 stepf(f,[tt]):= /*11.1.2019 */
1250 (block([ratprint:false,T:1,index:1],
1251   if listp(tt) and length(tt)>0 then T:tt[1],
1252   if listp(f) then
1253   (
1254     if listp(f[1])  then (
1255       if t<f[1][1] then 0
1256       else (
1257         while t>f[index][1] and index<length(f) do index:index+1,
1258         if t>f[index][1] then f[index][2] else f[max(index-1,1)][2]))
1259     else f[max(1,min(floor(t/T)+1,length(f)))]
1260   )
1261   else subst(t=T*floor(t/T),f)))$
1262
1263 /* ----------------------------------------------------------------------- */
1264 /* algorithmp(f) - checks whether f is a digital algorithm                 */
1265 /* ----------------------------------------------------------------------- */
1266 algorithmp(f):=if  listp(f) and not mapatom(f[1]) and
1267   op(f[1])="=" and lhs(f[1])='(x[k])
1268   then true else false$
1269
1270 /* ----------------------------------------------------------------------- */
1271 /* eulerb(f) - Euler algorithm (with backward difference equations)        */
1272 /* ----------------------------------------------------------------------- */
1273 eulerb(F,[xx]):=
1274 block([ratprint:false,_Fz,_T,_zz,zz,nn,zz1,nn1,h,x,u,k,n],
1275   _T: if length(xx)>0 then first(xx) else T,
1276   _Fz:ratsimp(ev(F,s=(1/ _zz-1)/ _T*_zz)),
1277   zz:coefficient_list(num(_Fz),_zz),
1278   nn:coefficient_list(denom(_Fz),_zz),
1279   h:max(1,float(lmax(numbers(zz)))),
1280   zz:float(expand(zz/h)),
1281   nn:float(expand(nn/h)),
1282   zz1:apply("+",makelist(u[k-n],n,0,length(zz)-1)*zz),
1283   nn1:apply("+",makelist(x[k-n],n,0,length(nn)-1)*nn),
1284   h:coefficient_list(zz1-nn1,x[k]),
1285   [x[k]=-h[1]/h[2],_T]
1286 )$
1287
1288 /* ----------------------------------------------------------------------- */
1289 /* eulerf(f) - Euler algorithm (with forward difference equations)         */
1290 /* ----------------------------------------------------------------------- */
1291 eulerf(F,[xx]):=
1292 block([ratprint:false,_Fz,_T,_zz,zz,nn,zz1,nn1,h,x,u,k,n],
1293   _T: if length(xx)>0 then first(xx) else T,
1294   _Fz:ratsimp(ev(F,s=(1/ _zz-1)/ _T)),
1295   zz:coefficient_list(num(_Fz),_zz),
1296   nn:coefficient_list(denom(_Fz),_zz),
1297   h:max(1,float(lmax(numbers(zz)))),
1298   zz:float(expand(zz/h)),
1299   nn:float(expand(nn/h)),
1300   zz1:apply("+",makelist(u[k-n],n,0,length(zz)-1)*zz),
1301   nn1:apply("+",makelist(x[k-n],n,0,length(nn)-1)*nn),
1302   h:coefficient_list(zz1-nn1,x[k]),
1303   [x[k]=-h[1]/h[2],_T]
1304 )$
1305
1306 /* ----------------------------------------------------------------------- */
1307 /* tustin(f) - Tustin algorithm                                            */
1308 /* ----------------------------------------------------------------------- */
1309 tustin(F,[xx]):=
1310 block([ratprint:false,_Fz,_T,_zz,zz,nn,zz1,nn1,h,x,u,k,n],
1311   _T: if length(xx)>0 then first(xx) else T,
1312   _Fz:ratsimp(ev(F,s=2*(1/ _zz-1)/ _T/(1/ _zz+1))),
1313   zz:coefficient_list(num(_Fz),_zz),
1314   nn:coefficient_list(denom(_Fz),_zz),
1315   h:max(1,float(lmax(numbers(zz)))),
1316   zz:float(expand(zz/h)),
1317   nn:float(expand(nn/h)),
1318   zz1:apply("+",makelist(u[k-n],n,0,length(zz)-1)*zz),
1319   nn1:apply("+",makelist(x[k-n],n,0,length(nn)-1)*nn),
1320   h:coefficient_list(zz1-nn1,x[k]),
1321   [x[k]=-h[1]/h[2],_T]
1322 )$
1323
1324 /* ----------------------------------------------------------------------- */
1325 /* algorithm_step_response - Test function for an algorithm                */
1326 /* ----------------------------------------------------------------------- */
1327 algorithm_step_response(f,[nn]):=
1328 block([varli,alg:rhs(f[1]),n:100,x,u],
1329   local(x,u),
1330   if length(nn)>0 then n:nn[1],
1331   varli:sublist(ev(listofvars(alg),k=0),lambda([v],not atom(v) and mapatom(v)
1332     and (op(v)=x or op(v)=u))),
1333   for k:lmin(flatten(map(args,varli))) thru -1 do (u[k]:0,x[k]:0), /* Init. */
1334   for k:0 thru n do (u[k]:1,x[k]:ev(alg)),
1335   makelist(float([k*f[2],x[k]]),k,0,n)
1336 )$
1337
1338
1339 /* ======================================================================= */
1340 /* State Space                                                             */
1341 /* ======================================================================= */
1342
1343 /* ----------------------------------------------------------------------- */
1344 /* systemp(x) - checks whether x is a linear system with state matrices    */
1345 /* ----------------------------------------------------------------------- */
1346 systemp(x) :=
1347 block([nx,nu,ny],
1348 if not listp(x) then return(false),
1349 if length(x)<3 or length(x)>4 then return(false),
1350 if not matrixp(x[1]) then return(false),
1351 if not matrixp(x[2]) then return(false),
1352 if not matrixp(x[3]) then return(false),
1353 [nx,nu,ny]:[length(x[1]),length(x[2][1]),length(x[3])],
1354 if length(x[1][1]) # nx then return(false),
1355 if length(x[2]) # nx then return(false),
1356 if length(x[3][1]) # nx then return(false),
1357 if length(x)=3 then return(true),
1358 D: if not matrixp(x[4]) then matrix([x[4]]) else x[4],
1359 if length(D) # ny then return(false),
1360 if length(D[1]) # nu then return(false),
1361 return(true)
1362 )$
1363
1364 /* ----------------------------------------------------------------------- */
1365 /* nsystemp(x) - checks whether x is a linear system with numeric coeffs.  */
1366 /* ----------------------------------------------------------------------- */
1367 nsystemp(x) :=
1368 if systemp(x) and freeof(false,fullmap(numberp,x)) then true else false$
1369
1370 /* ----------------------------------------------------------------------- */
1371 /* controller_canonical_form                                               */
1372 /* ----------------------------------------------------------------------- */
1373 controller_canonical_form(f) :=
1374 block(
1375 [ratprint:false,zz,nn,zli,nli,A,B,C,D,o],
1376 [zz,nn]:[num(xthru(f)),denom(xthru(f))],
1377 zli:coefficient_list(zz,s), nli:coefficient_list(nn,s),
1378 zli:zli/last(nli), nli:nli/last(nli),
1379 o:length(nli)-1,
1380 zli:append(zli,makelist(0,k,1,length(nli)-length(zli))),
1381 A:genmatrix(lambda([u,v],if u=o then -nli[v] else if v-u=1 then 1 else 0),o,o),
1382 B:apply(matrix,endcons([1],makelist([0],k,1,o-1))),
1383 C:matrix(rest(map(lambda([u,v],u-last(zli)*v),zli,nli),-1)),
1384 D:last(zli),
1385 [A,B,C,D]
1386 );
1387
1388 /* ----------------------------------------------------------------------- */
1389 /* observer_canonical_form                                                 */
1390 /* ----------------------------------------------------------------------- */
1391 observer_canonical_form(f) :=
1392 block(
1393 [ratprint:false,zz,nn,zli,nli,A,B,C,D,o],
1394 [zz,nn]:[num(xthru(f)),denom(xthru(f))],
1395 zli:coefficient_list(zz,s), nli:coefficient_list(nn,s),
1396 zli:zli/last(nli), nli:nli/last(nli),
1397 o:length(nli)-1,
1398 zli:append(zli,makelist(0,k,1,length(nli)-length(zli))),
1399 A:genmatrix(lambda([u,v],if v=o then -nli[u] else if u-v=1 then 1 else 0),o,o),
1400 B:apply(matrix,rest(map(lambda([u,v],[u-last(zli)*v]),zli,nli),-1)),
1401 C:matrix(endcons(1,makelist(0,k,1,o-1))),
1402 D:last(zli),
1403 [A,B,C,D]
1404 );
1405
1406 /* ----------------------------------------------------------------------- */
1407 /* controllability_matrix                                                  */
1408 /* ----------------------------------------------------------------------- */
1409 controllability_matrix(x,[y]) :=
1410 block([A,B],
1411   [A,B]:if systemp(x) then [x[1],x[2]] else [x,first(y)],
1412   transpose(apply(matrix,makelist(
1413     flatten(args(transpose(A^^n.B))),n,0,length(A)-1))));
1414
1415 /* ----------------------------------------------------------------------- */
1416 /* observability_matrix                                                    */
1417 /* ----------------------------------------------------------------------- */
1418 observability_matrix(x,[y]) :=
1419 block([A,C],
1420   [A,C]:if systemp(x) then [x[1],x[3]] else [x,first(y)],
1421   apply(matrix,makelist(flatten(args(C.A^^n)),n,0,length(A)-1)));
1422
1423 /* ======================================================================= */
1424 /* Data Acquisition and Processing                                         */
1425 /* ======================================================================= */
1426
1427 /* ----------------------------------------------------------------------- */
1428 /* moving_average(x,n) - Moving Average of n elements in a list x          */
1429 /* ----------------------------------------------------------------------- */
1430 moving_average(x,n):=
1431 block([h1,h2:[],h],
1432   h1:append(makelist(first(x),k,n/2),x,makelist(last(x),k,n/2)),
1433   for i:1 thru length(x) do
1434   (
1435     h:makelist(h1[k],k,i,i+n-1),
1436     h2:endcons(float(apply("+",h)/n),h2)
1437   ),
1438   h2
1439 )$
1440
1441 /* ----------------------------------------------------------------------- */
1442 /* regan - regression analysis                  (Wilhelm Haager 8.11.2014) */
1443 /* ----------------------------------------------------------------------- */
1444 /* f ... list of functions                                                 */
1445 /* p ... list of points                                                    */
1446 /* ----------------------------------------------------------------------- */
1447 regan(p,f):=
1448   lambda([x],f.(transpose(x).x)^^(-1).transpose(x).matrix(map(second,p)))
1449   (apply(matrix,map(lambda([u],ev(f,'x=u)),map(first,p))))$
1450
1451 /* ----------------------------------------------------------------------- */
1452 /* add_noise(f,sd) - add a noise with standard deviation sd to f           */
1453 /* ----------------------------------------------------------------------- */
1454 add_noise(f,[opts]):=
1455 block([sd,res,resu,zli,nli],
1456   local(nrandom,sli),
1457   nrandom(x):=float(x*cos(random(1.0)*2*%pi)*sqrt(-2*log(random(1.0)))),
1458   sli(n):=makelist(s**k,k,0,n),
1459   sd:if length(opts)>0 then opts[1] else 0.1,
1460   if not listp(f) then f:[f],
1461   res:map(lambda([ff],
1462     if numberp(ff) or mapatom(ff) then resu:ff+nrandom(sd)
1463     else
1464     (
1465       zli:map(lambda([v],v+nrandom(sd)),coefficient_list(num(ff))),
1466       nli:map(lambda([v],v+nrandom(sd)),coefficient_list(denom(ff))),
1467       resu:apply("+",zli*sli(length(zli)-1))/apply("+",nli*sli(length(nli)-1))
1468     ),
1469   resu),f),
1470   if length(res)=1 then res[1] else res
1471 )$
1472
1473 /* ----------------------------------------------------------------------- */
1474 /* mult_noise(f,sd) - multiply a noise with standard deviation sd to f     */
1475 /* ----------------------------------------------------------------------- */
1476 mult_noise(f,[opts]):=
1477 block([sd,res,resu,zli,nli],
1478   local(nrandom,sli),
1479   nrandom(x):=float(x*cos(random(1.0)*2*%pi)*sqrt(-2*log(random(1.0)))),
1480   sli(n):=makelist(s**k,k,0,n),
1481   sd:if length(opts)>0 then opts[1] else 0.1,
1482   if not listp(f) then f:[f],
1483   res:map(lambda([ff],
1484     if numberp(ff) or mapatom(ff) then resu:ff*(1+nrandom(sd))
1485     else
1486     (
1487       zli:map(lambda([v],v*(1+nrandom(sd))),coefficient_list(num(ff))),
1488       nli:map(lambda([v],v*(1+nrandom(sd))),coefficient_list(denom(ff))),
1489       resu:apply("+",zli*sli(length(zli)-1))/apply("+",nli*sli(length(nli)-1))
1490     ),
1491   resu),f),
1492   if length(res)=1 then res[1] else res
1493 )$
1494
1495 /* ======================================================================= */
1496 /* Various Useful Functions                                                */
1497 /* ======================================================================= */
1498
1499 /* ----------------------------------------------------------------------- */
1500 /* sigma(t) - unit step function                                           */
1501 /* ----------------------------------------------------------------------- */
1502 /* sigma(t):=if t < 0 then 0 else 1$   */
1503
1504 /* ----------------------------------------------------------------------- */
1505 /* coefficient_list(p,s) - calculates a list of the coefficients of        */
1506 /*                         the polynomial p(s) in ascending order          */
1507 /* ----------------------------------------------------------------------- */
1508 /*coefficient_list(p,s) := lambda([v],map(lambda([u],coeff(v,s,u)),        */
1509 /*   makelist(i,i,0,hipow(v,s))))(expand(p))$                              */
1510 coefficient_list(p,[x]) :=
1511 block([var,varlist:listofvars(p)],
1512   if length(x)>0 then var:first(x)
1513   else if length(varlist)=1 then var:first(varlist)
1514   else if member('s,varlist) then var:s
1515   else if member('z,varlist) then var:z,
1516   lambda([v],map(lambda([u],coeff(v,var,u)),
1517   makelist(i,i,0,hipow(v,var))))(expand(p))
1518 )$
1519
1520 /* ----------------------------------------------------------------------- */
1521 /* get_option(o,l) - returns the option o from the list l (if existing)    */
1522 /* option_exists(o,l) - tests, whether the option o exists in list l       */
1523 /* list_option_exists(o,l) - tests, whether o exists and is a list         */
1524 /* delete_option(o,l) - returns a list with option o deleted from l        */
1525 /* set_option(o=v,l) - sets an option o with the value v                   */
1526 /* ----------------------------------------------------------------------- */
1527 get_option(o,l,[d]) := assoc(o,sublist(l,lambda([v],
1528        not atom(v) and op(v)="=")),if d=[] then false else first(d))$
1529 delete_option(o,l) ::= buildq([o,l],l:if option_exists(o,l) then
1530        sublist(l,lambda([u], atom(u) or op(u)#"=" or first(u)#o)) else l)$
1531 option_exists(o,l) := if member(o,map(first,sublist(l,lambda([v],
1532        not atom(v) and op(v)="=")))) then true$
1533 list_option_exists(o,l) :=
1534    if option_exists(o,l) then if listp(assoc(o,l)) then true$
1535 set_option(o,l) ::= buildq([o,l],l:endcons(o,delete_option(first(o),l)))$
1536
1537 chop(f):=
1538 block([epsilon:1.0e-10],
1539   local(test),
1540   test(u):=if numberp(u) and cabs(u)<epsilon then 0 else u,
1541   if mapatom(f) then test(f) else fullmap(test,f)
1542 )$
1543
1544 /* ----------------------------------------------------------------------- */
1545 /* // - parallel connection of resistances                 (WH, 26.9.2014) */
1546 /* ----------------------------------------------------------------------- */
1547 "//"([x]):=xthru(1/apply("+",1/x));
1548 nary("//",115);
1549
1550 /* ----------------------------------------------------------------------- */
1551 /* /_ - cis operator and cf operator for AC calculation    (WH, 10.1.2016) */
1552 /* ----------------------------------------------------------------------- */
1553 "/_"(a,b):=rectform(float(a*(cos(b*%pi/180)+%i*sin(b*%pi/180))))$
1554 infix("/_",150,150)$
1555 cisform(_z)::=buildq(
1556     [_a:ev(cabs(_z),numer,eval),_b:ev(carg(_z)*180/%pi,numer,eval)],
1557     '(_a /_ _b))$
1558 kill("//_");
1559 "//_"(_z)::=buildq(
1560     [_a:ev(cabs(_z),numer,eval),_b:ev(carg(_z)*180/%pi,numer,eval)],
1561     '(_a /_ _b))$
1562 postfix("//_",20);
1563
1564 /* ----------------------------------------------------------------------- */
1565 /* spartfrac - partial fraction expansion with numeric and symbolic poles  */
1566 /* ----------------------------------------------------------------------- */
1567 /* STILL UNDER DEVELOPMENT !
1568 spartfrac(G,[var]):=
1569 block([ratprint:false,variable,varli:listofvars(G)],
1570   if length(var)>0 then variable:var[1]
1571   else if member(s,varli) then variable:s
1572   else if member(z,varli) then variable:z
1573   else variable:first(varli),
1574   map(lambda([u],sum_form(u,3)),
1575   partfrac(num(G)/sallroots(denom(G)),variable))
1576 )$
1577 */
1578
1579 /* ----------------------------------------------------------------------- */
1580 /* sallroots - like allroots, but can contain symbolic factors and roots   */
1581 /* ----------------------------------------------------------------------- */
1582 /* STILL UNDER DEVELOPMENT !
1583 sallroots(p,[var]):=
1584 block([polyfactor:true,p1,pli,pli1,pli2,variable,varli:listofvars(G)],
1585   if length(var)>0 then variable:var[1]
1586   else if member(s,varli) then variable:s
1587   else if member(z,varli) then variable:z
1588   else variable:first(varli),
1589   p1:factorout(p,1),
1590   if op(p1)="*"
1591   then (
1592     pli:args(p1),
1593     pli1:sublist(pli,lambda([u],freeof(variable,u) or length(listofvars(u))>1)),
1594     pli2:sublist(pli,lambda([u],not freeof(variable,u) and not length(listofvars(u))>1)),
1595     apply("*",cons(allroots(apply("*",pli2)),pli1)))
1596   else allroots(p)
1597 )$
1598 */
1599
1600 /* ----------------------------------------------------------------------- */
1601 /* silt(f) - inverse Laplace transform with numeric and symbolic poles     */
1602 /* ----------------------------------------------------------------------- */
1603 /* STILL UNDER DEVELOPMENT !
1604 silt(f):=
1605 block([polyfactor:true,ratprint:false,ft,res],
1606   if not listp(f) then f:[f],
1607   res:map(lambda([ff],
1608     ft:ilt(ff,s,t),
1609     if string(op(ft))="ilt" then
1610     (
1611       ft:ilt(float(num(ff)/sallroots(denom(ff))),s,t),
1612       ft:ev(ft,float,expand)
1613     )
1614     else ft),
1615     xthru(f)),
1616     if length(res)=1 then res[1] else res
1617 )$
1618 */
1619
1620 /* ----------------------------------------------------------------------- */
1621 /* operators - list of all operators in an expr.     (Wilhelm Haager 2013) */
1622 /* ----------------------------------------------------------------------- */
1623 operators(x):=
1624 block([opli:[]],
1625 local(ops),
1626 ops(x):=if atom(x)
1627         then opli
1628         else (if not member(op(x),opli)
1629               then opli:cons(op(x),opli),map(ops,args(x))),ops(x),opli);
1630 contains(op,x):=if member(op,operators(x)) then true;
1631 contains_any(op,x):=
1632   if member(true,map(lambda([u],contains(u,x)),flatten([op])))
1633   then true;
1634 contains_all(op,x):=
1635   if member(false,map(lambda([u],contains(u,x)),flatten([op])))
1636   then false
1637   else true;
1638
1639 /* ----------------------------------------------------------------------- */
1640 /* numbers - list of all numbers in an expr.         (Wilhelm Haager 2019) */
1641 /* ----------------------------------------------------------------------- */
1642 numbers(x):=block([li:[]],
1643         fullmap(lambda([u],if numberp(u) then li:cons(u,li),u),x),li)$
1644
1645 /* ----------------------------------------------------------------------- */
1646 /* fullermap(f,x) - like fullmap, works also on atoms                      */
1647 /* ----------------------------------------------------------------------- */
1648 fullermap(f,x):=if mapatom(x) then f(x) else fullmap(f,x);
1649
1650 /* ----------------------------------------------------------------------- */
1651 /* periodic(f,x1,x2) - makes the function f periodic       (W.Haager 2019) */
1652 /* ----------------------------------------------------------------------- */
1653 periodic(f,[x]):=
1654 block([par,x1,x2,ratprint:false],
1655   if length(x)>=2 then (x1:x[1],x2:x[2])
1656      else if length(x)=1 then (x1:0,x2:x[1])
1657      else (x1:0,x2:1),
1658   par:first(listofvars(f)),
1659   subst(par=par-(x2-x1)*(floor((par-x1)/(x2-x1))),f)
1660 )$
1661
1662
1663
1664 /* ======================================================================= */
1665 /* Auxiliary Functions (only for internal use)                             */
1666 /* ======================================================================= */
1667
1668 /* ----------------------------------------------------------------------- */
1669 /* _COMA_bpr(f) - calculation of an appropriate plot range                 */
1670 /* ----------------------------------------------------------------------- */
1671 _COMA_bpr(f):=
1672 block(
1673   [frequencies,omin,omax,g,dt],
1674   [f,dt]:_COMA_strip(f),
1675   f:map(lambda([u],if not mapatom(u) and op(u)='asymptotic then first(args(u)) else u),f),
1676   g:gain_crossover(f),
1677   frequencies:flatten([_COMA_cf(sublist(f,lambda([u],not _COMA_go(u)))),
1678                   if length(g)>0 then map(last,flatten([g])) else []]),
1679   if length(frequencies)=0 then frequencies:[1],
1680   omin:apply(min,frequencies)/3,   /* ADAPT RANGE HERE */
1681   omax:apply(max,frequencies)*3,   /* ADAPT RANGE HERE */
1682   [10**floor(log(omin)/log(10)), 10**ceiling(log(omax)/log(10))]
1683 )$
1684
1685 /* ----------------------------------------------------------------------- */
1686 /* _COMA_srat - calculates an appropriate display time for step_response   */
1687 /* ----------------------------------------------------------------------- */
1688 _COMA_srat(f):=
1689 block(
1690   [pl,l,dt],
1691   [f,dt]:_COMA_strip(f),
1692   pl:flatten(poles(f)),
1693   pl:delete(0.0,pl), pl:delete(0,pl),
1694   if length(pl)=0 then 1 else
1695   block(
1696   l:map(lambda([u],float(5/sqrt(realpart(u)**2+imagpart(u)**2))),pl),
1697   apply(max,l))+lmax(-log(dt)/s)
1698 )$
1699
1700 /* ----------------------------------------------------------------------- */
1701 /* _COMA_npv - next proper value for a plot range                          */
1702 /* ----------------------------------------------------------------------- */
1703 _COMA_npv(f,li):=
1704 block(
1705   [n:length(li)],
1706   while f>last(li) do li:append(li,rest(li,length(li)-n)*10),
1707   li:append(li,rest(li,length(li)-n)*10), /* zur Vermeidung eines numer. Bugs */
1708   while f<first(li) do li:append(rest(li,n-length(li))/10,li),
1709   first(sublist(li,lambda([u],is(f<u))))
1710 )$
1711
1712 /* ----------------------------------------------------------------------- */
1713 /* _COMA_cf(f) - cutoff frequencies                                        */
1714 /* ----------------------------------------------------------------------- */
1715 _COMA_cf(f):=
1716 block([pl],
1717   if not listp(f) then f:[f],
1718   f:sublist(f, lambda([u],not _COMA_go(u))),
1719   pl:flatten([poles(f),zeros(f)]),
1720   pl:sort(map(lambda([u],sqrt(realpart(u)**2+imagpart(u)**2)),pl)),
1721   pl:sublist(pl,lambda([u],u>0))
1722 )$
1723
1724 /* ----------------------------------------------------------------------- */
1725 /* _COMA_mx - "map extended", like map, adapts all lists according         */
1726 /*   to the last   one, then performs a map over all lists                 */
1727 /* ----------------------------------------------------------------------- */
1728 _COMA_mx([lists]) :=
1729 block(
1730   [res:["["],i],
1731   for i:1 step 1 thru length(lists)-1 do
1732      res:endcons(
1733        block([n1:length(lists[i]),n2:length(last(lists)),res:[],j],
1734        for j:1 step 1 thru n2 do res:endcons(lists[i][mod(j-1,n1)+1],res),
1735        res),res),
1736   res:endcons(last(lists),res),
1737   apply(map,res)
1738 )$
1739
1740 /* ----------------------------------------------------------------------- */
1741 /* _COMA_strip(f) - strips off time delays from transfer functions.        */
1742 /* _COMA_sdt(f) - newr, better version (2019-02-18)                        */
1743 /* ----------------------------------------------------------------------- */
1744 /* Result is a list with 2 elements:                                       */
1745 /* 1. the transfer functions without time delay                            */
1746 /* 2. the according time delays                                            */
1747 /* ----------------------------------------------------------------------- */
1748 _COMA_strip(f):=
1749 block([h,d,fli:[],dli:[]],
1750   f:if listp(f) then f else[f],
1751   for e in f do
1752   (
1753      d:1,
1754      h:if not atom(e) then scanmap(lambda([u],if not atom(u) and op(u)="^" and first(args(u))=%e
1755         then (d:u, 1) else u),e)
1756         else (e),
1757      fli:endcons(h,fli),
1758      dli:endcons(d,dli)),
1759      [fli,dli]
1760 )$
1761 _COMA_sdt(f):= /* strip off time delays */
1762 block([df,dfli,element,fli:[],dli:[],ratprint:false],
1763   df:expandwrt(f,%e),
1764   dfli:if mapatom(df) or op(df)#"+" then [df] else args(df),
1765   fli:subst(%e=1,dfli),
1766   dli:ratsimp(dfli/fli),
1767   [fli,dli]
1768 )$
1769
1770 /* ----------------------------------------------------------------------- */
1771 /* _COMA_go(f) - checks whether f is a graphic object                      */
1772 /* ----------------------------------------------------------------------- */
1773 _COMA_go(f) :=
1774 if (not atom(f) and member(op(f),['explicit,'points,'implicit,'parametric,
1775              'region,'polar,'polygon,'rectangle,'ellipse,'label,'stepf]))
1776    then true else false$
1777
1778 /* ----------------------------------------------------------------------- */
1779 /* _COMA_ppo(f,defs,opts) - process plot options                           */
1780 /* ----------------------------------------------------------------------- */
1781 /* f ... List of functions, objects (or transfer functions)                */
1782 /*       f = 0 ... only global options are processed                       */
1783 /*       f = 1 ... options are processed, but no functions                 */
1784 /* defs ... default options for a routine                                  */
1785 /* opts ... explicitly given options for a routine                         */
1786 /* ----------------------------------------------------------------------- */
1787 _COMA_ppo(f,defs,opts):=
1788 block(
1789   [x,hli,fli,ar,preamble:[],list_options, list_exceptions:['dimensions,
1790      'xrange,'yrange,'trange,'user_preamble,'ip_grid,'ip_grid_in,'allocation]],
1791   if listp(coma_defaults)=false then coma_defaults:['grid=true],
1792   hli:append(coma_defaults,defs,opts),
1793 /* Gnuplot Praeambeln in Liste zusammenfassen */
1794   set_option('user_preamble=flatten(map(lambda([u],rhs(u)),
1795      sublist(hli,lambda([u],lhs(u)='user_preamble)))),hli),
1796   preamble:assoc(user_preamble,hli,[]),
1797 /* Mehrfache Optionen loeschen */
1798   for i:length(hli) step -1 thru 1 do hli:set_option(hli[i],hli),
1799 /* wx-Entscheidung */
1800   if option_exists('wx,hli) then
1801      block(if wxx#true and wxx#false then wxx:assoc('wx,hli),
1802         delete_option('wx,hli)),
1803   if option_exists('terminal,opts) then wxx:false,
1804   if option_exists('dimensions,opts) then wxplot_size:assoc('dimensions,hli,[]),
1805   preamble:assoc(user_preamble,hli,[]),
1806 /* Aspect Ratio in die Gnuplot Praeambel einbauen */
1807   if option_exists('aspect_ratio,hli) then block(
1808      ar:float(assoc('aspect_ratio,hli)),
1809      delete_option('aspect_ratio,hli),
1810      preamble:endcons(concat("set size ratio ",ar),flatten([preamble]))),
1811   if length(preamble)>0 then set_option(user_preamble=preamble,hli),
1812 /* Listenoptionen verarbeiten */
1813   list_options:sublist(hli,lambda([u],listp(rhs(u)))),
1814   for i:1 thru length(list_exceptions) do
1815       delete_option(list_exceptions[i],list_options),
1816   for i:1 thru length(list_options) do
1817       delete_option(lhs(list_options[i]),hli),
1818   if f=0 then return(hli),
1819   if f=1 then return(append(hli,list_options)),
1820   list_options:map(lambda([v],map(lambda([u],lhs(v)=u),rhs(v))),list_options),
1821   append(hli,apply(_COMA_mx,append(list_options,[f])))
1822 )$
1823
1824 /* ======================================================================= */
1825 /*   END of COMA                                                           */
1826 /* ======================================================================= */