asymptote/Mechanics.asy

   1 /* Useful functions for drawing Physics 101 figures.
   2  *
   3  * Copyright (C) 2008-2011 W. Trevor King <wking@drexel.edu>
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 2 of the License, or
   8  * (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License along
  16  * with this program; if not, write to the Free Software Foundation, Inc.,
  17  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  18  */
  19
  20 import geometry;
  21 import three;
  22
  23 // ----------------- Labeled circle --------------------
  24
  25 void label_path(picture pic=currentpicture, Label L, path g, real margin=0,
  26                 pair rdir=0) {
  27   align a = L.align;
  28   embed e = L.embed;
  29   real m = labelmargin(L.p);
  30   real scale = (m + margin)/m;
  31   if (L.align.is3D) {
  32     L.align.dir3 *= scale;
  33   } else {
  34     L.align.dir *= scale;
  35   }
  36   if (L.embed == Rotate) {
  37     L.embed = Rotate(rdir);
  38   }
  39   label(pic=pic, L=L, g=g);
  40   L.align = a;
  41   L.embed = e;
  42 }
  43
  44 void label_path(picture pic=currentpicture, Label L, path3 g, real margin=0,
  45                 pair rdir=0) {
  46   align a = L.align;
  47   embed e = L.embed;
  48   real m = labelmargin(L.p);
  49   real scale = (m + margin)/m;
  50   if (L.align.is3D) {
  51     L.align.dir3 *= scale;
  52   } else {
  53     L.align.dir *= scale;
  54   }
  55   if (L.embed == Rotate) {
  56     L.embed = Rotate(rdir);
  57   }
  58   label(pic=pic, L=L, g=g);
  59   L.align = a;
  60   L.embed = e;
  61 }
  62
  63 struct LabeledCircle {
  64   pair center;
  65   real radius;
  66   pen outline;
  67   pen fill;
  68   Label label;
  69
  70   void operator init(pair center=(0,0), real radius=2mm,
  71                      pen outline=currentpen, pen fill=grey, Label L="") {
  72     this.center = center;
  73     this.radius = radius;
  74     this.outline = outline;
  75     this.fill = fill;
  76     this.label = L;
  77   }
  78
  79   void draw_label(picture pic=currentpicture, Label L=null) {
  80     align a;
  81     if (L == null) {
  82       L = this.label;
  83     }
  84     a = L.align;
  85     if (L.align != NoAlign && L.align != Align) {
  86       real m = labelmargin(L.p);
  87       real scale = (m + this.radius)/m;
  88       if (L.align.is3D) {
  89         L.align.dir3 *= scale;
  90       } else {
  91         L.align.dir *= scale;
  92       }
  93     }
  94     label(pic=pic, L=L, position=this.center);
  95     L.align = a;
  96   }
  97
  98   void draw(picture pic=currentpicture) {
  99     path p = shift(this.center)*scale(this.radius)*unitcircle;
 100     filldraw(pic, p, this.fill, this.outline);
 101     this.draw_label(pic=pic);
 102   }
 103 }
 104
 105 // ---------------------- Mass -------------------------
 106
 107 struct Mass {
 108   LabeledCircle lc;
 109   real m;
 110
 111   void operator init(pair center=(0,0), real m=1, real radius=2mm,
 112                      pen outline=currentpen, pen fill=grey, Label L="") {
 113     this.lc.operator init(center=center, radius=radius, outline=outline,
 114                           fill=fill, L=L);
 115     this.m = m;
 116   }
 117
 118   pair center() { return this.lc.center; }
 119   void set_center(pair center) { this.lc.center = center; }
 120   void draw(picture pic=currentpicture) = this.lc.draw;
 121 }
 122
 123 struct Block {
 124   pair center;
 125   real m;
 126   real width;
 127   real height;
 128   real direction;
 129   pen outline;
 130   pen fill;
 131   Label L;
 132
 133   void operator init(pair center=(0,0), real m=1, real width=5mm, real height=-1, real direction=0, pen outline=currentpen, pen fill=grey, Label L="") {
 134     this.center = center;
 135     this.m = m;
 136     this.width = width;
 137     if (height == -1)
 138       this.height = width;
 139     else
 140       this.height = height;
 141     this.direction = direction;
 142     this.outline = outline;
 143     this.fill = fill;
 144     this.L = L;
 145   }
 146
 147   void draw(picture pic=currentpicture) {
 148     picture picF;
 149     path c = rotate(direction)*scale(width,height)*shift(-.5,-.5)*unitsquare;
 150     filldraw(picF, c, fill, outline);
 151     label(pic=picF, L=L, position=(0,0));
 152     add(pic, picF, center);
 153   }
 154 }
 155
 156 // ---------------------- Vectors -------------------------
 157
 158 struct Vector {
 159   pair center;
 160   real mag;
 161   real dir; // angle in the plane of the drawing.
 162   real phi; // angle with the plane of the drawing, 90 is out of the page.
 163   pen outline;
 164   Label label;
 165   real out_of_plane_radius;
 166   real out_of_plane_tolerance;
 167
 168   void operator init(pair center=(0,0), real mag=5mm, real dir=0, real phi=0, pen outline=currentpen, Label L="") {
 169     this.center = center;
 170     this.mag = mag;
 171     this.dir = dir;
 172     this.phi = phi;
 173     this.outline = outline;
 174     this.label = L;
 175     this.out_of_plane_radius = 1mm;
 176     this.out_of_plane_tolerance = 0.01;
 177   }
 178
 179   Vector copy() {
 180     Vector v = Vector(center=this.center, mag=this.mag, dir=this.dir,
 181                       phi=this.phi, outline=this.outline, L=this.label);
 182     v.out_of_plane_radius = this.out_of_plane_radius;
 183     v.out_of_plane_tolerance = this.out_of_plane_tolerance;
 184     return v;
 185   }
 186
 187   pair dTip() {  // offset from center to tip
 188     pair p = (0,0);
 189     real phi_e = this.phi % 360; // effective phi
 190     if (Tan(phi_e) == 0 || abs(1.0/Tan(phi_e)) > this.out_of_plane_tolerance) {
 191       return this.mag*Cos(this.phi)*dir(this.dir);
 192     }
 193     return (0, 0);
 194   }
 195
 196   pair pTip() {
 197     return this.dTip() + this.center;
 198   }
 199
 200   void draw(picture pic=currentpicture) {
 201     pair p = this.dTip();
 202     path P;
 203     real phi_e = this.phi % 360; // effective phi
 204     if (this.mag < 0) (phi_e + 180) % 360;
 205     if (Tan(phi_e) == 0 || abs(1.0/Tan(phi_e)) > this.out_of_plane_tolerance) {
 206       // draw arrow in the plane of the drawing
 207       // TODO: thickening for phi?
 208       P = shift(this.center)*((0,0)--p);
 209       draw(pic, P, this.outline, Arrow);
 210     } else if (phi_e > 0 && phi_e < 180) {
 211       // draw a circled dot for out-of-the-page
 212       P = shift(this.center)*scale(this.out_of_plane_radius)*unitcircle;
 213       draw(pic, P, outline);
 214       dot(pic, this.center, this.outline);
 215     } else {
 216       // draw a circled cross for into-the-page
 217       real a = 0.8*sqrt(2.0)/2.0;
 218       P = shift(this.center)*scale(this.out_of_plane_radius)*unitcircle;
 219       draw(pic, P, this.outline);
 220       draw(pic, shift(this.center)*scale(this.out_of_plane_radius
 221                                          )*((-a,-a)--(a,a)), this.outline);
 222       draw(pic, shift(this.center)*scale(this.out_of_plane_radius
 223                                          )*((-a,a)--(a,-a)), this.outline);
 224     }
 225     label_path(pic=pic, L=this.label, g=P);
 226   }
 227 }
 228
 229 Vector operator +(Vector a, Vector b) {
 230   Vector c = a.copy();
 231   pair p = a.mag*dir(a.dir) + b.mag*dir(b.dir);
 232   c.mag = length(p);
 233   c.dir = degrees(p);
 234   return c;
 235 }
 236
 237 Vector operator -(Vector a, Vector b) {
 238   Vector c = a.copy();
 239   pair p = a.mag*dir(a.dir) - b.mag*dir(b.dir);
 240   c.mag = length(p);
 241   c.dir = degrees(p);
 242   return c;
 243 }
 244
 245 void vector_field(pair center=(0,0), real width=2cm, real height=2cm,
 246                   real dv=0.5cm, real buf=2pt, Vector v=null,
 247                   pen outline=invisible) {
 248   /* There will be a buffer of at least buf on each side */
 249   if (v == null) {
 250     v = Vector();  // unlikely to be what they want, but it will draw something
 251   }
 252
 253   pair ovcenter = v.center;
 254   real ovmag = v.mag;
 255   path bufsq = shift(center)*xscale(width-2*buf)*yscale(height-2*buf)
 256     *shift((-.5,-.5))*unitsquare;  // buffered bounding box
 257   pair uv = dir(v.dir);  // unit vector in the direction of v
 258   pair d = dv * dir(v.dir+90);
 259   real dx = d.x;
 260   real dy = d.y;
 261   int nx = 1;  // x steps
 262   int ny = 1;  // y steps
 263   bool diag = false;
 264
 265   if (abs(fmod(v.phi, 180)) == 90) {  // pure in/out, make a 2D grid
 266     dx = dy = dv;  // v.dir was meaninless, reset dx and dy
 267     nx = abs((int)((width-2*buf) / dx)) + 1;
 268     ny = abs((int)((height-2*buf) / dy)) + 1;
 269   } else if (abs(fmod(v.dir, 180)) == 0) {  // pure left/right, vert. border
 270     ny = abs((int)((height-2*buf) / dy)) + 1;
 271     dx = 0;
 272   } else if (abs(fmod(v.dir, 180)) == 90) {  // pure up/down, horiz. border
 273     nx = abs((int)((width-2*buf) / dx)) + 1;
 274     dy = 0;
 275   } else {  // diagonal, draw along a vertical an horizontal border
 276     diag = true;
 277     // this requires enough special handling that we break it out below
 278   }
 279
 280   if (!diag) {  // square grid
 281     real xx=buf, xy=buf;  // buffer distace per side
 282     if (dx != 0)
 283       xx = (width-(nx-1)*fabs(dx))/2.0;  // "extra" left over after division
 284     if (dy != 0)
 285       xy = (height-(ny-1)*fabs(dy))/2.0;
 286
 287     real xstart = center.x - width/2 + xx;
 288     real ystart = center.y - height/2 + xy;
 289     if (dx < 0 || (dx == 0 && dot(uv, dir(0)) < 0))
 290       xstart += width - 2*xx;
 291     if (dy < 0 || (dy == 0 && dot(uv, dir(90)) < 0))
 292       ystart += height - 2*xy;
 293
 294     for (int i=0; i<nx; i+=1) {
 295       for (int j=0; j<ny; j+=1)  {
 296         v.center = (xstart+i*dx, ystart+j*dy);
 297         if (abs(fmod(v.phi, 180)) != 90) {  // pure left/right/up/down
 298           // set length so that v just touches far bufsq face
 299           path vlong = (v.center+uv*min(width, height)/2)
 300             --(v.center+uv*max(width, height));
 301           //dot(v.center); draw(vlong); continue;  // intersection debugging
 302           pair xs[] = intersectionpoints(vlong, bufsq);
 303           assert(xs.length == 1, format("%d", xs.length));
 304           v.mag = length(xs[0] - v.center);
 305         }
 306         v.draw();
 307       }
 308     }
 309   } else {  // v is diagonal
 310     pair cc = (sgn(d.x)*width, sgn(d.y)*height);   // catty-corner vector
 311     real ccbuf = buf*max(fabs(1/cos(angle(cc))),   // buf away in y
 312                          fabs(1/sin(angle(cc))));  // buf away in x
 313     int n = abs((int)((dot(cc, unit(d))-2*ccbuf)/length(d))) + 1;  // lines
 314     pair pstart = center - (n-1)*d/2.0;
 315
 316     for (int i=0; i<n; i+=1) {
 317       // project along +/- v until you hit the edge of bufsq
 318       pair p = (pstart + i*d);
 319       path vlong = (p-uv*length(cc))--(p+uv*length(cc));  // extends out of box
 320       //dot(p); draw(vlong); continue;  // intersection debugging
 321       pair xs[] = intersectionpoints(vlong, bufsq);
 322       assert(xs.length == 2, format("%d", xs.length));
 323       pair start = xs[0];
 324       if (dot((start - xs[1]), uv) > 0)
 325         start = xs[1];
 326       v.center = start;
 327       v.mag = length(xs[1]-xs[0]);
 328       v.draw();
 329     }
 330   }
 331
 332   v.center = ovcenter;  // restore original center
 333   v.mag = ovmag;  // restore original magnitude
 334
 335   // draw bounding box
 336   draw(shift(center)*xscale(width)*yscale(height)*shift((-.5,-.5))*
 337        unitsquare, outline);
 338 }
 339
 340 Vector Velocity(pair center=(0,0), real mag=5mm, real dir=0, real phi=0, Label L="")
 341 {
 342   Vector v = Vector(center=center, mag=mag, dir=dir, phi=phi, L=L, outline=rgb(1,0.1,0.2)); // red
 343   return v;
 344 }
 345
 346 // ---------------------- Forces -------------------------
 347
 348 Vector Force(pair center=(0,0), real mag=5mm, real dir=0, real phi=0, Label L="")
 349 {
 350   Vector v = Vector(center=center, mag=mag, dir=dir, phi=phi, L=L, outline=rgb(0.1,0.2,1)); // blue
 351   return v;
 352 }
 353
 354 // ---------------------- Measures -------------------------
 355
 356 // Distance derived from CAD.MeasuredLine
 357 struct Distance {
 358   pair pFrom;
 359   pair pTo;
 360   real offset;
 361   real scale;
 362   pen outline;
 363   Label label;
 364
 365   void operator init(pair pFrom=(0,0), pair pTo=(5mm,0), real offset=0, real scale=5mm, pen outline=currentpen, Label L="") {
 366     this.pFrom = pFrom;
 367     this.pTo = pTo;
 368     this.offset = offset;
 369     this.scale = scale;
 370     this.outline = outline;
 371     this.label = L;
 372   }
 373
 374   void draw(picture pic=currentpicture) {
 375     pair o = this.offset*unit(rotate(-90)*(this.pTo - this.pFrom));
 376     path p = (this.pFrom+o) -- (this.pTo+o);
 377     draw(pic, p, outline, Arrows);
 378     label_path(pic=pic, L=this.label, g=p, rdir=this.pTo - this.pFrom);
 379   }
 380 }
 381
 382 struct Angle {
 383   pair B;
 384   pair A; // center of angle
 385   pair C;
 386   real radius; // radius < 0 for exterior angles.
 387   pen outline;
 388   pen fill;
 389   Label label;
 390
 391   void operator init(pair B, pair A, pair C, real radius=5mm, pen outline=currentpen, pen fill=invisible, Label L="") {
 392     this.B = B;
 393     this.A = A;
 394     this.C = C;
 395     this.radius = radius;
 396     this.outline = outline;
 397     this.fill = fill;
 398     this.label = L;
 399   }
 400
 401   void draw(picture pic=currentpicture) {
 402     bool direction;
 403
 404     real ccw_angle = degrees(C-A)-degrees(B-A);
 405     bool direction = CCW;
 406     if (ccw_angle < 0) ccw_angle += 360.0;
 407     if (ccw_angle > 180)
 408       direction = CW;
 409     if (radius < 0)
 410       direction = !direction;
 411     path p = arc(this.A, fabs(radius), degrees(B-A), degrees(C-A), direction);
 412     if (this.fill != invisible) {
 413       path pcycle = this.A -- p -- cycle;
 414       fill(pic, pcycle, this.fill);
 415     }
 416     draw(pic, p, this.outline);
 417     if (direction == CW) {
 418       p = reverse(p);
 419     }
 420     label_path(pic=pic, L=this.label, g=p);
 421   }
 422 }
 423
 424 Vector hatVect (string name,  pair center=(0,0), real dir=0, real phi=0) {
 425   string s = replace("$\mathbf{\hat{X}}$", "X", name);
 426   Label L = Label(s, position=EndPoint, align=RightSide);
 427   Vector v = Vector(
 428       center=center, mag=5mm, dir=dir, phi=phi, L=L, outline=rgb(0,0,0));
 429   return v;
 430 }
 431
 432 Vector ihat (pair center=(0,0), real dir=0, real phi=0) {
 433   Vector v = hatVect(name="i", center=center, dir=dir, phi=phi);
 434   return v;
 435 }
 436
 437 Vector jhat (pair center=(0,0), real dir=90, real phi=0) {
 438   Vector v = hatVect(name="j", center=center, dir=dir, phi=phi);
 439   return v;
 440 }
 441
 442 void draw_ijhat(pair center=(0,0), real idir=0) {
 443   Vector ihat = ihat(center, idir);
 444   Vector jhat = jhat(center, idir+90);
 445   ihat.draw();
 446   jhat.draw();
 447 }
 448
 449 // ---------------------- Shapes -------------------------
 450
 451 struct Wire {
 452   pair pFrom;
 453   pair pTo;
 454   pen outline;
 455   Label label;
 456
 457   void operator init(pair pFrom=(0,0), pair pTo=(5mm,0), pen outline=currentpen, Label L="") {
 458     this.pFrom = pFrom;
 459     this.pTo = pTo;
 460     this.outline = outline;
 461     this.label = L;
 462   }
 463
 464   void draw(picture pic=currentpicture) {
 465     path p = this.pFrom--this.pTo;
 466     draw(pic, p, outline);
 467     label_path(pic=pic, L=this.label, g=p, rdir=this.pTo - this.pFrom);
 468   }
 469 }
 470
 471 struct Surface {
 472   pair pFrom;
 473   pair pTo;
 474   real thickness;
 475   pen outline;
 476   pen filla;
 477   pen fillb;
 478   Label label;
 479
 480   void operator init(pair pFrom=(0,0), pair pTo=(5mm,0), real thickness=5mm, pen outline=currentpen, pen filla=rgb(.5,.5,.5), pen fillb=rgb(.7,.7,.7), Label L="") {
 481     this.pFrom = pFrom;
 482     this.pTo = pTo;
 483     this.thickness = thickness;
 484     this.outline = outline;
 485     this.filla = filla;
 486     this.fillb = fillb;
 487     this.label = L;
 488   }
 489
 490   void draw(picture pic=currentpicture) {
 491     pair pDiff = pTo - pFrom;
 492     pair pDepth = rotate(-90)*unit(pDiff)*thickness;
 493     path p = (0, 0) -- pDiff -- (pDiff + pDepth) -- pDepth -- cycle;
 494     p = shift(this.pFrom) * p;
 495     axialshade(
 496         pic=pic, g=p, pena=filla, a=this.pFrom,
 497         penb=fillb, b=this.pFrom+pDepth);
 498     draw(pic, p, outline);
 499     label_path(pic=pic, L=this.label,
 500         g=shift(pDepth/2)*(this.pFrom -- this.pTo), rdir=this.pTo - this.pFrom);
 501   }
 502 }
 503
 504 struct Spring {
 505   pair pFrom;
 506   pair pTo;
 507   real k;
 508   real width;
 509   real dLength; // length of a single loop (when unstretched)
 510   real deadLength; // length before loops start on each end
 511   real unstretchedLength;
 512   int nLoops;
 513   pen outline;
 514   Label label;
 515
 516   void operator init(pair pFrom=(0,0), pair pTo=(5mm,0), real k=1, real width=3mm, real dLength=1mm, real deadLength=3mm, real unstretchedLength=12mm, pen outline=currentpen, Label L="") {
 517     this.pFrom = pFrom;
 518     this.pTo = pTo;
 519     this.k = k;
 520     this.width = width;
 521     this.dLength = dLength;
 522     this.unstretchedLength = unstretchedLength;
 523     this.outline = outline;
 524     this.label = L;
 525
 526     this.nLoops = floor((unstretchedLength-2*deadLength)/dLength);
 527     if (this.nLoops == 0)
 528       this.nLoops = 1;
 529     this.deadLength = (unstretchedLength-this.nLoops*dLength)/2;
 530   }
 531
 532   void draw(picture pic=currentpicture) {
 533     pair pDiff = pTo - pFrom;
 534
 535     real working_dLength = (length(pDiff) - 2*deadLength) / nLoops;
 536     path p = (0,0)--(deadLength,0);
 537     path loop = (0,0)
 538       --(working_dLength/4, width/2)
 539       --(working_dLength*3/4, -width/2)
 540       --(working_dLength, 0);
 541     pair loopStart;
 542     for (int i=0; i<nLoops; ++i) {
 543       loopStart = point(p, length(p));
 544       p = p--(shift(loopStart)*loop);
 545     }
 546     loopStart = point(p, length(p));
 547     p = p--(loopStart+(deadLength,0));
 548     p = shift(this.pFrom)*rotate(degrees(pDiff)) * p;
 549     draw(pic, p, outline);
 550     label_path(pic=pic, L=this.label, g=this.pFrom--this.pTo,
 551         margin=this.width/2, rdir=this.pTo - this.pFrom);
 552   }
 553 }
 554
 555 struct Pendulum {
 556   pair pivot;
 557   Mass mass;
 558   Angle angle;
 559   Wire str;
 560
 561   void operator init(pair pivot=(0,0), Mass mass=Mass()) {
 562     this.pivot = pivot;
 563     this.mass = mass;
 564     this.angle = Angle(mass.center(), pivot, pivot-(0,1));
 565     this.str = Wire(pivot, mass.center());
 566   }
 567
 568   void draw(picture pic=currentpicture, bool drawVertical=false) {
 569     str.draw(pic=pic);
 570     if (drawVertical == true) {
 571       pair final = pivot + realmult((0,0.5),(mass.center()-pivot));
 572       draw(pic=pic, pivot--final, p=currentpen+dashed);
 573     }
 574     draw(pic=pic, pivot);
 575     mass.draw(pic=pic);
 576     angle.draw(pic=pic);
 577   }
 578 }
 579
 580 // The angle argument is deflection from straight down (i.e. 0 degrees = plumb)
 581 Pendulum makePendulum(pair pivot=(0,0), Mass mass=Mass(), real length=15mm, real angleDeg=0, Label angleL="", Label stringL="") {
 582   mass.set_center(pivot + length*dir(angleDeg-90));
 583   Pendulum p = Pendulum(pivot=pivot, mass=mass);
 584   p.angle.label = angleL;
 585   p.str.label = stringL;
 586   return p;
 587 }
 588
 589 struct Ring {
 590   triple center;
 591   triple normal;
 592   real radius;
 593   real width;
 594   real axis_pre;
 595   real axis_post;
 596   pen outline;
 597   pen fill;
 598   pen axis;
 599   Label label;
 600   Label axis_label;
 601
 602   void operator init(triple center=(0,0,0), triple normal=(0,0,1),
 603                      real radius=5mm, real width=1mm, real axis_pre=0,
 604                      real axis_post=0, pen outline=currentpen,
 605                      pen fill=invisible, pen axis=invisible, Label L="",
 606                      Label axis_label="") {
 607     this.center = center;
 608     this.normal = normal;
 609     this.radius = radius;
 610     this.width = width;
 611     this.axis_pre = axis_pre;
 612     this.axis_post = axis_post;
 613     this.outline = outline;
 614     this.fill = fill;
 615     this.axis = axis;
 616     this.label = L;
 617     this.axis_label = axis_label;
 618   }
 619
 620   void draw(picture pic=currentpicture) {
 621     path3 a = shift(this.center)*(
 622         (-this.axis_pre*this.normal)--(this.axis_post*this.normal));
 623     draw(pic, a, this.axis, Arrow3);
 624     path3 c = circle(c=this.center, r=this.radius, normal=this.normal);
 625     //tube t = tube(c, width=this.width);  // slow and ugly!
 626     //draw(t.s, this.fill);
 627     draw(pic, c, this.outline);
 628     label_path(pic=pic, L=this.label, g=c);
 629     label_path(pic=pic, L=this.axis_label, g=a);
 630   }
 631 }
 632
 633 // TODO: plate, cylinder, table