share/diffequations/drawdf.mac

   1 /*
   2  * Copyright (C) 2010 Mark H. Weaver <mhw@netris.org>
   3  *
   4  * drawdf: Direction field plotting package for Maxima
   5  *
   6  * This program is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU General Public License
   8  * as published by the Free Software Foundation; either version 2
   9  * of the License, or (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  */
  20
  21 /*
  22  * If you encounter stack overflows when drawing complex plots,
  23  * try redefining draw-transform (from draw.lisp) as follows:
  24  *
  25  * (defun draw-transform (expr transform)
  26  *   (reduce #'append
  27  *           (mapcar #'(lambda (x) (draw-transform-one x transform))
  28  *                   expr)
  29  *           :from-end t))
  30  */
  31
  32 if not get('draw,'version) then load("draw") $
  33
  34 default_aspect_ratio: 4/3 $
  35 default_horiz_grid_steps: 30 $
  36 default_horiz_arrow_spots: 30 $
  37
  38 vec_mag(v) := sqrt(v.v)$
  39 normalize_vec(v) := v/sqrt(v.v)$
  40 if drawdf_translate='drawdf_translate then drawdf_translate:false;
  41 if drawdf_translate then translate(normalize_vec);
  42
  43 make_f_rk_4(dxdt_expr, dydt_expr, x_var, y_var, x_min, x_max, y_min, y_max, partition_expr) := (
  44   define(funmake('f_rk_4,vars),
  45     buildq(
  46       [ dxdt_expr, dydt_expr, dxdt_var:?gensym(), dydt_var:?gensym(),
  47         x_var, y_var, x_min, x_max, y_min, y_max, partition_expr ],
  48       block([dxdt_var, dydt_var],
  49         if x_var<x_min or x_var>x_max or y_var<y_min or y_var>y_max then error(),
  50         if partition_expr#current_partition then (
  51           if current_partition='init then current_partition:partition_expr
  52           else error()),
  53         dxdt_var:dxdt_expr,
  54         dydt_var:dydt_expr,
  55         if imagpart(dxdt_var)^2 > 0 or imagpart(dydt_var)^2 > 0 then error(),
  56         [dxdt_var,dydt_var]))),
  57   if drawdf_translate then translate(f_rk_4)) $
  58
  59 /* Derived from rk in dynamics.mac by Jaime E. Villate <villate@fe.up.pt> */
  60 calc_soln_curve(derivs, vars, init, mins, maxs, max_step, max_tstep, min_tstep, duration, max_nsteps, partition_expr) :=
  61 block (
  62   [ xxx,ddd,k1,k2,k3,k4,t0,t1,r,dxxx,nxxx,nk1,tstep,
  63     current_partition:'init,
  64     numer:true,display2d:false ],
  65
  66   make_f_rk_4(first(derivs), second(derivs), first(vars), second(vars),
  67     first(mins), first(maxs), second(mins), second(maxs), partition_expr),
  68   t0:0,
  69   xxx: init,
  70   r: errcatch ( k1: apply('f_rk_4,xxx) ),
  71   if length(r)=0 then return([]),
  72   ddd: [xxx],
  73   tstep: max_tstep,
  74   while max_nsteps > 0 and t0 < duration do (
  75     max_nsteps: max_nsteps-1,
  76     r: errcatch (
  77       t1: t0+tstep,
  78       k2: apply('f_rk_4,xxx+k1*tstep/2),
  79       k3: apply('f_rk_4,xxx+k2*tstep/2),
  80       k4: apply('f_rk_4,xxx+k3*tstep),
  81       dxxx: tstep*(k1+2*k2+2*k3+k4)/6,
  82       nxxx: xxx+dxxx,
  83       nk1: apply('f_rk_4,nxxx)),
  84     if length(r)=0 then (
  85       tstep: tstep/2,
  86       if tstep < min_tstep then return())
  87     else if (abs(first(dxxx)) > first(max_step) or
  88              abs(second(dxxx)) > second(max_step)) and
  89            tstep/2 >= min_tstep then (
  90       tstep: tstep/2)
  91     else (
  92       /* print("t:", t1, "dt:", tstep, "x:", nxxx, "dx:", dxxx), */
  93       tstep: max_tstep,
  94       k1: nk1,
  95       t0: t1,
  96       xxx: nxxx,
  97       ddd: cons(xxx, ddd))),
  98   reverse(ddd))$
  99 if drawdf_translate then translate(calc_soln_curve);
 100
 101 /*
 102  * Note: if dir='reverse, this will return the curve in
 103  *   reverse time order (starting from the initial point),
 104  *   in order to help solution_arrows() draw the arrows
 105  *   around saddle points more nicely.  In general, you
 106  *   should not rely on the direction of the returned curve.
 107  */
 108 calc_dir_soln_curve(derivs, vars, init, mins, maxs, max_step, dir, max_tstep, min_tstep, duration, max_nsteps, partition_expr) := (
 109   if elementp(dir, {'forward,'right}) then (
 110     calc_soln_curve(derivs, vars, init, mins, maxs, max_step, max_tstep, min_tstep, duration, max_nsteps, partition_expr))
 111   elseif elementp(dir, {'reverse,'left}) then (
 112     calc_soln_curve(-derivs, vars, init, mins, maxs, max_step, max_tstep, min_tstep, duration, max_nsteps, partition_expr))
 113   elseif dir = 'both then block(
 114     [ rev:calc_soln_curve(-derivs, vars, init, mins, maxs, max_step, max_tstep, min_tstep, duration, max_nsteps, partition_expr),
 115       fwd:calc_soln_curve(derivs, vars, init, mins, maxs, max_step, max_tstep, min_tstep, duration, max_nsteps, partition_expr) ],
 116     if fwd = [] then reverse(rev)
 117     else if rev = [] then fwd
 118     else append(reverse(rev), rest(fwd)))
 119   else error("Invalid direction:",dir))$
 120
 121 disallow_arrows_in(x_lo, y_lo, x_hi, y_hi) := block(
 122   [ x_step, y_step, a_spot_x, a_spot_y, b_spot_x, b_spot_y,
 123     a_spot_x_hi, a_spot_y_hi, b_spot_x_hi, b_spot_y_hi ],
 124
 125   x_step: (x_max-x_min) / horiz_arrow_spots,
 126   y_step: (y_max-y_min) / vert_arrow_spots,
 127
 128   calc_point_vars([x_hi,y_hi]),
 129   a_spot_x_hi:a_spot_x, a_spot_y_hi:a_spot_y,
 130   b_spot_x_hi:b_spot_x, b_spot_y_hi:b_spot_y,
 131   calc_point_vars([x_lo,y_lo]),
 132
 133   for y:a_spot_y thru a_spot_y_hi do (
 134     for x:a_spot_x thru a_spot_x_hi do (
 135       soln_arrow_a_spots[x,y]:2 )),
 136
 137   for y:b_spot_y thru b_spot_y_hi do (
 138     for x:b_spot_x thru b_spot_x_hi do (
 139       soln_arrow_b_spots[x,y]:2 )))$
 140
 141 solution_arrows(curve) := block(
 142   [ arrows:[], margin:0.5,
 143     spacing:12, start_space:6, end_space:4, edge_space:4,
 144     possible_arrow_dist:0, possible_arrow_pt:false,
 145     last_arrow_dist:0, new_dist, last_pt:first(curve),
 146     gr_x, gr_y, a_spot_x, a_spot_y, b_spot_x, b_spot_y,
 147     x_step, y_step, grid_step_vec, dv,
 148     numer:true ],
 149
 150   x_step: (x_max-x_min) / horiz_arrow_spots,
 151   y_step: (y_max-y_min) / vert_arrow_spots,
 152   grid_step_vec: [x_step,y_step] * 0.3,
 153
 154   define(funmake('f_calc_dv,[x_var,y_var]),
 155     float(grid_step_vec * normalize_vec([dxdt,dydt] / grid_step_vec))),
 156   if drawdf_translate then translate(f_calc_dv),
 157
 158   last_arrow_dist:spacing - start_space,
 159   for pt in curve do (
 160     calc_point_vars(pt),
 161     if (gr_x < margin or gr_x > horiz_arrow_spots - margin or
 162         gr_y < margin or gr_y > vert_arrow_spots - margin or
 163         soln_arrow_a_spots[a_spot_x,a_spot_y] = 2 or
 164         soln_arrow_b_spots[b_spot_x,b_spot_y] = 2) then (
 165       if (possible_arrow_pt # false and possible_arrow_dist >= edge_space) then (
 166         confirm_solution_arrow(possible_arrow_pt) ),
 167       last_pt:pt, possible_arrow_pt:false, last_arrow_dist:spacing - edge_space )
 168     else (
 169       new_dist:vec_mag((pt - last_pt) / grid_step_vec),
 170       last_pt:pt,
 171       possible_arrow_dist:possible_arrow_dist + new_dist,
 172       last_arrow_dist:last_arrow_dist + new_dist,
 173       if (soln_arrow_a_spots[a_spot_x,a_spot_y] = 1 or
 174           soln_arrow_b_spots[b_spot_x,b_spot_y] = 1) then (
 175         if (possible_arrow_pt # false and possible_arrow_dist < spacing) then (
 176           possible_arrow_pt:false ),
 177         last_arrow_dist:0 )
 178       elseif (possible_arrow_pt = false and last_arrow_dist >= spacing) then (
 179         possible_arrow_pt:pt, possible_arrow_dist:0 )
 180       elseif (possible_arrow_pt # false and possible_arrow_dist >= spacing) then (
 181         confirm_solution_arrow(possible_arrow_pt),
 182         possible_arrow_pt:pt, possible_arrow_dist:0 ))),
 183   if (possible_arrow_pt # false and possible_arrow_dist >= end_space) then (
 184     confirm_solution_arrow(possible_arrow_pt) ),
 185   arrows )$
 186
 187 calc_point_vars(pt) := (
 188   gr_x: (first(pt)-x_min)/x_step,
 189   gr_y: (second(pt)-y_min)/y_step,
 190   a_spot_x: fix(gr_x + 1),
 191   a_spot_y: fix(gr_y + 1),
 192   b_spot_x: fix(gr_x + 0.5),
 193   b_spot_y: fix(gr_y + 0.5) )$
 194
 195 confirm_solution_arrow(pt) := block (
 196   [ gr_x, gr_y, a_spot_x, a_spot_y, b_spot_x, b_spot_y ],
 197   calc_point_vars(pt),
 198   last_arrow_dist:possible_arrow_dist,
 199   soln_arrow_a_spots[a_spot_x,a_spot_y]:1,
 200   soln_arrow_b_spots[b_spot_x,b_spot_y]:1,
 201   dv:apply('f_calc_dv, pt),
 202   arrows:cons(funmake('vector, [pt-dv/2, dv]), arrows) )$
 203
 204 solution_plot(pt) := block( [ curve, duration, numer:true ],
 205   duration: (
 206     if soln_duration#false then soln_duration
 207     else if soln_nsteps#false then soln_nsteps*soln_tstep
 208     else 10.0),
 209   curve: calc_dir_soln_curve(
 210     [dxdt,dydt], [x_var,y_var], float(pt),
 211     [x_min-(x_max-x_min), y_min-(y_max-y_min)],
 212     [x_max+(x_max-x_min), y_max+(y_max-y_min)],
 213     [x_max-x_min, y_max-y_min]/100,
 214     soln_direction, soln_tstep, soln_tstep/300, duration, 100000, partition_expr),
 215   if curve=[] then curve:[pt],
 216   [ 'point_type='none, funmake('soln_curve, [curve, 'soln_arrows = soln_arrows]) ] )$
 217
 218 /*
 219  * Ideally, this should probably detect and remove overlaps in the curve.
 220  */
 221 equipotential_plot(pt) := block(
 222   [ dxdt:dydt, dydt:-dxdt, soln_arrows:false, soln_direction:'both ],
 223   solution_plot(pt) )$
 224
 225 linearize_at(pt, derivs, vars) := block(
 226   [ subs: map("=", vars, pt) ],
 227   subst(subs, apply('matrix, outermap('diff, derivs, vars))))$
 228
 229 /*
 230  * WARNING: this currently assumes you are running a new enough version
 231  * of Maxima to have the new improved eigenvectors function!
 232  */
 233 saddle_plot(pt) := block(
 234   [ eigen_info, eigvals, eigval, eigvects, the_graphics:[], soln_direction ],
 235
 236   eigen_info: block( [numer:false, keepfloat:false],
 237     float(eigenvectors(linearize_at(pt, [dxdt,dydt], [x_var,y_var])))),
 238
 239   eigvals: inpart(eigen_info,1,1),
 240   for i:1 thru length(eigvals) do (
 241     eigval: inpart(eigvals,i),
 242     if (imagpart(eigval) = 0 and eigval # 0) then (
 243       soln_direction: (if eigval > 0 then 'forward else 'reverse),
 244       eigvects: inpart(eigen_info,2,i),
 245       for eigvect in eigvects do (
 246         the_graphics: append(
 247           solution_plot(pt + 0.001 * normalize_vec(eigvect)),
 248           solution_plot(pt - 0.001 * normalize_vec(eigvect)),
 249           the_graphics)))),
 250   the_graphics)$
 251
 252 valid_var_lo_hi_spec(spec) := (
 253   is(listp(spec) and length(spec)=3 and symbolp(first(spec)) and
 254     numberp(float(second(spec))) and numberp(float(third(spec)))))$
 255
 256 list_of_points_from(params) := block([],
 257   params: args(params),
 258   if length(params) # 1 then return(params),
 259   params: first(params),
 260   if (length(params) # 2) or listp(first(params)) then return(params),
 261   return([params]))$
 262
 263 df_graphics(derivs_spec, [params]) := block(
 264   [ dxdt:1, dydt, x_var:'x, y_var:'y, x_min:-10, x_max:10, y_min:-10, y_max:10,
 265     aspect_ratio:'auto, horiz_grid_steps:'auto, vert_grid_steps:'auto,
 266     horiz_arrow_spots:'auto, vert_arrow_spots:'auto,
 267     soln_tinitial:0, soln_tstep:0.1, soln_nsteps:false, soln_duration:false,
 268     field_tstep:0.1, field_nsteps:false, field_duration:false,
 269     soln_direction:'both, soln_arrows:false,
 270     show_field:true, field_arrows:'auto, field_degree:1,
 271     field_horiz_offset:0, field_vert_offset:0, partition_expr:false,
 272     field_color:'auto, init_soln_color:'auto,
 273     field_head_size:0.5, field_head_type:'filled_circle,
 274     temp_params:[], param, pvar, pvar_table,
 275     x_step, y_step, grid_step_vec, the_graphics:[], pre_graphics:[],
 276     soln_arrow_a_spots, soln_arrow_b_spots,
 277     numer:true, keepfloat:true ],
 278   local(f_calc_dv),
 279
 280   if listp(derivs_spec) then (
 281     if length(derivs_spec) = 1 then [dydt]: derivs_spec
 282     elseif length(derivs_spec) = 2 then [dxdt,dydt]: derivs_spec
 283     else error("Invalid derivs_spec", derivs_spec))
 284   else dydt: derivs_spec,
 285
 286   if (length(params) >= 1 and listp(first(params)) and
 287       length(first(params)) = 2 and every('symbolp, first(params))) then (
 288     [x_var,y_var]: first(params),
 289     params: rest(params))
 290   elseif (length(params) >= 2 and valid_var_lo_hi_spec(first(params)) and
 291           valid_var_lo_hi_spec(second(params))) then (
 292     [x_var,x_min,x_max]: float(first(params)),
 293     [y_var,y_min,y_max]: float(second(params)),
 294     params: rest(params,2)),
 295
 296   while params#[] do (
 297     param: first(params),
 298     params: rest(params),
 299     if (operatorp(param, "=") and elementp(lhs(param), {'aspect_ratio,
 300         'horiz_grid_steps,'vert_grid_steps,'horiz_arrow_spots,'vert_arrow_spots}))
 301     then lhs(param) :: float(rhs(param))
 302     else if (operatorp(param, "=") and lhs(param)='field_grid and
 303              listp(rhs(param)) and length(rhs(param))=2)
 304     then [horiz_grid_steps,vert_grid_steps]: float(rhs(param))
 305     else if (listp(param) and first(param) = x_var and length(param) = 3)
 306     then [x_min,x_max]: float(rest(param))
 307     else if (listp(param) and first(param) = y_var and length(param) = 3)
 308     then [y_min,y_max]: float(rest(param))
 309     else if (listp(param) and not atom(first(param)))
 310     then params: append(param, params)
 311     else temp_params: cons(param, temp_params) ),
 312   params: reverse(temp_params),
 313
 314   if aspect_ratio = 'auto then aspect_ratio:(
 315     if (horiz_grid_steps # 'auto and vert_grid_steps # 'auto)
 316     then horiz_grid_steps / vert_grid_steps
 317     elseif (horiz_arrow_spots # 'auto and vert_arrow_spots # 'auto)
 318     then horiz_arrow_spots / vert_arrow_spots
 319     else default_aspect_ratio ),
 320
 321   if horiz_grid_steps = 'auto then horiz_grid_steps:(
 322     if vert_grid_steps # 'auto
 323     then fix(vert_grid_steps * aspect_ratio)
 324     else default_horiz_grid_steps ),
 325   if vert_grid_steps = 'auto
 326   then vert_grid_steps:fix(horiz_grid_steps / aspect_ratio),
 327
 328   if horiz_arrow_spots = 'auto then horiz_arrow_spots:(
 329     if vert_arrow_spots # 'auto
 330     then fix(vert_arrow_spots * aspect_ratio)
 331     else default_horiz_arrow_spots ),
 332   if vert_arrow_spots = 'auto
 333   then vert_arrow_spots:fix(horiz_arrow_spots / aspect_ratio),
 334
 335   soln_arrow_a_spots: make_array(fixnum, horiz_arrow_spots+2, vert_arrow_spots+2),
 336   soln_arrow_b_spots: make_array(fixnum, horiz_arrow_spots+2, vert_arrow_spots+2),
 337
 338   pvar_table: '[ direction        = soln_direction,
 339                  dir              = soln_direction,
 340                  tinitial         = soln_tinitial,
 341                  tstep            = soln_tstep,
 342                  nsteps           = soln_nsteps,
 343                  duration         = soln_duration,
 344                  field_tstep      = field_tstep,
 345                  field_nsteps     = field_nsteps,
 346                  field_duration   = field_duration,
 347                  soln_arrows      = soln_arrows,
 348                  show_field       = show_field,
 349                  partition_expr   = partition_expr,
 350                  field_horiz_offset = field_horiz_offset,
 351                  field_vert_offset = field_vert_offset,
 352                  field_color      = field_color,
 353                  field_arrows     = field_arrows,
 354                  field_head_size  = field_head_size,
 355                  field_head_type  = field_head_type,
 356                  field_degree     = field_degree ],
 357
 358   for param in params do (
 359     if listp(param) then (
 360       if (pvar:assoc(first(param),pvar_table)) # false then pvar :: second(param)
 361       elseif (first(param) = 'trajectory_at and length(param) >= 3)
 362       then pre_graphics: append(reverse(solution_plot(rest(param))), pre_graphics)
 363       else error("Invalid parameter:", param))
 364     else if operatorp(param, "=") then (
 365       if (pvar:assoc(lhs(param),pvar_table)) # false then pvar :: rhs(param)
 366       else pre_graphics: cons(param, pre_graphics))
 367     else if operatorp(param, 'point_at)
 368     then pre_graphics: append(reverse([ 'point_type = 'filled_circle, funmake('points, [[args(param)]]) ]), pre_graphics)
 369     else if operatorp(param, 'points_at)
 370     then for pt in list_of_points_from(param) do (
 371       pre_graphics: append(reverse([ 'point_type = 'filled_circle, funmake('points, [[pt]]) ]), pre_graphics))
 372     else if operatorp(param, 'soln_at)
 373     then pre_graphics: append(reverse(solution_plot(args(param))), pre_graphics)
 374     else if operatorp(param, 'solns_at)
 375     then for pt in list_of_points_from(param) do pre_graphics: append(reverse(solution_plot(pt)), pre_graphics)
 376     else if operatorp(param, 'equipot_at)
 377     then pre_graphics: append(reverse(equipotential_plot(args(param))), pre_graphics)
 378     else if operatorp(param, 'equipots_at)
 379     then for pt in list_of_points_from(param) do pre_graphics: append(reverse(equipotential_plot(pt)), pre_graphics)
 380     else if operatorp(param, 'saddle_at)
 381     then pre_graphics: append(reverse(saddle_plot(args(param))), pre_graphics)
 382     else if operatorp(param, 'saddles_at)
 383     then for pt in list_of_points_from(param) do pre_graphics: append(reverse(saddle_plot(pt)), pre_graphics)
 384     else if operatorp(param, 'no_arrows_in)
 385     then apply('disallow_arrows_in, args(param))
 386     else pre_graphics: cons(param, pre_graphics)),
 387
 388   for element in pre_graphics do (
 389     if operatorp(element, 'soln_curve) then (
 390       the_graphics: cons(funmake('points, [first(element)]), the_graphics),
 391       if assoc('soln_arrows, rest(args(element)), false) then (
 392         the_graphics: append(solution_arrows(first(element)), the_graphics)))
 393     else the_graphics: cons(element, the_graphics)),
 394
 395   if soln_arrows or not show_field then (
 396     if field_arrows = 'auto then field_arrows: false,
 397     if field_color = 'auto then field_color: red,
 398     if init_soln_color = 'auto then init_soln_color: black )
 399   else (
 400     if field_arrows = 'auto then field_arrows: true,
 401     if field_color = 'auto then field_color: black,
 402     if init_soln_color = 'auto then init_soln_color: red ),
 403
 404   x_step: (x_max-x_min) / horiz_grid_steps,
 405   y_step: (y_max-y_min) / vert_grid_steps,
 406   grid_step_vec: [x_step,y_step] * (if field_degree#1 then 0.7 else if field_arrows then 0.8 else 0.5),
 407
 408   the_graphics: append(
 409     [ 'point_type   = 'none,
 410       'point_size   = 1,
 411       'color        = init_soln_color,
 412       'head_length  = 0.6 * x_step,
 413       'head_angle   = 30 ],
 414     the_graphics ),
 415
 416   if show_field then block ( [numer:true],
 417     if field_degree=1 then (
 418       define(funmake('f_calc_dv,[x_var,y_var]),
 419         float(grid_step_vec * normalize_vec([dxdt,dydt] / grid_step_vec))),
 420       if drawdf_translate then translate(f_calc_dv),
 421       for x: (x_min + x_step/2 + field_horiz_offset) step x_step thru x_max do (
 422         for y: (y_min + y_step/2 + field_vert_offset) step y_step thru y_max do (
 423           the_graphics: append(
 424             errcatch(
 425               block([v:[x,y], dv:apply('f_calc_dv,[x,y]), rdv, idv],
 426                 rdv: realpart(dv), idv: imagpart(dv),
 427                 if idv.idv > ratepsilon then error()
 428                 else if field_arrows then funmake('vector, [v-rdv/2, rdv])
 429                 else funmake('points, [[v-rdv/2, v+rdv/2]]))),
 430             the_graphics))))
 431     else if field_degree=2 then block (
 432       [ dxdt2: diff(dxdt,x_var)*dxdt + diff(dxdt,y_var)*dydt,
 433         dydt2: diff(dydt,x_var)*dxdt + diff(dydt,y_var)*dydt,
 434         t_var:?gensym() ],
 435       define(funmake('f_calc_dv,[x_var,y_var]),
 436         float([dxdt,dydt])),
 437       define(funmake('f_calc_dv2,[x_var,y_var]),
 438         float([dxdt2,dydt2])),
 439       if drawdf_translate then translate(f_calc_dv),
 440       if drawdf_translate then translate(f_calc_dv2),
 441       for x: (x_min + x_step/2 + field_horiz_offset) step x_step thru x_max do (
 442         for y: (y_min + y_step/2 + field_vert_offset) step y_step thru y_max do (
 443           the_graphics: append(
 444             apply('append,
 445               errcatch(
 446                 block(
 447                   [ v:[x,y], dv:apply('f_calc_dv,[x,y]), rdv, idv,
 448                     dv2:apply('f_calc_dv2,[x,y]), rdv2, idv2,
 449                     curve, t_min:false, t_max:false, programmode:true ],
 450                   rdv: realpart(dv), idv: imagpart(dv),
 451                   rdv2: realpart(dv2), idv2: imagpart(dv2),
 452                   if idv.idv > ratepsilon or idv2.idv2 > ratepsilon then error(),
 453                   curve: rdv2*t_var^2/2 + rdv*t_var,
 454                   for unit_vec in [[1,0],[-1,0],[0,1],[0,-1]] do (
 455                     for root in map('float,map('rhs,realroots(((unit_vec.(curve/grid_step_vec))-0.5)))) do (
 456                       if root>0 then (
 457                         if t_max=false or t_max>root then t_max:root
 458                       ) else if root<0 then (
 459                         if t_min=false or t_min<root then t_min:root
 460                       ))),
 461                   append(
 462                     [ funmake('parametric, append(curve + v, [t_var,t_min,t_max])) ],
 463                     if field_arrows then [ point_type = field_head_type, funmake('points, [[at(curve+v, [t_var=t_max])]]) ] else [])))),
 464             the_graphics))))
 465     else if field_degree='solns then block (
 466       [ curve, duration ],
 467       duration: (
 468         if field_duration#false then field_duration
 469         else if field_nsteps#false then field_nsteps*field_tstep
 470         else 10.0),
 471       for x: (x_min + x_step/2 + field_horiz_offset) step x_step thru x_max do (
 472         for y: (y_min + y_step/2 + field_vert_offset) step y_step thru y_max do (
 473           the_graphics: append(
 474             apply('append,
 475               errcatch(
 476                 for fudge in [[0,0],[1,1],[-1,1],[-1,-1],[1,-1],'done] do (
 477                   if fudge='done then error(),
 478                   curve:reverse(
 479                     calc_dir_soln_curve([dxdt,dydt], [x_var,y_var], [x,y]+fudge*grid_step_vec*1e-8,
 480                       [x,y]-grid_step_vec/2, [x,y]+grid_step_vec/2,
 481                       grid_step_vec/5, 'both,
 482                       field_tstep, field_tstep/300, duration, 100, partition_expr)),
 483                   if curve#[] then return(curve)),
 484                 append(
 485                   [ point_type = 'none, funmake('points, [curve]) ],
 486                   if field_arrows then [ point_type = field_head_type, funmake('points, [[first(curve)]]) ] else []))),
 487             the_graphics))))
 488     else error("field_degree must be either 1, 2, or 'solns")),
 489
 490   the_graphics: append(
 491     [ 'head_length   = 0.3 * x_step,
 492       'head_angle    = 20,
 493       'point_type    = 'none,
 494       'point_size    = field_head_size,
 495       'points_joined = true,
 496       'xrange        = [x_min, x_max],
 497       'yrange        = [y_min, y_max],
 498       'color         = field_color ],
 499     the_graphics),
 500
 501   apply('gr2d, the_graphics))$
 502
 503 drawdf([params]) := (draw(apply('df_graphics, params)), 0)$
 504 wxdrawdf([params]) := (wxdraw(apply('df_graphics, params)), 0)$
 505
 506 /*
 507
 508 /* Examples */
 509
 510 drawdf(x^2+y^2, [x,-5,5], [y,-5,5], solns_at([1,1], [2,1], [3,4]))$
 511
 512 drawdf([y,-9*sin(x)-y/5], [x,1,5], [y,-2,2])$
 513 drawdf([y,-9*sin(x)-y/5], [x,1,5], [y,-2,2], field_degree=2)$
 514 drawdf([y,-9*sin(x)-y/5], [x,1,5], [y,-2,2], tstep=0.01, soln_at(3.14,0))$
 515
 516 drawdf(y, [x,-2,2], [y,-10,10], tstep=0.05, duration=3, solns_at([0,1], [0,-1], [0,0]))$
 517 drawdf(y, [x,-2,2], [y,-10,10], tstep=0.05, duration=3, soln_arrows=true, solns_at([0,1], [0,-1], [0,0]))$
 518
 519 drawdf([x-x*y/2, (x*y - 3*y)/4], [x,2.5,3.5], [y,1.5,2.5], solns_at([3,1.8],[3,1.85],[3,1.9],[3,1.95])) $
 520
 521 drawdf(exp(-t)+y, [t,y], [t,-5,5], [y,-20,20], soln_at(0,-0.5))$
 522
 523 drawdf(3*sin(t)+1+y, [t,y], solns_at([0,-2.6],[0,-2.4]), color=blue, soln_at(0,-2.5))$
 524 drawdf(3*sin(t)+1+y, [t,y], soln_arrows=true, solns_at([0,-2.6],[0,-2.4],[0,-2.5]))$
 525
 526 /* competing species */
 527 drawdf([x*(1-x-y), y*(3/4-y-x/2)], [x,0,1.1], [y,0,1], duration=40, soln_arrows=true, point_at(1/2,1/2), solns_at([0.1,0.2], [0.2,0.1], [1,0.8], [0.8,1], [0.1,0.1], [0.6,0.05], [0.05,0.4], [1,0.01], [0.01,0.75])) $
 528
 529 /* damped pendulum phase portraits */
 530 drawdf([y,-9*sin(x)-y/5], tstep=0.05, soln_arrows=true, saddles_at([%pi,0], [-%pi,0]))$
 531 drawdf([y,-9*sin(x)-y/5], tstep=0.05, show_field=false, soln_arrows=true, saddles_at([3*%pi,0], [-3*%pi,0], [%pi,0], [-%pi,0]))$
 532
 533 /* direction fields sampled quadratically and numerically */
 534 drawdf([x*(1-x-y), y*(3/4-y-x/2)], [x,0,1.1], [y,0,1], duration=40, soln_arrows=true, field_degree=2, point_at(1/2,1/2), solns_at([0.1,0.2], [0.2,0.1], [1,0.8], [0.8,1], [0.1,0.1], [0.6,0.05], [0.05,0.4], [1,0.01], [0.01,0.75])) $
 535 drawdf([x*(1-x-y), y*(3/4-y-x/2)], [x,0,1.1], [y,0,1], duration=40, soln_arrows=true, field_degree='solns, point_at(1/2,1/2), solns_at([0.1,0.2], [0.2,0.1], [1,0.8], [0.8,1], [0.1,0.1], [0.6,0.05], [0.05,0.4], [1,0.01], [0.01,0.75])) $
 536
 537 drawdf([y,-9*sin(x)-y/5], tstep=0.05, soln_arrows=true, field_degree=2, saddles_at([%pi,0], [-%pi,0]))$
 538 drawdf([y,-9*sin(x)-y/5], tstep=0.05, soln_arrows=true, field_degree='solns, saddles_at([%pi,0], [-%pi,0]))$
 539
 540 */