From e3a2db7bb6a8a4b238350f6a7e71120366848ab8 Mon Sep 17 00:00:00 2001
From: "W. Trevor King" <wking@drexel.edu>
Date: Fri, 27 May 2011 15:19:03 -0400
Subject: [PATCH] Add vector_field() to Mechanics.  Also do some whitespace
 cleanups.

---
 asymptote/Mechanics-test.asy |  50 ++++++++++++-
 asymptote/Mechanics.asy      | 135 +++++++++++++++++++++++++++++------
 2 files changed, 164 insertions(+), 21 deletions(-)

diff --git a/asymptote/Mechanics-test.asy b/asymptote/Mechanics-test.asy
index acbc701..a63e7d0 100644
--- a/asymptote/Mechanics-test.asy
+++ b/asymptote/Mechanics-test.asy
@@ -50,7 +50,7 @@ for (int i=0; i<vs.length; i+=1)
   vs[i].draw();
 
 Pendulum p = makePendulum(pivot=(3.5u,-3u), mass=b, length=4u, angleDeg=-20,
-			  angleL="$\rho$", stringL="r");
+                          angleL="$\rho$", stringL="r");
 
 real len = abs(p.pivot.x-b.center.x);
 Spring s = Spring(pFrom=b.center-(2u,0), pTo=b.center, L="$k_1$");  s.draw();
@@ -60,3 +60,51 @@ p.draw(drawVertical=true);
 
 Surface s = Surface(pFrom=(0,-7u), pTo=(3.5u, -7u), L="Table");
 s.draw();
+
+
+Vector v = Velocity();
+pair vfc = (5u,0);  // vector field center
+real vfwidth = 2u;
+real vfheight = 2u;
+real vfdv = u/3;
+real vfbuf = 2pt;
+
+draw(shift(vfc)*xscale(vfwidth-2*vfbuf)*yscale(vfheight-2*vfbuf)
+     *shift((-.5,-.5))*unitsquare, grey);
+vector_field(vfc, width=vfwidth, height=vfheight, dv=vfdv, buf=vfbuf, v=v,
+             outline=green+dashed);
+
+vfc -= (0, vfheight + u);
+v.dir = 90;
+draw(shift(vfc)*xscale(vfwidth-2*vfbuf)*yscale(vfheight-2*vfbuf)
+     *shift((-.5,-.5))*unitsquare, grey);
+vector_field(vfc, width=vfwidth, height=vfheight, dv=vfdv, buf=vfbuf, v=v,
+             outline=green+dashed);
+
+vfc -= (0, vfheight + u);
+v.dir = -180;
+draw(shift(vfc)*xscale(vfwidth-2*vfbuf)*yscale(vfheight-2*vfbuf)
+     *shift((-.5,-.5))*unitsquare, grey);
+vector_field(vfc, width=vfwidth, height=vfheight, dv=vfdv, buf=vfbuf, v=v,
+             outline=green+dashed);
+
+vfc += (vfwidth + u, 2*(vfheight + u));
+v.dir = -90;
+draw(shift(vfc)*xscale(vfwidth-2*vfbuf)*yscale(vfheight-2*vfbuf)
+     *shift((-.5,-.5))*unitsquare, grey);
+vector_field(vfc, width=vfwidth, height=vfheight, dv=vfdv, buf=vfbuf, v=v,
+             outline=green+dashed);
+
+vfc -= (0, vfheight + u);
+v.dir = -10;
+draw(shift(vfc)*xscale(vfwidth-2*vfbuf)*yscale(vfheight-2*vfbuf)
+     *shift((-.5,-.5))*unitsquare, grey);
+vector_field(vfc, width=vfwidth, height=vfheight, dv=vfdv, buf=vfbuf, v=v,
+             outline=green+dashed);
+
+vfc -= (0, vfheight + u);
+v.phi = 90;
+draw(shift(vfc)*xscale(vfwidth-2*vfbuf)*yscale(vfheight-2*vfbuf)
+     *shift((-.5,-.5))*unitsquare, grey);
+vector_field(vfc, width=vfwidth, height=vfheight, dv=vfdv, buf=vfbuf, v=v,
+             outline=green+dashed);
diff --git a/asymptote/Mechanics.asy b/asymptote/Mechanics.asy
index 1753a88..2836784 100644
--- a/asymptote/Mechanics.asy
+++ b/asymptote/Mechanics.asy
@@ -121,7 +121,7 @@ struct Vector {
       path P = (0,0)--p;
       draw(picF, P, outline, Arrow);
       if (rotateLabel == true)
-	label_rotate = p;
+        label_rotate = p;
       label_align = unit(rotate(90)*p);
     } else if (phi_e > 0 && phi_e < 180) {
       // draw a circled dot for out-of-the-page
@@ -137,13 +137,108 @@ struct Vector {
       draw(picF, scale(this.out_of_plane_radius)*((-a,a)--(a,-a)), outline);
     }
     label(pic = picF,
-	  L = rotate(degrees(label_rotate)) * L,
-	  position = p+labelOffset,
-	  align = label_align);
+          L = rotate(degrees(label_rotate)) * L,
+          position = p+labelOffset,
+          align = label_align);
     add(pic, picF, center);
   }
 }
 
+void vector_field(pair center=(0,0), real width=2cm, real height=2cm,
+                  real dv=0.5cm, real buf=2pt, Vector v=null,
+                  pen outline=invisible) {
+  /* There will be a buffer of at least buf on each side */
+  if (v == null) {
+    v = Vector();  // unlikely to be what they want, but it will draw something
+  }
+
+  pair ovcenter = v.center;
+  real ovmag = v.mag;
+  path bufsq = shift(center)*xscale(width-2*buf)*yscale(height-2*buf)
+    *shift((-.5,-.5))*unitsquare;  // buffered bounding box
+  pair uv = dir(v.dir);  // unit vector in the direction of v
+  pair d = dv * dir(v.dir+90);
+  real dx = d.x;
+  real dy = d.y;
+  int nx = 1;  // x steps
+  int ny = 1;  // y steps
+  bool diag = false;
+
+  if (abs(fmod(v.phi, 180)) == 90) {  // pure in/out, make a 2D grid
+    dx = dy = dv;  // v.dir was meaninless, reset dx and dy
+    nx = abs((int)((width-2*buf) / dx)) + 1;
+    ny = abs((int)((height-2*buf) / dy)) + 1;
+  } else if (abs(fmod(v.dir, 180)) == 0) {  // pure left/right, vert. border
+    ny = abs((int)((height-2*buf) / dy)) + 1;
+    dx = 0;
+  } else if (abs(fmod(v.dir, 180)) == 90) {  // pure up/down, horiz. border
+    nx = abs((int)((width-2*buf) / dx)) + 1;
+    dy = 0;
+  } else {  // diagonal, draw along a vertical an horizontal border
+    diag = true;
+    // this requires enough special handling that we break it out below
+  }
+
+  if (!diag) {  // square grid
+    real xx=buf, xy=buf;  // buffer distace per side
+    if (dx != 0)
+      xx = (width-(nx-1)*fabs(dx))/2.0;  // "extra" left over after division
+    if (dy != 0)
+      xy = (height-(ny-1)*fabs(dy))/2.0;
+
+    real xstart = center.x - width/2 + xx;
+    real ystart = center.y - height/2 + xy;
+    if (dx < 0 || (dx == 0 && dot(uv, dir(0)) < 0))
+      xstart += width - 2*xx;
+    if (dy < 0 || (dy == 0 && dot(uv, dir(90)) < 0))
+      ystart += height - 2*xy;
+
+    for (int i=0; i<nx; i+=1) {
+      for (int j=0; j<ny; j+=1)  {
+        v.center = (xstart+i*dx, ystart+j*dy);
+        if (abs(fmod(v.phi, 180)) != 90) {  // pure left/right/up/down
+          // set length so that v just touches far bufsq face
+          path vlong = (v.center+uv*min(width, height)/2)
+            --(v.center+uv*max(width, height));
+          //dot(v.center); draw(vlong); continue;  // intersection debugging
+          pair xs[] = intersectionpoints(vlong, bufsq);
+          assert(xs.length == 1, format("%d", xs.length));
+          v.mag = length(xs[0] - v.center);
+        }
+        v.draw();
+      }
+    }
+  } else {  // v is diagonal
+    pair cc = (sgn(d.x)*width, sgn(d.y)*height);   // catty-corner vector
+    real ccbuf = buf*max(fabs(1/cos(angle(cc))),   // buf away in y
+                         fabs(1/sin(angle(cc))));  // buf away in x
+    int n = abs((int)((dot(cc, unit(d))-2*ccbuf)/length(d))) + 1;  // lines
+    pair pstart = center - (n-1)*d/2.0;
+
+    for (int i=0; i<n; i+=1) {
+      // project along +/- v until you hit the edge of bufsq
+      pair p = (pstart + i*d);
+      path vlong = (p-uv*length(cc))--(p+uv*length(cc));  // extends out of box
+      //dot(p); draw(vlong); continue;  // intersection debugging
+      pair xs[] = intersectionpoints(vlong, bufsq);
+      assert(xs.length == 2, format("%d", xs.length));
+      pair start = xs[0];
+      if (dot((start - xs[1]), uv) > 0)
+        start = xs[1];
+      v.center = start;
+      v.mag = length(xs[1]-xs[0]);
+      v.draw();
+    }
+  }
+
+  v.center = ovcenter;  // restore original center
+  v.mag = ovmag;  // restore original magnitude
+
+  // draw bounding box
+  draw(shift(center)*xscale(width)*yscale(height)*shift((-.5,-.5))*
+       unitsquare, outline);
+}
+
 Vector Velocity(pair center=(0,0), real mag=5mm, real dir=0, real phi=0, Label L="")
 {
   Vector v = Vector(center=center, mag=mag, dir=dir, phi=phi, L=L, outline=rgb(1,0.1,0.2)); // red
@@ -179,7 +274,7 @@ struct Distance {
   }
   
   void draw(picture pic=currentpicture, bool rotateLabel=true,
-	    real labelangle=90, real labeloffset=-1) {
+            real labelangle=90, real labeloffset=-1) {
     picture picF;
     pair pDiff = pTo - pFrom;
     if (labeloffset == -1) {
@@ -194,9 +289,9 @@ struct Distance {
       label_rotate=(1,0);
     draw(picF, p, outline, Arrows);
     label(pic = picF,
-	  L = rotate(degrees(label_rotate)) * L,
-	  position =  pDiff/2
-	    + unit(rotate(labelangle)*pDiff) * labeloffset);
+          L = rotate(degrees(label_rotate)) * L,
+          position =  pDiff/2
+            + unit(rotate(labelangle)*pDiff) * labeloffset);
     //label(pic=picF, L = rotate(degrees(label_rotate)) format("%g", pDiff/scale), position = TODO);
     add(pic, picF, pFrom+offset*unit(rotate(-90)*pDiff));
   }
@@ -220,7 +315,7 @@ struct Angle {
   }
   
   void draw(picture pic=currentpicture, bool rotateLabel=false,
-	    real labelOffsetAdjustment=0) {
+            real labelOffsetAdjustment=0) {
     picture picF;
     picture picL;
     bool direction;
@@ -246,8 +341,8 @@ struct Angle {
     
     draw(picF, p, outline);    
     label(pic = picF,
-	  L = rotate(degrees(label_rotate)) * L,
-	  position = P + unit(P) * (pLabelSize.y / 2 + labelOffsetAdjustment));
+          L = rotate(degrees(label_rotate)) * L,
+          position = P + unit(P) * (pLabelSize.y / 2 + labelOffsetAdjustment));
     add(pic, picF, A);
   }
 }
@@ -302,9 +397,9 @@ struct Wire {
       label_rotate=(1,0);
     draw(picF, p, outline);
     label(pic = picF,
-	  L = rotate(degrees(label_rotate)) * L,
-	  position =  pDiff/2
-	    + unit(rotate(90)*pDiff) * pLabelSize.y / 2);
+          L = rotate(degrees(label_rotate)) * L,
+          position =  pDiff/2
+            + unit(rotate(90)*pDiff) * pLabelSize.y / 2);
     add(pic, picF, pFrom);
   }
 }
@@ -340,11 +435,11 @@ struct Surface {
     if (rotateLabel == false)
       label_rotate=(1,0);
     axialshade(pic=picF, g=p, pena=filla, a=(0,0), penb=fillb,
-	       b=pDepth);
+               b=pDepth);
     draw(picF, p, outline);
     label(pic = picF,
-	  L = rotate(degrees(label_rotate)) * L,
-	  position = (pDiff+pDepth)/2);
+          L = rotate(degrees(label_rotate)) * L,
+          position = (pDiff+pDepth)/2);
     add(pic, picF, pFrom);
   }
 }
@@ -404,9 +499,9 @@ struct Spring {
       label_rotate=(1,0);
     draw(picF, p, outline);
     label(pic = picF,
-	  L = rotate(degrees(label_rotate)) * L,
-	  position =  pDiff/2
-	  + unit(rotate(90)*pDiff) * (width + pLabelSize.y) / 2);
+          L = rotate(degrees(label_rotate)) * L,
+          position =  pDiff/2
+          + unit(rotate(90)*pDiff) * (width + pLabelSize.y) / 2);
     add(pic, picF, pFrom);
   }
 }
-- 
2.26.2