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