Pencil pushers

A couple of days ago, Matt Parker posted this video to his Standup Maths YouTube channel.

In it, he and Hugh Hunt, reader in engineering dynamics and vibration at Cambridge, stand a pencil on its end at the edge of a table, give the bottom end of it a horizontal whack, and see what happens. They also calculate what happens… sort of.

See, the thing is, the title of the video is “Will a falling pencil hit the table? We do the maths!” They certainly do some math, but not enough to answer the question in the title. Kind of disappointing. So let’s do a little more math and finish the problem for them. In the following, we’ll be referring to their work, so you really should watch the video before proceeding.

Here’s the setup. A pencil, which we are modeling as a uniform rod, is standing upright at the edge of a table. We’re going to apply a sharp, short duration horizontal force at the bottom of the rod and work out its subsequent motion. The sketch below shows the layout and the coordinate system I want to use in the solution.

Rod setup

The [g] represents the acceleration due to gravity. The rod is characterized by its length, [2a], and its mass, [m]. Because the mass is uniformly distributed along its length, the center of mass, [C], is halfway up the rod, and we can calculate its moment of inertia about [C] as

[I_C = \frac{1}{3} m a^2]

(If you’re looking at that formula and your first thought is “Why isn’t that a 12 in the denominator?” join the club. It took me a few minutes to remember that Hugh had defined [a] as half the length of the rod, and that’s why there’s a 3 instead of a 12. We’re used to seeing the formula written in terms of the full length

[I_C = \frac{1}{12} m L^2]

but since [L = 2a],

[I_C = \frac{1}{12} m (2a)^2 = \frac{1}{12} m (4a^2) = \frac{1}{3} m a^2]

so the 3 is correct.)

[A] is the end that’s initially at the top, and we’ll be tracking its position to find out if hits the table or flies past it.

[\hat{\:F}] is the impulse with which we strike the bottom of the rod. It’s shown as acting a little above the bottom, but I drew it that way so the arrow wouldn’t get lost in the top of the table. Impulse is defined as the integral of the force over the duration of its application. We’re going to assume the duration of the impulse is short, [\Delta t], and that it ends at time [t = 0]. Thus, the impulse is

[\hat{\:F} = \int_{-\Delta t}^0 F(t)\,dt]

The impulse-momentum equation is essentially just Newton’s second law integrated over time. For this situation, in the horizontal direction, we have

[\hat{\:F} = m v_0 - m v_{-}]

where [v_0] is the velocity of C at [t = 0] (just after the impulse), and [v_{-}] is the velocity of C at [t = -\Delta t] (just before the impulse). Because the rod starts out stationary, [v_{-} = 0] and

[\hat{\:F} = m v_0]

There’s a rotational analog to the linear impulse-momentum equation

[\hat{\:M}_C = \hat{\:F} a = I_C \,\omega_0 - I_C \,\omega_{-}]

where [\omega_0] is the angular velocity of the rod at [t = 0] and [\omega_{-}] is the angular velocity at [t = - \Delta t]. Because the rod is stationary before we hit it, [\omega_{-} = 0], so

[\hat{\:F} a = I_C \omega_0 = \frac{1}{3} m a^2 \omega_0]

There’s some simplification here that the video doesn’t mention. The full definition of angular impulse about [C] is

[\hat{\:M}_C = \int_{- \Delta t}^0 F(t) \, r_{\perp}(t) \, dt]

where [r_{\perp}(t)] is the perpendicular distance from the line of action of [\hat{\:F}] to [C]. It’s because the impulse is assumed to act over a very short length of time that we can take [r_{\perp}(t) = a] throughout the impulse. The idea is that impulse is so short, the rod doesn’t have a chance to rotate significantly. And because [a] is not a function of [t], we can bring it out of the integral, leaving us with [\hat{\:M}_C = \hat{\:F} a].

Combining our two impulse-momentum equations to eliminate [\hat{\:F}] gives us

[v_0 = \frac{1}{3} a \omega_0]

which is exactly what the video got.

Now let’s track the motions of [C] and [A]. After the impulse, the bottom of the rod will be just off the edge of the table, and the only force acting on the rod will be gravity (in the grand tradition of all elementary mechanics problems, we’re ignoring air resistance). So [C] will trace out a parabola, with

[x_C = v_0 \,t = \frac{1}{3} a \,\omega_0\, t]


[y_C = a - \frac{1}{2} g \,t^2]

being the parabola’s parametric equations.

With the only force acting through the rod’s center of mass, the angular velocity will remain constant at [\omega_0]. So

[\theta = \omega_0 \,t]

and we can use that to write the equations for the coordinates of [A]:

[x_A = x_C - a \sin \omega_0 t = \frac{1}{3} a\, \omega_0\, t - a \sin \omega_0 t] [y_A = y_C + a \cos \omega_0 t = a\,(1 + \cos \omega_0 t) - \frac{1}{2} g\, t^2]

Now, taking a hint from the video, we’re going to consider the conditions at which the tip of the rod, [A], is even, horizontally, with the edge of the table. In other words, when [x_A = 0].

[x_A = \frac{1}{3} a\, \omega_0\, t - a \sin \omega_0 t = 0]

That means

[\frac{1}{3} a\, \omega_0\, t = a \sin \omega_0 t]


[\frac{1}{3} \theta = \sin \theta]

This is the formula they got in the video, but they were looking only at the special case for which the tip of the pencil just barely hits the corner of the table. What we’ve shown is that the same formula works whenever the tip of the pencil passes a vertical line even with the edge of the table. And instead of stopping here, as the video does, we can use the solution to this formula to go further and answer the question in the video’s title: Under what conditions will the pencil hit the table?

The solution of the formula, to an unnecessarily precise seven digits, is [\theta = 2.278863], which is equivalent to about 130.57°. What we need to do is consider the value of [y_A] when [\theta] is equal to this value. There are three conditions of interest:

  1. If [y_A > 0] when [\theta = 2.278863], then the tip of the pencil will pass above the table and there will be no hit.
  2. If [y_A = 0] when [\theta = 2.278863], then the tip of the pencil will just graze the corner of the table. This is the boundary between hit and miss and is the condition they focused on in the video.
  3. If [y_A < 0] when [\theta = 2.278863], then the pencil hits the table.

Rod positions

Yes, [y_A] being less than zero for [\theta = 2.278863] is physically impossible because the pencil will hit the table before [\theta] can get to that value. But that’s OK. All we really need to know is the condition for which the pencil will hit the table. We can’t expect to track the motion of the pencil after it bounces off the table.

So let’s do some algebra on our expression for [y_A]. We can use the [\theta = \omega_0 t] relationship to eliminate [t] from the right-hand side:

[y_A = a\, (1 + \cos \theta) - \frac{1}{2}\frac{g}{\omega_0^2} \theta^2]

Dividing through by [a], we get a formula with nondimensional terms, which will allow us to do some numerical work.

[\frac{y_A}{a} = (1 + \cos \theta) - \frac{1}{2}\frac{g}{a\, \omega_0^2} \theta^2]

Substituting in [\theta = 2.278863] and rearranging, we get these three conditions:

  1. [y_A > 0] when [\frac{g}{a\, \omega_0^2} < 0.134650] (misses table)
  2. [y_A = 0] when [\frac{g}{a\, \omega_0^2} = 0.134650] (grazes table)
  3. [y_A < 0] when [\frac{g}{a\, \omega_0^2} > 0.134650] (hits table)

So for a given rod length and a given gravitational field, it’s the angular velocity imparted to the rod by the impulse that determines whether the rod strikes the table or not. If the angular velocity is above this limit,

[\omega_0 = \sqrt{\frac{g}{0.134650\, a}}]

the rod will not hit the table. If the angular velocity is below this value, it will.

We can now work out what this means in terms of the magnitude of the impulse. Recall that

[\hat{\:F} = \frac{1}{3} m\, a\, \omega_0]

so we can solve for the value of [\hat{\:F}] that is the boundary between hit and miss:

[\hat{\:F} = 0.908396\, m\, \sqrt{g a}]

Impulses higher than this will be misses, impulses lower will be hits.

If you stop for a minute and think about this result, you should be able to convince yourself that it makes sense. Clearly, a larger mass will require a stronger impulse to get the rod to fly off the table. And a higher gravitational field would mean the rod drops more quickly, also requiring a stronger impulse to get it to miss the table. The relationship between [\hat{\:F}] and [a] is perhaps not as obvious, but think of it this way: The longer the rod, the more slowly it rotates for a given impulse; therefore, the only way to get it to swing around to the required angle of [\theta = 2.278863] before striking the table is to hit it harder.

Now let’s run some numbers for what we see in the video. We’ll assume the pencil is of length 18 cm. With [g = 9.81\:\mathrm{m/s^2}] and [a = 0.09\:\mathrm{m}], we get [\omega_0 = 28.45\:\mathrm{rad/s}] (about 270 rpm) as the boundary between hit and miss. That corresponds to

[v_0 = \frac{1}{3}a\, \omega_0 = 0.85\:\mathrm{m/s}]

as the horizontal component of velocity at [C]. The velocity of the base of the pencil just after the impulse needs to be

[v_0 + a\, \omega_0 = \frac{4}{3} a\, \omega_0 = 3.4\:\mathrm{m/s}]

These values seem pretty easy to achieve, which is probably why most of the hits in the video sent their pencils flying off without hitting the table.

Larry Tesler

You probably read a few articles this past week commemorating Larry Tesler, the father of cut/copy-paste and the most prominent advocate of modeless computing. It’s good to read what other have to say about him, but if you have the time you should read some of his own writing.1

His ACM article, “A Personal History of Modeless Text Editing and Cut/Copy-Paste,” is probably the best distillation of the ideas he’s known for and the connection between them. The main thread is a fun historical narrative, but don’t skip over the sidebars in the article. In “How Modes Degrade Usability,” he ties modes to the once-predominant verb/object method of telling a computer what to do. In verb/object computing, you give the computer an action command and follow that with the item(s) on which to perform the action. On a computer of today, that would be like selecting Copy from the Edit menu and then selecting the text to copy.

That may sound insane, but one of the things that made the Mac initially puzzling to people who were heavy users of PCs and mainframes back in the mid-80s was that it flipped verb/object around. I remember having to explain to several people that the main thing to remember on a Mac was that you never go to the menu bar for a command until you’ve selected the thing that the command would act on. The verb/object ordering was hard to shake for them because it seemed very English-like and conversational. Those of us who use the command line (or Vim2) still do a lot of verb/object computing, but nowadays we do it within a modeless windowing environment. Modal computing tends to be the exception rather than the rule.

Tesler’s ideas about how computing should be done has become so pervasive I imagine it’s almost inconceivable to younger readers that there was ever a different way. When Steve Jobs stepped down as CEO of Apple, I wrote “In a sense, you’re using a Steve Jobs product whether it has an Apple logo or not.” Much the same can be said of Larry Tesler.

  1. I really like his annotated version of the PUB manual, but I have a thing for markup languages. 

  2. Kakoune is a modal text editor that turns verb/object around. By highlighting selections (like Vim’s visual mode) all the time, it allows you to make selections first, see what they are, and then issue the commands that act on them. It’s an interesting exception to the connection Tesler draws between modal and verb/object ordering. 

A name display bug in Files

I’ve written in the past about how I’ve had to use tags to work around a misfeature of Files: that it doesn’t display full file names, including the extensions. Today I learned that there’s a bug in that misfeature.

I was working from my iPad on some files that were stored in iCloud Drive. Here’s Files showing the contents of the folder I was working in:

Files inconsistent in showing extensions

The file with the blue dot has a .tex extension. The extension isn’t shown, which is the normal (and frustrating) way Files behaves. But the file with the green dot has a .py extension and its extension is shown. What’s going on?

My first concern was that the green file somehow had two extensions—that its full name was—so I opened a Terminal in Textastic,1 navigated to the appropriate iCloud folder, and ran ls to see the real, full names of both files.

Terminal session showing full file names

As you can see, there was no double extension.2 The problem wasn’t that I had somehow misnamed the Python file, it was that Files has a bug in the name display. It happens to be a bug I approve of, and one that I wish were consistent across all file types at all times, but it’s a bug nonetheless.

It says something about the state of software quality at Apple that it can’t even get its design errors right.

  1. Normally, I’d do this in Prompt, but I was in the middle of editing one of the files in Textastic when I noticed the extension anomaly in Files, so it was quicker to start a terminal session there. I really like Textastic. 

  2. If you’re looking at the path in the terminal session and wondering why an iCloud folder isn’t in ~/Library/Mobile Documents/com~apple~CloudDocs, its because I made a symbolic link to the projects folder in iCloud from a projects directory in my home directory. Saves time when using the cd command, whether I’m working directly on my Mac or remotely through terminal emulation. 

Floor plans with Python and Shapely

In yesterday’s post, we outlined portions of a bitmapped floor plan drawing of the the US Capitol, saved the resulting graphic as an SVG file, and used one of Python’s XML libraries to extract the coordinates of the various boundaries. We ended up with seven CSV files, outline.csv, scale.scv, and five files named cutout-01.csv through cutout-05.csv. Today, we’ll use the Shapely library to calculate some properties of the shapes we drew.

Annotated Capitol floor plan

Recall that the blue lines are the overall outline of the building’s third floor interior, the red lines are the outlines of the openings in the third floor that accommodate various multistory rooms in the building, and the green box is set to match the length of the 64-foot scale below the title.

The feature of the Shapely library that does the lion’s share of the calculating for us is the Polygon object, which has some very nice properties for handling real-world shapes. First, it doesn’t require the polygon to be convex, so we can model the reentrant corners that are common in buildings, landscaping, and property boundaries. Second, it allows for holes within the outer shell of the polygon, which is perfect for handling floor openings, atriums, and other items that make shapes multiply connected.

Our goal is to work out the floor area of the third floor. That’s the area enclosed by the outer walls minus the openings. Here’s the code that does it:

 1:  from shapely.geometry import Polygon
 2:  import pandas as pd
 4:  # Work out the scale of the drawing
 5:  dfscale = pd.read_csv('scale.csv')
 6:  scale = 64/(dfscale.x.max() - dfscale.x.min())
 8:  # Get the outline and scale it
 9:  dfoutline = pd.read_csv('outline.csv')
10:  dfoutline.x = dfoutline.x*scale
11:  dfoutline.y = dfoutline.y*scale
13:  # Do the same for the cutouts
14:  dfcutouts = []
15:  for i in range(5):
16:    dfcutouts.append(pd.read_csv(f'cutout-{i+1:02}.csv'))
17:    dfcutouts[i].x = dfcutouts[i].x*scale
18:    dfcutouts[i].y = dfcutouts[i].y*scale
20:  # Make a polygon without holes
21:  ocoords = list(zip(dfoutline.x, dfoutline.y))
22:  outline = Polygon(ocoords)
23:  print(f'{"Area of outline":>20s}: {outline.area:7,.0f} ft²')
25:  # Make a polygon with holes to represent the floor
26:  ccoords = []
27:  for i in range(5):
28:    ccoords.append(list(zip(dfcutouts[i].x, dfcutouts[i].y)))
29:  floor = Polygon(ocoords, ccoords)
30:  print(f'{"Floor area":>20s}: {floor.area:7,.0f} ft²')
32:  # Check all the holes individually
33:  holes = ['House Chamber', 'Statuary Hall', 'Great Rotunda', 'Old Senate Chamber', 'Senate Chamber']
34:  for i in range(5):
35:    hole = Polygon(ccoords[i])
36:    print(f'{holes[i]:>20s}: {hole.area:7,.0f} ft²')

And here are the results:

     Area of outline: 134,107 ft²
          Floor area:  95,848 ft²
       House Chamber:  13,493 ft²
       Statuary Hall:   5,128 ft²
       Great Rotunda:   7,010 ft²
  Old Senate Chamber:   2,869 ft²
      Senate Chamber:   9,759 ft²

The script starts by using the Pandas library, to read in the CSV files and arrange them in to dataframes. Pandas isn’t really necessary, as the Python standard library has a CSV module, but I’m used to Pandas and I like how quickly you can do calculations on dataframes.

Lines 5 and 6 use the coordinates of the green box saved in scale.csv to calculate the scale of the drawing in feet per pixel. That value, saved in the variable scale, is later used to convert the coordinates of the outline and the cutouts from pixels to feet.

When Lines 9–11 are done, the dfoutline dataframe contains the coordinates in feet of all the vertices of the outer boundary. Similarly, when Lines 14–18 are done, dfcutouts is a list of dataframes with coordinates in feet, one for each of the floor openings.

Now we start using Shapely. Line 21 puts the x and y values of dfoutline into a list of coordinate pairs. Line 22 constructs a Polygon, named outline, from that list, and Line 23 uses the area property to print out the area. This is the gross area enclosed by the boundary. Our next step is to figure out the net area.

Lines 26–28 set up a list of lists of coordinate pairs for all the floor openings. The floor is then defined in Line 29 using a variant on the Polygon constructor that includes holes. This means floor is a multiply connected polygon, and the area property used in Line 30 accounts for the holes.

Note that if all we wanted was the net floor area, we didn’t have to create the outline polygon. I just did that so we could see the difference.

Another unnecessary bit is the set of calculations in Lines 33–36, which tells us the areas of each of the openings. It just lets us check that the sum of the holes is equal to the difference between the outline and floor areas.

While Shapely would be nice if all it did was calculate the areas of irregular shapes, that isn’t all it does. One of the things I use it for is to generate random sampling spots for floor tile, grout, slab concrete, etc. I can generate (x, y) coordinates over the enclosing rectangle of a floor and then use the contains method to filter out those that aren’t within the irregular floor shape of interest. There are also some very helpful set-theoretic functions like intersection, union, and difference, for handling relationships between several shapes.

I should mention that breaking the work into two scripts, one for creating the vertex CSVs and the other for calculating the area, is inefficient. A single script that reads the SVG and puts the vertex coordinates directly into lists would have eliminated all the CSV stuff. But sometimes, especially when you’re first learning the details of a library, it’s nice to have the intermediate steps saved out to files so you can satisfy yourself that things are working the way you expect.

A couple of years ago, I wrote a library for calculating section properties. It does a lot of things that Shapely doesn’t, but it also repeats a lot of what Shapely does and isn’t as versatile in handling input data. On my list of future projects is merging the two, so I can get all the important section properties from a Shapely Polygon.