share/draw/drawutils.mac

   1 /*               COPYRIGHT NOTICE
   2
   3 Copyright (C) 2010 Donald J Bindner
   4
   5 This program is free software; you can redistribute
   6 it and/or modify it under the terms of the
   7 GNU General Public License as published by
   8 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
  12 will be useful, but WITHOUT ANY WARRANTY;
  13 without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 GNU General Public License for more details at
  16 http://www.gnu.org/copyleft/gpl.html
  17 */
  18
  19
  20 /*
  21 This is a set of user contributed plotting routines
  22 based on package draw.
  23
  24 See documentation in drawutils.texi
  25 */
  26
  27 ?defvar(drawutils_version, 1) $
  28
  29 if not $draw_version then load("draw") $
  30
  31
  32
  33
  34
  35 /*********************/
  36 /* Vector fields     */
  37 /*       by          */
  38 /* Donald J. Bindner */
  39 /*      2010         */
  40 /*********************/
  41
  42
  43 plot_vector_field( F, X, Y, [args] ) := block(
  44  [vars, G,P,Q, points, grid_d,u,scale, range,xrange,yrange, v],
  45  scale : assoc('scale,args,1), args:delete( 'scale=scale, args ),
  46  /* create function versions of each componont and of the field itself */
  47  vars : [ X[1],Y[1] ],
  48  P : G( vars, ev(F[1])), P : subst( lambda, G, P ),
  49  Q : G( vars, ev(F[2])), Q : subst( lambda, G, Q ),
  50  G : lambda( [z], [apply(P,z), apply(Q,z)] ),
  51  /* create a list of points to base arrows at */
  52  points : listify( cartesian_product(
  53       setify( makelist(  X[2] + (X[3]-X[2])*i/10, i, 0, 10 )),
  54       setify( makelist(  Y[2] + (Y[3]-Y[2])*i/10, i, 0, 10 )) )),
  55  /* diagonal length of a grid square */
  56  grid_d : sqrt((X[3]-X[2])^2 + (Y[3]-Y[2])^2)/10.0,
  57  /* u is the divisor that shortens arrows to fit in grid squares */
  58  u : max( apply( max, makelist( abs(P(k[1],k[2])), k, points )) / (X[3]-X[2]) * 10,
  59           apply( max, makelist( abs(Q(k[1],k[2])), k, points )) / (Y[3]-Y[2]) * 10 ),
  60  if scale#0 and u>grid_d/50.0 then u:scale/u else u:1,
  61  /* xrange and yrange need to be large enough to contain head
  62   * and tail of each arrow (or some won't appear in output) */
  63  range : lambda( [z], [ apply(min,z), apply(max,z) ] ),
  64  xrange : range( flatten(makelist( [k[1],k[1]+P(k[1],k[2])], k, points ))),
  65  yrange : range( flatten(makelist( [k[2],k[2]+Q(k[1],k[2])], k, points ))),
  66  /* generate the vectors to draw and set an appropriate arrow
  67   * head length for each one */
  68  v : flatten( makelist( [head_length=u*sqrt(G(k).G(k))/5+grid_d/100, vector(k, u*G(k) )], k, points )),
  69  /* draw! */
  70  draw2d(color=blue,line_width=1,head_type='nofilled,
  71         head_angle=20,'xrange=xrange,'yrange=yrange,
  72         xlabel=string(X[1]),ylabel=string(Y[1]), args, v )   )$
  73
  74
  75 plot_vector_field3d( F, X, Y, Z, [args] ) := block(
  76  [vars, G,P,Q,R, points, grid_d,u,scale, range,xrange,yrange,zrange, v],
  77  scale : assoc('scale,args,1), args:delete( 'scale=scale, args ),
  78  /* create function versions of each componont and of the field itself */
  79  vars : [ X[1],Y[1],Z[1] ],
  80  P : G( vars, ev(F[1])), P : subst( lambda, G, P ),
  81  Q : G( vars, ev(F[2])), Q : subst( lambda, G, Q ),
  82  R : G( vars, ev(F[3])), R : subst( lambda, G, R ),
  83  G : lambda( [z], [apply(P,z), apply(Q,z), apply(R,z)] ),
  84  /* create a list of points to base arrows at */
  85  points : listify( cartesian_product(
  86       setify( makelist(  X[2] + (X[3]-X[2])*i/5, i, 0, 5 )),
  87       setify( makelist(  Y[2] + (Y[3]-Y[2])*i/5, i, 0, 5 )),
  88       setify( makelist(  Z[2] + (Z[3]-Z[2])*i/5, i, 0, 5 )) )),
  89  /* the diagonal length of a grid square */
  90  grid_d : sqrt((X[3]-X[2])^2 + (Y[3]-Y[2])^2 + (Z[3]-Z[2])^2)/5.0,
  91  /* u is the divisor that shortens arrows to fit in grid squares */
  92  u : max( apply( max, makelist( abs(P(k[1],k[2],k[3])), k, points )) / (X[3]-X[2]) * 5,
  93           apply( max, makelist( abs(Q(k[1],k[2],k[3])), k, points )) / (Y[3]-Y[2]) * 5,
  94           apply( max, makelist( abs(R(k[1],k[2],k[3])), k, points )) / (Z[3]-Z[2]) * 5 ),
  95  if scale#0 and u>grid_d/50.0 then u:scale/u else u:1,
  96  /* xrange,yrange,zrange need to be large enough to contain head
  97   * and tail of each arrow (or some won't appear in output) */
  98  range : lambda( [z], [ apply(min,z), apply(max,z) ] ),
  99  xrange : range( flatten(makelist( [k[1],k[1]+P(k[1],k[2],k[3])], k, points ))),
 100  yrange : range( flatten(makelist( [k[2],k[2]+Q(k[1],k[2],k[3])], k, points ))),
 101  zrange : range( flatten(makelist( [k[3],k[3]+R(k[1],k[2],k[3])], k, points ))),
 102  /* generate the vectors to draw and set an appropriate arrow
 103   * head length for each one */
 104  v : flatten(makelist( [head_length=u*sqrt(G(k).G(k))/10+grid_d/100, vector(k, u*G(k))], k, points )),
 105  /* draw! */
 106  draw3d(color=blue,line_width=1,head_type='nofilled,
 107         head_angle=10,'xrange=xrange,'yrange=yrange,
 108         'zrange=zrange,xlabel=string(X[1]),ylabel=string(Y[1]),zlabel=string(Z[1]), args, v ) )$
 109
 110