I had a bit of shock this afternoon when I opened my RSS feed reader to see if anything was new.

Not much new, but a lot that’s old. Over 1400 posts from Kieran Healy, holder of the Krzyzewski Chair in Sociological R at the second best basketball university in North Carolina and author of a much-anticipated forthcoming book on how to make good graphs.

What happened? I don’t know for sure, but something in Kieran’s site generation software decided to include every post he’s written in his blog’s RSS feed. It’s an impressive body of work, going back to 2002, but I didn’t have time during my lunch hour to read it all.

My homemade feed reader works like this. For every site I subscribe to, it

2. checks each article against a SQLite database of articles I’ve already read; and

After going through all the subscriptions, the script sorts the unread articles in alphabetical order and arranges them in a static HTML page on my server, adding a table of contents to the top of the page. The script runs via a cron job a few times an hour from 6:00 am until midnight.

So many of Kieran’s posts appeared today because my database of read posts is relatively young and only the last dozen or so of his articles are in it. It was all the earlier ones that were on my feed reader page.

This is my fault, not Kieran’s. I knew perfectly well when I wrote my script that blogging software will sometimes regenerate its feed with all new GUIDs for each article. When this happens, it makes the articles look new to the feed reader. I’d seen this happen even back when I was using professionally written feed reading apps. What made this especially troublesome for my definitely-not-professionally-written feed reading system was that it’s not equipped with a “Mark all as read” button. Which gave me three choices:

1. Do the programming to add a “Mark all as read” button, something I will almost never use.
2. Go through and individually mark all 1400 old posts as read so they get entered into the database and don’t appear again. Fat chance.
3. Figure out another way to add all these posts to the database.
4. Change my feed reading script to just ignore articles that are more than a few days old, regardless of whether they’re in the database.

I chose #4 because it was the quickest to implement and should protect me against this kind of thing happening again. Kieran’s older posts disappeared from my feed reading page, and my blog reading went back to normal. Afterward, though, I realized that I could have implemented #3 in combination with #4, ignoring the older articles for the purposes of assembing the feed reading page but adding them to the database of read articles to give me added protection against seeing them pop up again.

I’ll try to get that working in the next day or two and then post the script in its final form. I doubt that many people really want to set up their own feed reading system, but you never know.

# Subplots, axes, Matplotlib, OmniGraffle, and LaTeXiT

When I learn something new in Matplotlib, I usually write a short post about it to reinforce what I’ve learned and to give me a place to look it up when I need to do it again. In my section properties post from last week, I had a 2×2 set of plots that helped explain which arctangent result I wanted to choose under different circumstances.

Here’s the plot:

And here’s the code that made most of it:

python:
1:  #!/usr/bin/env python
2:
3:  import matplotlib.pyplot as plt
4:  import numpy as np
5:
6:  x = np.linspace(-3, 3, 101)
7:  y1 = (10+6)/2 - (10-6)/2*np.cos(2*x) - 3*np.sin(2*x)
8:  y2 = (10+6)/2 - (6-10)/2*np.cos(2*x) - 3*np.sin(2*x)
9:  y3 = (10+6)/2 - (6-10)/2*np.cos(2*x) + 3*np.sin(2*x)
10:  y4 = (10+6)/2 - (10-6)/2*np.cos(2*x) + 3*np.sin(2*x)
11:
12:  f, axarr = plt.subplots(2, 2, figsize=(8, 8))
13:  axarr[0, 0].plot(x, y2, lw=2)
14:  axarr[0, 0].axhline(y=2, color='k', lw=1)
15:  axarr[0, 0].axvline(x=0, color='k')
16:  axarr[0, 0].set_ylim(0, 12)
17:  axarr[0, 0].set_xticks([])
18:  axarr[0, 0].set_yticks([])
19:  axarr[0, 0].set_frame_on(False)
20:
21:  axarr[0, 1].plot(x, y1, lw=2)
22:  axarr[0, 1].axhline(y=2, color='k', lw=1)
23:  axarr[0, 1].axvline(x=0, color='k')
24:  axarr[0, 1].set_ylim(0, 12)
25:  axarr[0, 1].set_xticks([])
26:  axarr[0, 1].set_yticks([])
27:  axarr[0, 1].set_frame_on(False)
28:
29:  axarr[1, 0].plot(x, y3, lw=2)
30:  axarr[1, 0].axhline(y=2, color='k', lw=1)
31:  axarr[1, 0].axvline(x=0, color='k')
32:  axarr[1, 0].set_ylim(0, 12)
33:  axarr[1, 0].set_xticks([])
34:  axarr[1, 0].set_yticks([])
35:  axarr[1, 0].set_frame_on(False)
36:
37:  axarr[1, 1].plot(x, y4, lw=2)
38:  axarr[1, 1].axhline(y=2, color='k', lw=1)
39:  axarr[1, 1].axvline(x=0, color='k')
40:  axarr[1, 1].set_ylim(0, 12)
41:  axarr[1, 1].set_xticks([])
42:  axarr[1, 1].set_yticks([])
43:  axarr[1, 1].set_frame_on(False)
44:


What was new to me was the use of the pyplot.subplots function to generate both the overall figure and the grid of subplots in one fell swoop. It’s possible that this technique was new to me because the documentation for Matplotlib’s Pyplot API doesn’t contain an entry for subplots.1 I don’t remember where I first learned about it—Stack Overflow would be a good guess—but I’ve since learned that pyplot.subplots is basically a combination of pyplot.figure and Figure.subplots.

Lines 6–10 define the four functions to be plotted. The x values are the same for each and the y values are named according to the quadrant they’re going to appear in. The y values are defined so the moments and product of inertia match the annotations shown in the graph. The actual numbers used in these definitions are less important than their signs and their relative magnitudes, as the plots are intended to be generic.

Line 12 then defines the figure and the array of “axes,” where you have to remember that Matplotlib unfortunately uses that word in a way that doesn’t fit the rest of the world’s usage. In Matplotlib, “axes” is usually treated as a singular noun and refers to the area of an individual plot. After Line 12, the axarr variable is a 2×2 array of Matplotlib axes.

Lines 13–19 then define the subplot in the upper left quadrant (what you learned as Quadrant II in analytic geometry class). Line 19 turns off the usual plot frame, and Lines 17–18 ensure there are no tick marks or labels. Lines 14—15 draw the $x$ and $y$ axes (here I’m using the normal definition of the word). You’ll notice that I’ve drawn the $x$ axis at $y = 2$ instead of $y = 0$. I didn’t like the way the graphs looked with the $x$ axis lower, so I moved it up. Again, this doesn’t change the meaning behind the graph because it’s generic.

The rest of the lines down through 43 are just repetitions for the the other quadrants. Finally, Line 45 saves the figure to a PDF file that looks like this:

Now it’s time to annotate the figure. In theory, I could do this in Matplotlib, but that’s a lot of programming for something that’s more visual than algorithmic. If I were making dozens of these figures, I’d probably invest the time in annotating them in Matplotlib, but for a one-off it’s much faster to do it in OmniGraffle.

I can open the PDF directly in OmniGraffle and start editing. First, I select the white background rectangle that’s usually included in files like this and delete it. It doesn’t add anything, and it’s too easy to select by mistake. Then I select all the axes (again, the usual definition) and add the arrowheads.

The command is very helpful in selecting repeated elements like this.

After placing red circles at the maxima, it was time to label the axes (yes, usual definition; we’re out of Matplotlib now) and add the annotations. I made the annotations in LaTeXiT, a very nice little program for generating equations to be pasted into graphics programs. I’ve been using it for ages.

LaTeXiT cleverly ties into your existing LaTeX installation, so you can take advantage of all the packages you’re used to having available. I usually have LaTeXiT use the Arev package because I like its sans-serif look in figures.

After adding all the annotations, I export the figure from OmniGraffle as a PNG, run it through OptiPNG to save a little bandwidth, and upload it to the server. If this were a figure for a report instead of the blog, I’d export it as a PDF.

1. I’ve complained about Matplotlib’s documentation before, so I’ll spare you the rant this time.

# Canvas and my remote iPad

My older son’s notebook computer, an Asus bought a couple of years ago, has developed a hinge problem that’s reached the point where he doesn’t want to take it to class for fear of it falling apart. After talking over his needs, we decided he could get through the semester with my old MacBook Air. So I set it up for a new user, moved all of my files to an external disk, and delivered it to him yesterday. Coincidentally, on the drive back up through central Illinois, I listened to an episode of Canvas that gave decent explanation of why I could give up my notebook computer.

You could, of course, argue the every episode of Canvas is an explanation of how you can give up your notebook computer. It’s the podcast in which Federico Viticci and Fraser Speirs cover the software and work habits that allow you to use your iOS devices (especially the iPad) to accomplish things you might otherwise think you need a “real computer” to do. But Episode 52 was especially apropos because it covered SSH clients for iOS, which are the reason I feel comfortable in my current state, without a laptop computer for the first time in maybe 25 years or more.

I held off getting an iPad for several years, not because I thought it was a toy or a “consumption only” device, but because my work habits—lots of scripting and command-line use in a multi-window environment—weren’t aligned with the iPad’s strengths. I like to think I wasn’t an anti-iPad zealot during this time. I saw it as the perfect computer for many people, including my wife. I got her an iPad 2 back in 2011; she hasn’t touched a “real computer” since.

So when Split Screen and the iPad Pro were introduced, my ears pricked up. I got the 9.7″ model in late 2016 and have been slowly figuring out how to work with it. Panic’s Prompt and, more recently, the [mosh]1 client Blink Shell are my key apps. My typical setup is to have one of them on the right in Split Screen, connected to my iMac, while I edit in Textastic on the left. This edit/test system on my iPad is very similar to the BBEdit/Terminal window arrangement I use when working on a Mac.

The irony of using a modern, highly graphical device like the iPad to handle a remote, command-line connection to another computer is not lost on me. I often think back to using the Hazeltine terminal that was in a room around the corner from my graduate school office to connect to a Cyber 175 mainframe. And when my iPad is tethered to my iPhone, it’s not unlike using the Hazeltine’s acoustic coupler.

1. Mosh, the mobile shell, is a secure connection like SSH, but is designed to handle more gracefully the interrupted communications common to mobile connections.

# Transforming section properties and principal directions

The Python section module has one last function we haven’t covered: the determination of the principal moments of inertia and the axes associated with them. We’ll start by looking at how the moments and products of inertia change with our choice of axes.

The formula for area,

will give the same answer regardless of where we put the origin of the $x\text{-}y$ coordinate system or how we orient them. You can see this if you think of the area as being the sum of all the little $dx\, dy$ squares in the cross-section.

The formulas for the location of the centroid,

will give different answers for different positions and orientations of the $x$ and $y$ axes, but those answers will all correspond to same physical point of the cross-section.

The moments and product of inertia as we defined them, relative to the centroid,

do not depend on the position of the $x\text{-}y$ origin (because $x - x_c$ and $y - y_c$ measure the horizontal and vertical distances away from the centroid, which is the same for any origin), but do depend on the orientation of the axes. We’ll show how this works by putting the origin at the centroid (which simplifies the math but does not make the results any less general) and comparing the moments and product of inertia for two coordinate systems, one of which is rotated relative to the other.

Note that $\theta$ is the angle from the $x$ axis to the $\xi$ axis and is positive in the counterclockwise direction.1

Because our origin is at the centroid, $x_c = y_c = \xi_c = \eta_c = 0$, and we can write the equations for the moments and products of inertia in a more compact form:

In the $x\text{-}y$ system,

and in the $\xi\text{-}\eta$ system,

We can go back and forth between the two coordinate systems by noting that

Thus,

The $\theta$ terms can come out of the integrals, leaving us with

or

Similarly,

So far, this is just a bunch of algebra that could’ve been done quickly in SymPy. Now it’s time to start thinking.

Looking at the expression for $I_{\xi\eta}$, you might notice that each term includes parts from the double angle formulas. So we can rewrite it this way:

Note that $I_{\xi\eta} = 0$ when

or

Because the tangent function repeats itself every 180°, this expression can be solved with an infinite number of values of $\theta$ that are 90° apart from one another. These orientations all look basically the same, except the $\xi$ and $\eta$ axes swap positions and flip around. For each of them, $I_{\xi\eta} = 0$.

Since we’ve written the expression for $I_{\xi\eta}$ in terms of $2\theta$, lets’s do the same for $I_{\xi\xi}$ and $I_{\eta\eta}$. We start by recognizing the double angle formula for sine in each equation:

Then we use the thoroughly non-obvious identity,2

and use the usual trig identities to get

Therefore,3

and

Now let’s look at how these moments of inertia change with $\theta$. Suppose we wanted to find the $\theta$ that maximized (or minimized) the value of $I_{\xi\xi}$? We’d take the derivative of the expression for $I_{\xi\xi}$ and set it to zero:

This should look familiar. It’s solution is

the same thing we got setting $I_{\xi\eta} = 0$. And if we took the derivative of $I_{\eta\eta}$ with respect to $\theta$ and set it to zero to find the maxima and minima of $I_{\eta\eta}$, we’d get the same thing.

These orientations of the axes are known as the principal directions of the cross-section. They give us both a product of inertia of zero and the largest and smallest values of the moments of inertia. (If $I_{\xi\xi}$ is at a maximum, then $I_{\eta\eta}$ is at a minimum, and vice versa.)

The largest and smallest moments of inertia are commonly called $I_1$ and $I_2$, respectively. They can be calculated by substituting our solution for $\theta$ back into the expressions for $I_{\xi\xi}$ and $I_{\eta\eta}$, but there’s some messy math along the way. It’s easier to recognize that the maximum and minimum moments of inertia are determined entirely by the second and third terms of $I_{\xi\xi}$ and $I_{\eta\eta}$, which are in the form

This expression can be thought of as the horizontal projection of a pair of vectors, one of length $A$ at an angle $\alpha$ to the horizontal and the other of length $B$ at right angles to $A$.

The largest value of this expression will come when the hypotenuse of the triangle, $\sqrt{A^2 + B^2}$, is itself horizontal and pointing to the right. The algebraically smallest value will come when the hypotenuse is horizontal and pointing to the left.

Applying this idea to our expressions for $I_{\xi\xi}$ and $I_{\eta\eta}$, the larger principal moment of inertia will be

and the smaller will be

The axis associated with the larger principal moment of inertia is called the major principal axis and the axis associated with the smaller principal moment of inertia is called the minor principal axis. These are sometimes called the strong and weak axes, respectively. Whatever you call them, they’ll be 90° apart.

Now let’s look at the principal function from the section module and see how these formulas were used.

python:
def principal(Ixx, Iyy, Ixy):
'Principal moments of inertia and orientation.'

avg = (Ixx + Iyy)/2
diff = (Ixx - Iyy)/2      # signed
I1 = avg + sqrt(diff**2 + Ixy**2)
I2 = avg - sqrt(diff**2 + Ixy**2)
theta = atan2(-Ixy, diff)/2
return I1, I2, theta


Looks like I was careless with the $(I_{yy} - I_{xx})$ term and got it backward in the expression for diff, doesn’t it? Also, there seems to be a stray negative sign in the expression for theta. But the principal function does work despite these apparent errors. What we’re running into is the sometimes vexing difference between math and computation.

First, in the formulas for $I_1$ and $I_2$, the diff term gets squared, so flipping its sign doesn’t matter. Second, the numerical calculation of the arctangent isn’t as straightforward as you might think.

There are two arctangent functions in Python’s math library (and in the libraries of many languages):

• atan takes a single argument and returns a result between $-\pi/2$ and $\pi/2$ (-90° and 90°, but in radians instead of degrees).
• atan2 takes two arguments, the $y$ and $x$ components of a vector directed out from the origin at the angle of interest, and returns a result between $-\pi$ and $\pi$ (-180° and 180°), depending on which quadrant the vector points toward.

We can’t use atan in our code because it isn’t robust for some inputs. If we tried

python:
theta = atan(2*Ixy/(Iyy - Ixx))/2


as our formula suggests, we’d get divide-by-zero errors whenever $I_{xx} = I_{yy}$. We can’t have that because there are real cross sections of practical importance for which that’s the case. Any equal-legged angle, for example.

But atan2 can be a problem, too, because we need to distinguish between the major and minor principal axes. In particular, I decided that theta should be the angle between the $x$ axis and the major principal axis. Using atan2 directly from the formula like this

python:
theta = atan2(2*Ixy, Iyy - Ixx)/2


can return an angle 90° away from what we want.

Using the inertia function developed earlier, the moments and product of inertia of the equal-legged angle we just looked at are

Ixx = 9.4405
Iyy = 9.4405
Ixy = -5.1429


Plopping these numbers into the naive formula above, we get

theta = -0.7854


or -45°. This is the angle from the $x$ axis to the weak axis, not the strong axis. The correct answer is 45°, just like the blue line in the figure.

To figure out a way around this problem, let’s plot $I_{\xi\xi}$ for the four cases of interest:

1. $I_{xy} > 0 \quad \text{and} \quad I_{yy} > I_{xx}$
2. $I_{xy} > 0 \quad \text{and} \quad I_{yy} < I_{xx}$
3. $I_{xy} < 0 \quad \text{and} \quad I_{yy} < I_{xx}$
4. $I_{xy} < 0 \quad \text{and} \quad I_{yy} > I_{xx}$

This will let us see what we need for all four quadrants of the atan2 function.

In each of the subplots, successive peaks and valleys of $I_{\xi\xi}$ are 90° apart.

We’re looking for the maximum values of $I_{\xi\xi}$ that are closest to $\theta = 0$, which I’ve marked with the red dots. That means

1. When $I_{xy} > 0 \quad \text{and} \quad I_{yy} > I_{xx}$ (upper right), we want the negative $\theta$ with an absolute value greater than 45°.
2. When $I_{xy} > 0 \quad \text{and} \quad I_{yy} < I_{xx}$ (upper left), we want the negative $\theta$ with an absolute value less than 45°.
3. When $I_{xy} < 0 \quad \text{and} \quad I_{yy} < I_{xx}$ (lower left), we want the positive $\theta$ with an absolute value less than 45°.
4. When $I_{xy} < 0 \quad \text{and} \quad I_{yy} > I_{xx}$ (lower right), we want the positive $\theta$ with an absolute value greater than 45°.

The invocation of atan2 that gives us all of these is

python:
theta = atan2(-2*Ixy, Ixx - Iyy)/2


which we can visualize this way, where the curved arrows represent the angle $2\theta$ for each type of result:

By flipping the signs of both arguments of atan2, we get the sign and magnitude of theta we’re looking for. Note that the expression used in the principal function,

python:
theta = atan2(-Ixy, diff)/2


is equivalent to

python:
theta =  atan2(-2*Ixy, Ixx - Iyy)/2


because of the way we defined diff. And by flipping the signs of both the numerator and denominator, we’re not changing the quotient or the definition of $\theta$. We’re just choosing which solution of

is the most useful.

If you’re mathematically inclined, you may recognize the rotation of axes as a tensor transformation and the determination of principal moments of inertia and principal directions as a eigenvalue/eigenvector problem. But writing principal in those terms would have required me to use more libraries than just math. The formulas in principal are simple, even if their derivation can take us all over the map.

Now that we’ve figured out how principal works, what good is it? It can be shown4 that when the loads on a beam are aligned with one of the principal directions, the beam will bend in that direction only. If the loading is not aligned with a principal direction, the beam will bend both in the direction of the load and in a direction perpendicular to it.

For example, if we were using the equal-legged angle above as a beam and hung a vertical downward load off of it, it would bend both downward and to the left. Not the most intuitively obvious result, but true nonetheless.

Everyone who takes an advanced strength of materials class learns the formulas for the principal moments of inertia and their directions, but there’s usually a bit of hand waving to make the math go faster. And, because the strong and weak directions are typically easy to determine by inspection, the details of picking out the correct arctangent value aren’t discussed. But there’s a richness to even the simplest mechanics, I enjoy exploring it. And since computers can’t figure things out by inspection, you can’t gloss over the details when writing a program.

1. In case you don’t recognize them, the Greek letters $\xi$ and $\eta$ are xi and eta, respectively. They’re often used for coordinate directions when $x$ and $y$ are already taken. You’re probably more familiar with theta, $\theta$, usually the first choice to represent an angle.

2. Don’t believe it? Well, I told you it was non-obvious. But go ahead and multiply out the right hand side and see for yourself.

3. Pay close attention to the negative signs.

4. Don’t worry, I’m not going to show it (not here, anyway). We’re almost done.