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