share/draw/drawutils.mac

   1 /*               COPYRIGHT NOTICE
   2
   3 Copyright (C) 2010 Donald J Bindner
   4               2015 Pankaj Sejwal
   5
   6 This program is free software; you can redistribute
   7 it and/or modify it under the terms of the
   8 GNU General Public License as published by
   9 the Free Software Foundation; either version 2
  10 of the License, or (at your option) any later version.
  11
  12 This program is distributed in the hope that it
  13 will be useful, but WITHOUT ANY WARRANTY;
  14 without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 GNU General Public License for more details at
  17 http://www.gnu.org/copyleft/gpl.html
  18 */
  19
  20
  21 /*
  22 This is a set of user contributed plotting routines
  23 based on package draw.
  24
  25 See documentation in drawutils.texi
  26 */
  27
  28 drawutils_version :  1 $
  29
  30 if draw_version = 'draw_version then load("draw") $
  31
  32
  33
  34
  35
  36 /*********************/
  37 /* Vector fields     */
  38 /*       by          */
  39 /* Donald J. Bindner */
  40 /*      2010         */
  41 /*********************/
  42
  43
  44 plot_vector_field( F, X, Y, [args] ) := block(
  45  [vars, G,P,Q, points, grid_d,u,scale, range,xrange,yrange, v],
  46  scale : assoc('scale,args,1), args:delete( 'scale=scale, args ),
  47  /* create function versions of each componont and of the field itself */
  48  vars : [ X[1],Y[1] ],
  49  P : G( vars, ev(F[1])), P : subst( lambda, G, P ),
  50  Q : G( vars, ev(F[2])), Q : subst( lambda, G, Q ),
  51  G : lambda( [z], [apply(P,z), apply(Q,z)] ),
  52  /* create a list of points to base arrows at */
  53  points : listify( cartesian_product(
  54       setify( makelist(  X[2] + (X[3]-X[2])*i/10, i, 0, 10 )),
  55       setify( makelist(  Y[2] + (Y[3]-Y[2])*i/10, i, 0, 10 )) )),
  56  /* diagonal length of a grid square */
  57  grid_d : sqrt((X[3]-X[2])^2 + (Y[3]-Y[2])^2)/10.0,
  58  /* u is the divisor that shortens arrows to fit in grid squares */
  59  u : max( apply( max, makelist( abs(P(k[1],k[2])), k, points )) / (X[3]-X[2]) * 10,
  60           apply( max, makelist( abs(Q(k[1],k[2])), k, points )) / (Y[3]-Y[2]) * 10 ),
  61  if scale#0 and u>grid_d/50.0 then u:scale/u else u:1,
  62  /* xrange and yrange need to be large enough to contain head
  63   * and tail of each arrow (or some won't appear in output) */
  64  range : lambda( [z], [ apply(min,z), apply(max,z) ] ),
  65  xrange : range( flatten(makelist( [k[1],k[1]+P(k[1],k[2])], k, points ))),
  66  yrange : range( flatten(makelist( [k[2],k[2]+Q(k[1],k[2])], k, points ))),
  67  /* generate the vectors to draw and set an appropriate arrow
  68   * head length for each one */
  69  v : flatten( makelist( [head_length=u*sqrt(G(k).G(k))/5+grid_d/100, vector(k, u*G(k) )], k, points )),
  70  /* draw! */
  71  draw2d(color=blue,line_width=1,head_type='nofilled,
  72         head_angle=20,'xrange=xrange,'yrange=yrange,
  73         xlabel=string(X[1]),ylabel=string(Y[1]), args, v )   )$
  74
  75
  76 plot_vector_field3d( F, X, Y, Z, [args] ) := block(
  77  [vars, G,P,Q,R, points, grid_d,u,scale, range,xrange,yrange,zrange, v],
  78  scale : assoc('scale,args,1), args:delete( 'scale=scale, args ),
  79  /* create function versions of each componont and of the field itself */
  80  vars : [ X[1],Y[1],Z[1] ],
  81  P : G( vars, ev(F[1])), P : subst( lambda, G, P ),
  82  Q : G( vars, ev(F[2])), Q : subst( lambda, G, Q ),
  83  R : G( vars, ev(F[3])), R : subst( lambda, G, R ),
  84  G : lambda( [z], [apply(P,z), apply(Q,z), apply(R,z)] ),
  85  /* create a list of points to base arrows at */
  86  points : listify( cartesian_product(
  87       setify( makelist(  X[2] + (X[3]-X[2])*i/5, i, 0, 5 )),
  88       setify( makelist(  Y[2] + (Y[3]-Y[2])*i/5, i, 0, 5 )),
  89       setify( makelist(  Z[2] + (Z[3]-Z[2])*i/5, i, 0, 5 )) )),
  90  /* the diagonal length of a grid square */
  91  grid_d : sqrt((X[3]-X[2])^2 + (Y[3]-Y[2])^2 + (Z[3]-Z[2])^2)/5.0,
  92  /* u is the divisor that shortens arrows to fit in grid squares */
  93  u : max( apply( max, makelist( abs(P(k[1],k[2],k[3])), k, points )) / (X[3]-X[2]) * 5,
  94           apply( max, makelist( abs(Q(k[1],k[2],k[3])), k, points )) / (Y[3]-Y[2]) * 5,
  95           apply( max, makelist( abs(R(k[1],k[2],k[3])), k, points )) / (Z[3]-Z[2]) * 5 ),
  96  if scale#0 and u>grid_d/50.0 then u:scale/u else u:1,
  97  /* xrange,yrange,zrange need to be large enough to contain head
  98   * and tail of each arrow (or some won't appear in output) */
  99  range : lambda( [z], [ apply(min,z), apply(max,z) ] ),
 100  xrange : range( flatten(makelist( [k[1],k[1]+P(k[1],k[2],k[3])], k, points ))),
 101  yrange : range( flatten(makelist( [k[2],k[2]+Q(k[1],k[2],k[3])], k, points ))),
 102  zrange : range( flatten(makelist( [k[3],k[3]+R(k[1],k[2],k[3])], k, points ))),
 103  /* generate the vectors to draw and set an appropriate arrow
 104   * head length for each one */
 105  v : flatten(makelist( [head_length=u*sqrt(G(k).G(k))/10+grid_d/100, vector(k, u*G(k))], k, points )),
 106  /* draw! */
 107  draw3d(color=blue,line_width=1,head_type='nofilled,
 108         head_angle=10,'xrange=xrange,'yrange=yrange,
 109         'zrange=zrange,xlabel=string(X[1]),ylabel=string(Y[1]),zlabel=string(Z[1]), args, v ) )$
 110
 111
 112
 113
 114 /*****************/
 115 /* Venn diagrams */
 116 /*       by      */
 117 /* Pankaj Sejwal */
 118 /*      2015     */
 119 /*****************/
 120
 121
 122 load(basic)$
 123
 124 fun(n):=/* To plot n-circles at equal distance from each other, find coordinates using roots of complex function of order n */
 125 map(lambda([s],[realpart(s),imagpart(s)]),makelist(cos(2*k*%pi/n)+%i*sin(2*k*%pi/n),k,1,n))$
 126
 127
 128 liss(n):=/* Create equations for circles using the coordinates obtained from fun()) */
 129 block([pts:fun(n)],pts:map(lambda([w],(x-first(w))^2+(y-last(w))^2=n),pts))$
 130
 131
 132 collectatoms(rel):=/* Collect all atoms from the relation provided to be plotted,
 133 eg, (a and b or not(c))=>[a,b,c], calls findatoms() to get job done */
 134 block([final:[]],findatoms(rel),flatten(reverse(final)))$
 135
 136
 137 findatoms(rel):=block([temp,ntemp],ntemp:args(rel),/* Collects atoms in logical relation iteratively  */
 138 temp:sublist(ntemp,lambda([s],atom(s))),
 139 if(temp#[]) then push(temp,final),
 140 for item in temp do ntemp:delete(item,ntemp),
 141 map(findatoms,ntemp))$
 142
 143
 144 transform(eq,args,lis):=/* Substitutes the equations of circles into logical relation,
 145 maxima automatically handles it as (x^2+y^2<const) in case of inclusion and not(x^2+y^2<const)
 146 equal to (x^2+y^2>const) in case of exclusion */
 147 block([temp],
 148 eq:map(lambda([s],lhs(s)<rhs(s)),eq),
 149 temp:map("=",lis,eq),
 150 args:psubst(temp,args),args)$
 151
 152
 153 randomcolor():=/*create a random color for each circle plotted */
 154 block([temp:[],final:[]],for i:1 thru 6 do
 155 (temp:[],for j:1 thru 4 do push(random(2),temp),push(temp,final)),
 156 final:psubst([[0,0,0,0]=0,[0,0,0,1]=1,[0,0,1,0]=2,[0,0,1,1]=3,[0,1,0,0]=4,
 157 [0,1,0,1]=5,[0,1,1,0]=6,[0,1,1,1]=7,[1,0,0,0]=8,[1,0,0,1]=9,
 158 [1,0,1,0]=a,[1,0,1,1]=b,[1,1,0,0]=c,[1,1,0,1]=d,[1,1,1,0]=e,[1,1,1,1]=f],final),
 159 final)$
 160
 161
 162 vennplot(args):=/* Does plotting work for the logical relation and is the only function needed by user */
 163 block([pts,reg,form,temp,wee,colr,i:1],
 164 lis:collectatoms(args),
 165 form:liss(length(lis)),
 166 n:length(lis),
 167 pts:(fun(n)),
 168 temp:transform(form,args,lis),
 169 wee:[title=string(args),proportional_axes=xy,
 170 grid= true,line_type= solid,x_voxel = 50,y_voxel = 50],
 171 reg:apply(lambda([s],region(s,x,-(n+1),n+1,y,-(n+1),n+1)),[transform(form,args,lis)]),
 172 wee:endcons(reg,wee),
 173 wee:endcons(grid=false,wee),
 174 wee:endcons(font="Courier-Oblique",wee),
 175 wee:endcons(font_size=15,wee),
 176 wee:endcons(line_width=3,wee),
 177 for item in form do
 178       ( colr:apply(concat,cons("#",randomcolor())),
 179       wee:endcons(key=string(part(lis,i)),wee),
 180       wee:endcons(color= (colr),wee),
 181       wee:endcons(label([string(part(lis,i)),first(part(pts,i)),last(part(pts,i))-0.4]),wee),i:i+1,
 182       wee:endcons(implicit(item,x,-(n+1),n+1,y,-(n+1),n+1),wee)),
 183 apply(draw2d,wee))$
 184
 185
 186 /* Usage examples:
 187
 188 vennplot(a and b and not(c))$
 189 vennplot(not(d) and b);
 190
 191 */
 192
 193