Computational Data Science capstone class

In just over a week the CDS capstone class starts and I’m not nearly ready. It’s been a while since I taught a class (Spring 2021 in the pandemic) and I’m a little rusty. But I remembered how useful it is to do syllabus planning on here so I thought I’d give it a try.

[Brief aside to explain why I find this useful: Making little notes to myself can work great, and I can use shorthand for all the things I want to get done. Doing it here forces me to spell everything out, and what I’ve found is that it gets me thinking more deeply about everything. I guess it’s the prospect of you, dear reader, paying attention and not wanting to embarrass myself in front of you, but it really works to get my thoughts more concrete. Plus there’s the added bonus of getting great feedback from you. Finally, when I read these things in future years it helps me to have all those details spelled out, as my cryptic notes tend not to make any sense years later.]

Class context

Our CDS major is pretty new, with only a few graduates so far. This will be the second time this class is taught, but I didn’t teach it last time. This class is typically for seniors in their last semester. It finishes up the major and it meets a couple important general education requirements (oral intensive class and independent learning class). The main goal of the class is to have the students complete a major project using a computational lens on their chosen area of focus. They are required to take at least three courses in that area of focus, with at least one at the intermediate level.

Here are details about the program and the major. The required classes include a python programming class, my Intro to CDS courses (which is a lot more python including web scraping and api managment), a couple stats courses, and a few courses borrowed from the Math (discrete math) and Business Analytics majors.

Organization

Last year the instructor treated each day of the week very differently and I really like the general structure:

  • Mondays: hear from students about a particular data tool
  • Wednesdays: specific support for their projects
  • Fridays: Sample presentations on different data sets

Mondays (tool sharing)

There’s less than 10 students in the class (another reason I’m teaching it) so they’ll all be able to present at least one tool on Mondays. My hope is that we don’t get any repeats and we all develop a stronger tool set. One big question about those days is whether I should try to get them all to be in python, but I’m not sure how much I care about that. I’ll likely do the first day (8 days from now!) so I’ll need to think about what’s a good model.

Could more than one student go in a week? Should they get the whole class period? That’s a lot of pressure but they only have to do it roughly once. Maybe I could model the first few weeks before we really get it going.

I think it should be ~15 minutes of explaining what problems the tool helps with, ~15 minutes about the logistics of using it, ~15 minutes showing it off, and ~15 minutes brainstorming how people might make use of it.

Wednesdays (student project support)

These are the days where I’ll scaffold the work they’ll need to do to really get their projects going. Rather than saying “hey your big project is due at the end of the term, good luck!” I want to make sure that there are near-weekly benchmarks for them to achieve. We’ll not only have them due on those Wednesdays but we’ll make sure that those class times are filled with reflective activities.

My biggest challenge is getting them to be resources for each other, especially since likely no two students will have the same focus area.

Early in the semester we’ll spend time refining what the projects are. The idea of a computational lens on their focus area really needs to be tightly defined so that I’m not constantly making them change direction throughout the term.

Mid-semester will be about helping each other brainstorm deeper ideas, figuring out which things are important and which things are just showing off some tool.

Late semester will be all about refining the final product, namely a video presentation of their project.

Fridays (Data storytelling)

Each week I’ll give them a new data set on Monday. For each they have to come up with a 5-minute live presentation on Friday telling a story with that data. While this takes them away from their project (I likely won’t give them data they’re already working with) it helps them practice the work of computational storytelling and they get to see just how many diverse stories can be told.

One of my colleagues has a ton of useful data sets he’s letting me borrow. A lot of them have to do with the criminal justice system. I also plan to try to be flexible and use timely data that’s in the news. One example that I’m still thinking about is the list of school closings on the most recent snow day we had (last week). What would you do with that?

What else?

I really like that structure but there are some parts that need to be refined. For example, what exactly determines the grade in the class? I’m of course partial to a Standards-Based Grading approach, so let me flesh out what that would look like:

  • A standard (or several) on how to share a tool
  • A standard (or several) on how to learn about a new tool
  • A standard (or several) on how to tell a story with a data set
  • A standard for all the various parts of the scaffolding for the main projects

That last one, especially, is crap. If the project is a bunch of points, then the SBG approach means it doesn’t matter what they do during the scaffolding. I’ll need to think more about that.

Another approach is the contract- or specifications-grading where I lay out the expectations for A-level, B-level, etc work and they work at the pace/quality/etc they want to achieve. That works really well for a big project-based class so I’ll have to give that some thought over the next week as well.

The key, for me, with all of these types of decisions is finding a system that provides accountability and motivation for the students to learn as much as they can. That’s what I’ll be focusing on this week.

Your thoughts

I’d love to hear what you think. Here are some starters for you:

  1. Why has it been so long since you taught? Surely you haven’t moved to the dark side.
  2. Why didn’t you teach it last time? Did they know how much you’d stink at it?
  3. Wait, your general education requirements are embedded in majors courses?
  4. Why does the small class size explain why you’re teaching it?
  5. Of course Mondays should be all python all the time, why would you even consider other things?
  6. Here’s my definition of “computational lens on a focus area” . . .
  7. I’m going to be in this class. I’m most excited about . . .
  8. I’m supposed to be in this class. Now that I’ve read this, how do I drop?
  9. I think it’s dumb to try to get them to help each other out on Wednesdays. Instead you should . . .
  10. Here’s what I’d do with the school closing data . . .
  11. You seriously haven’t decided which grading scheme to use yet? Loser.
  12. Your plan for Mondays adds up to 60 minutes. Surely you actually only teach 50 minute classes!
Posted in computational data science, sbar, sbg, syllabus creation, teaching | 2 Comments

Rolling without slipping on curved surfaces

I’ve been trying to see if I can model balls rolling on curved surfaces and I think I’ve cracked it. Here’s a teaser to get you interested:

What you see is a sphere rolling on a curved surface. The blue line is the path of the contact point and the orange line is the path of the center of the ball. What that ball is doing is called “rolling without slipping” which just means that the contact point doesn’t slide at all. In fact, the part of the ball in contact with the surface is (momentarily) at rest!

Rolling without slipping is something that happens in nature a lot. If something isn’t doing that (like a bowling ball at the beginning of your throw), it often is brought to rolling-without-slipping by friction. Once it gets to that point it tends to stay in that condition. So I thought it would be cool to learn how to model it.

So why did I think this would be an interesting challenge? It seems that introductory physics courses have problems with rolling without slipping all the time. It’s how we learn that rolling balls get down hills slower than the point mass calculations we typically start students with. The reason is that the potential energy of the hill has to be split between the translational kinetic energy that we think of with point particles *and* the rotational kinetic energy it gets while rolling. That’s why it’s slower.

But the reason this is interesting is I wanted to do two things: figure out the angular rotation at all times so that I could make fun animations, and ask what happens if you don’t just do a boring inclined plane — specifically a surface that’s curved.

2D first

First I tackled an effective 2D problem, or a problem I can draw on a sheet of paper. The easiest to do is a non-curved surface that’s at an angle. The rolling without slipping condition in that case is pretty easy:

\omega \left(= \dot{\theta}\right) = \frac{v}{R}

where v is the speed of either the contact point (remember, the part of the ball touching the plane doesn’t move, where that point is does) or the speed of the center of the ball, since they’re the same if the surface is flat.

Now consider the Lagrangian approach for the problem. Let’s assume the angle the line makes with the horizontal is \alpha. Here’s my chicken scratch, barely legible approach:

If you like using x for horizontal and y for vertical, you use the first part. If you want to use a variable along the slope, you use the bottom part (where I used the variable “a”). Either way you get a pretty straightforward result showing that it moves down the slope at a constant acceleration that’s less than g, thanks to the rolling. Note that for the rest of this post I’ll be using the x and y approach because the slope is constantly changing making the “a” approach tricky.

Alright, so what happens when the surface is curved? I ran into this in my Brachistocrone for rolling things post, and the upshot is that the rolling without slipping condition changes to:

\omega=\frac{\rho-R}{\rho R}v_c

where v_c is the speed of the contact point.

Here’s my notes where I proved that to myself:

The key it to realize that while the yellow and orange stripes line up as it rolls, the sphere is frustrated from fully executing its rotation because the curve has risen up to meet it. The part that’s taken back is that d\phi-d\theta that’s labeled.

For that old post the more complicated expression was a major pain in the butt, and I was convinced for a few weeks it would be for this too until I went back to the paper I referenced in that post again a couple days ago and realized that you get a much cleaner version of the rolling without slipping condition:

\omega=\frac{v}{R}

which looks just like my simple one up at the top of this post! Yep, if you just switch to the speed of the center of the ball, instead of the speed of the contact point, you get to use our old tried-and-true rolling without slipping condition. This was the big breakthrough I needed!

Here’s the setup:

  1. Choose the generalized coordinate as x_c (the x-location of the contact point). Assume y_c=f(x_c).
  2. Find expressions for the coordinates and velocities of the center of the ball based on x_c
  3. Using the trick above for the rolling-without-slipping condition, determine the kinetic and potential energies of the ball, at first in terms of the center coordinates but ultimately in terms of x_c
  4. Do the usual Euler-Lagrange trick to find the equation of motion for x_c
  5. Give both the EOM (and the initial conditions) along with the rolling-without-slipping condition to the differential equation solver so that you can get both x_c and \theta as functions of time
  6. Make fun animations

First let’s consider 1 & 2. Here’s an image showing most of the important variables:

How do you get the center-of-ball coordinates from the contact coordinates? The key is the right angle at the contact point. The center is R-units in the direction perpendicular to the surface from the contact point. That direction can be found from the function of the curve:

<x,y>=<x_c-R\frac{f'(x_c)}{\sqrt{1+f'(x_c)^2}}, y_c+R\frac{1}{\sqrt{1+f'(x_c)^2}}>

where f’ indicates the derivative (or slope) of the curve.

Ok, then we literally have the same kinetic and potential energies as what’s in my notes above for the flat curve, only we now know that all the x’s and y’s in there are actually functions of x_c.

Then step four gives us the following differential equation for x_c:

That’s for the special case where f(x)=sin(x). Ugly right? But who cares? We just dump it into an ODE solver! Here’s the result (remember that we get both x_c and \theta back so we can make the motion and the rotation look right):

and here’s rolling down the function f(x)=sin(x)+x:

Cool, right?!

Now for 3D

As soon as I got that working I realized that I could get the 3D working, only I realized that I’d have to shift from also solving for as single angle to solving for all the Euler angles. But, I‘m good at dealing with those so I went for it. Here’s the modified setup:

  1. Choose both x_c and y_c as the generalized coordinates. Assume z_c=f(x_c,y_c)
  2. Find expressions for the coordinates and velocities of the center of the ball based on x_c and y_c
  3. Using the trick above for the rolling-without-slipping condition, determine the kinetic and potential energies of the ball, at first in terms of the center coordinates but ultimately in terms of x_c and y_c.
  4. Do the usual Euler-Lagrange trick to find the equations of motion for x_c and y_c
  5. Give both the EOMs (and the initial conditions) along with the rolling-without-slipping condition to the differential equation solver so that you can get both x_c and the Euler angles as functions of time
  6. Make fun animations

Step 2 needs some adjustment from the 2D case. To find the direction that is perpendicular to a surface at a particular location. The trick is to make a new function:

S(x_c,y_c,z_c)=z_c-f(x_c,y_c)

This function has a set of level surfaces where S is a constant. Our surface is when S=0. But the key is to recognize that the shortest path from one level surface to another is found by finding the gradient of the function. Therefore:

<x,y,z>=<x_c,y_c,z_c>+R\frac{\vec{\nabla}S}{\left|\vec{\nabla}S\right|}

Step 3 looks identical, just with 3 variables instead of 2.

Step 5, especially the Euler angle part, is a little tricky. What is the rolling without slipping condition? We have \omega=v/R but what is the magnitude of \omega? Well, we know what that is in terms of the Euler angles (and their time derivatives)!

\omega_x=-\dot{\theta}\sin\phi+\dot{\psi}\sin\theta\cos\phi

\omega_y=\dot{\theta}\cos\phi+\dot{\psi}\sin\theta\sin\phi

\omega_z=\dot{\phi}+\dot{\psi}\cos\theta

We could just grab their total amplitude, but we actually know something about the direction of \vec{\omega}. We know that it’s perpendicular to both the normal vector to the surface and to the direction of rolling.

\vec{\omega}=\omega\left(\frac{\vec{\nabla}S}{\left|\vec{\nabla}S\right|}\times \frac{<\dot{x},\dot{y},\dot{z}>}{\sqrt{\dot{x}^2+\dot{y}^2+\dot{z}^2}}\right)

Note that “the direction of rolling” is parallel to both the direction that the contact point is moving and the direction the center of the ball is moving, but you’ll see in a second that it’s much easier to use the ball for this (we’re normalizing that vector either way, so it doesn’t matter that they aren’t the same length, it only matters that they’re parallel). Out front, the magnitude of \omega can be found from our simple rolling-without-slipping condition. Putting it all together yields:

\vec{\omega}=<\omega_x,\omega_y,\omega_z>=\frac{1}{R}\frac{\vec{\nabla}S}{\left|\vec{\nabla}S\right|}\times <\dot{x},\dot{y},\dot{z}>

Here’s how all that looks in Mathematica:

One thing that took me a while was figuring out how to make a ball look like it was rolling in Mathematica. Here’s how I did that:

Once you have the prototype ball you just use GeometricTransformation[ball[[1]], {m, r}] on it where m is the Euler rotation matrices dotted together and r is the location of the center of the ball.

Whoo hoo! On to step 6:

Here’s a comparison of 4 balls rolling on the function f(x,y)=sin(x)+sin(y) (it looks somewhat like an egg carton). They’re all started at the same point. Two are small and two are larger. In each size one is a solid ball and one is a shell. They all have the same mass:

And here’s the biggest ball I could get to roll on that surface. Why would there be a limit? Because you can’t use a ball that has a larger radius than the concave up radius of curvature at any point on its journey (Mathematica actually craps out right when this happens):

Here’s a still frame zoomed in a little showing that biggest ball right at the smallest radius of curvature:

One cool thing is that I can have Mathematica calculate the radius of curvature at all points along the trajectory. Here’s that plot where the red horizontal line is the radius of the ball:

How do you calculate the radius? You get it from this equation from way at the top of this post:

\omega=\frac{\rho-R}{\rho R}v_c

Note that before I really got my head wrapped around how to do this I thought I would have to calculate that radius of curvature (which, by the way depends on both where you are on the surface and what direction you’re rolling) directly from the definition of the curve itself. That was an interesting rabbit hole that this approach didn’t need (sorry for bothering you, my friend (and occasional guest blogger Art) for how to do this).

Your thoughts?

I’m really excited about this new type of modeling I can do. What are your thoughts? Here are some starters for you:

  1. This is cool! What I especially like is . . .
  2. This is dumb. I could do this on a napkin.
  3. Why do you sometimes use hypens and sometimes not in the phrase ro-lli-ng—-wit-hou—-t-sli–ppin———-g?
  4. Can you model a ball rolling on a rotating turntable? (actually this is what started all of this, I just haven’t gotten back to it yet)
  5. What is up with that [[1]] you have to add to the GeometricTransformation command?
  6. Will this make you a better golfer?
  7. I like your handwriting better than \LaTeX. Do that from now on.
  8. I hate your handwriting.
  9. That level curve method is exactly what you did in 2D too, idiot.
  10. This is Art: Are you telling me you didn’t use any of those calculations I gave you?!
  11. What happens if the surface changes with time?
Posted in fun, general physics, mathematica, physics | Leave a comment

Situations that share equations of motion

Recently my friend Rhett Allain has been making some awesome videos showing how to solve complex problems with a Lagrangian approach. I love it when he posts a new video because it usually motivates me to try to model something similar.

Here’s the video he posted that inspired this post:

Rhett’s vid about two masses tied together through a hole in a table

I love this particular problem so I decided to dig in and see what I could learn from building my own model of it. One thing I was curious about was whether using a Lagrange multiplier technique could save me time. Basically I wondered if I could just model the two balls quite simply using Cartesian coordinates (xyz coordinates) and then directly imposing a constraint on the length of the string between them.

Unfortunately my first attempt had a huge mistake in it leading me to make an erroneous claim in this response video:

My initial video response to Rhett that has a big mistake in it (can you spot it?)

I got so excited that I had found an analogous system that shared the equations of motion, only upon reflection I realized I’d made a big mistake.

The length of the string should be given by the length above the hole:

\sqrt{x^2+y^2}

along with the length of the string below the hole:

\left|z\right|

My mistake was assuming that was equivalent to:

\sqrt{x^2+y^2+z^2}

Do you see the mistake? I hope it takes you more than a second or two to see it because I spent half a day thinking I’d discovered something really cool! Note that I still wrote this post so, trust me, there is still something cool, just not what I thought it was!

So what was I thinking? Well, let’s look at the code:

First let’s assume the masses are the same. Then the kinetic energy (the line starting with “T=”) would just be the kinetic energy of a free particle. In other words this would be just the model of a ball flying through a simple gravitational field. But that then gets corrected with the constraint line (the one that starts “cons = “). That says that the sum of squares of the variables needs to be a constant. That would imply that the particle has to stay on (or inside of) the surface of a sphere.

That’s what I got excited about in my video, because it seemed like the problem of the two balls connected through a hole with string was identical in form to a ball that has to stay on a sphere in a gravitational field. In hindsight I realize that a better way to describe that latter analogous problem is to call it what it is: a 3D pendulum!

But, ugh, I got the constraint wrong. It should have been that the length of string (which is held constant) is given by:

\sqrt{x^2+y^2}+\left|z\right|

Of course in my case I put the hole at z=0 with “up” representing positive z’s. In other words the z in this problem will always be negative. So that gives:

\sqrt{x^2+y^2}-z=\text{constant}

Making that correction to the “cons = ” line gives me the correct model of the two particles tied together by a string running through that hole.

Now, here’s the cool part. If you set the masses equal, you can interpret the kinetic and potential energies (along with the constraint) to be describing a system where a ball has to stay on a surface described by that constraint.

Here’s what that looks like side by side:

See. I told you. Cool!

Your thoughts?

So what do you think? Here are some starters for you:

  1. This IS cool. What I especially liked was . . .
  2. Um, this isn’t cool, thanks for wasting my time. What I especially disapprove of is . . .
  3. Let me get this straight: in an effort to look smart you posted a youtube vid that has a huge error in it and still thought you might as well post that mistake-ridden vid in this post. Loser.
  4. This is Rhett. You’re awesome.
  5. This is Rhett. Stop stalking me.
  6. Why do you use the Lagrange multiplier approach? Are you afraid to type sines and cosines or something?
  7. Lots of systems share a set of equations of motion. Why do you think it’s cool?
  8. It took me exactly 2.3 seconds to find your dumb mistake.
Posted in mathematica, physics | 3 Comments

Rigid bodies, formulation and examples

My friend Rhett Allain gave me a good challenge recently with this tweet:

I had been working on a problem that he posted about regarding a bead sliding freely on a hoop that is spinning about an axis in its plane that goes through the center. That’s a pretty typical Lagrangian Dynamics-type problem and I wondered what would happen if the hoop wasn’t driven to constantly go with a particular angular frequency but rather was set spinning on that axis with that initial angular frequency but after that it would respond dynamically (speed up or slow down its rate of spinning) according to the physics involved. Here’s what that looks like:

But then when Rhett asked about letting the hoop be even more free to rotate, in other words not just around that one axis, I realized I’d have to calculate things a different way.

That led me to remind myself how to do Lagrangian Dynamics with rigid bodies, mostly dealing with things like the Inertia Tensor and body axes. This post is to make it so I don’t have to start from scratch with that again. Note that this post lays out the broad strokes of what I’m doing here.

Lagrangian dynamics with rigid bodies

For those of you hoping I’ll derive the Euler rotation equations, I should let you know that I don’t explicitly do that here. I get close, but I really move in a direction where I can solve the equations of motion and make some fun vids of the motion.

Here’s the steps involved:

  1. Determine the generalized coordinates
  2. Determine the kinetic and potential energies
  3. Use the Principal Axes (actually you don’t have to do that)
  4. Do all the usual Lagrangian stuff

Generalized coordinates (Euler rotations)

If you need to describe the orientation of a rigid body, how would you do it? I’ve done this with my students by asking them to all grab their chairs and figure out a systematic way to orient them the way mine is (which I do when they have their eyes closed). At first students think you just need the polar angles \theta and \phi because they figure that all the r’s are fixed and so you don’t need that third polar coordinate. But then I show them two situations with my chair that share the polar angles but look different. Ultimately my students usually land on a version of the Euler angles, which are really just the polar angles with a final twist of the body around the original vertical axis. Most texts call that last angle \psi.

Well, if the Euler angles are the information that you need to describe the orientation of a rigid body, that means they’re just the perfect generalized coordinates for a Lagrangian Dynamics approach. All we have to be able to do is describe the potential and kinetic energy of the object as functions of the Euler angles and we’re good to go!

Kinetic and Potential energy as functions of the Euler angles

Early on physics students learn that if they’re dealing with rotating bodies they can replace \frac{1}{2} m v^2 with \frac{1}{2} I \omega^2. Essentially what happens is that you realize that all parts of the system share an angular speed (\omega) even if different parts are traveling with different linear speeds (v). Since they all share the same thing, you can pull that out of any sums or integrals you’re doing to add up all the kinetic energy and whatever’s left you just lump into something we call the moment of inertia (I).

The problem is that for oddly shaped things that simple approach doesn’t quite work. Sure all the pieces and parts still share an angular velocity (now we’re talking about the vector \vec{\omega} since the axis of rotation is something we need to know), but it’s not as simple as \frac{1}{2} I \omega^2. Probably the easiest way to see that is that the angular momentum is not necessarily parallel to the angular velocity (STOP: read that again, as it’s one of the single most important aspects of rigid body motion and it sounds really weird to most people). Instead there’s a relationship between the angular momentum vector and the angular velocity vector that’s more complicated. If you know the angular velocity you can find the angular momentum with a 3D rotation, and that’s done computationally with a 3×3 tensor which we call the inertia tensor. Then you have \vec{L}=\overleftrightarrow I\cdot\vec{\omega}. Essentially a 3×3 tensor dotted into a vector gives you a new vector that (often) points in a new direction (hence them not being parallel anymore).

Really the inertia tensor has just what we need. Just like we can do \frac{1}{2}\vec{p}\cdot\vec{v} for kinetic energy normally (try it!), we can do \frac{1}{2}\vec{L}\cdot\vec{\omega} or:

\text{Kinetic Energy}=\frac{1}{2}\left(\overleftrightarrow I \cdot \vec{\omega}\right) \cdot\vec{\omega}

Aha! So now if we know how both the inertia tensor and the angular velocity vector can be described by the Euler angles (and their time derivatives) we’re in business!

The inertia tensor is certainly calculate-able if you know the orientation of the object (which you do if you know the Euler angles). Essentially you can figure out where all the particles of the body are and then do all the integrals necessary for the inertia tensor (remember, the tensor is just made up of what’s left after you factor out the common terms that all the particles share – namely the angular velocity – so it’s just a collection of sums/integrals involving the locations of the particles). What really sucks, though, is that at every point in time there is a new orientation and hence a whole new set of integrals to do. Ugh. Don’t worry though, that gets a lot easier in the next section!

The angular velocity is certainly a function of the time derivatives of the Euler angles since as they change they rotate the body, and rotation is all that angular velocity is used to describe. Here’s the cool part: the angular velocity (\vec{\omega}) is equal to the sum of three vectors, one for each Euler angle. Each Euler angle contributes it’s time derivative (which is exactly the right units for angular velocity) in the direction of the instantaneous axis of rotation that Euler angle rotates around. I know, that was complicated, but it’s relatively straightforward to work out:

First consider the easiest one: \phi. When you think about it that polar vector, all it does is rotate everything around the vertical-, or z-, axis. Therefore \phi‘s contribution to \vec{\omega} is \dot{\phi} \langle 0, 0, 1\rangle.

Next consider \theta: When it changes, it tilts the object further from the vertical. But it does so in the plane determined by \phi. The instantaneous axis of rotation is the location of the new y-axis after the \phi rotation. So that means that the \theta contribution to \vec{\omega} is \dot{\theta}\langle -\sin\phi,\cos\phi,0\rangle

Finally we have \psi. It rotates the body around its original vertical axis. That means its instantaneous axis of rotation is wherever its vertical axis has gotten to after the polar angles have done their thing. So \psi‘s contribution is \dot{\psi}\langle \sin\theta\cos\phi, \sin\theta\sin\phi,\cos\theta\rangle.

All together then we have these for the components of the angular velocity vector:

\omega_x=-\dot{\theta}\sin\phi+\dot{\psi}\sin\theta\cos\phi

\omega_y=\dot{\theta}\cos\phi+\dot{\psi}\sin\theta\sin\phi

\omega_z=\dot{\phi}+\dot{\psi}\cos\theta

Only calculate the inertia tensor once!

Ok, we have the angular velocity all set to go, but we saw how the inertia tensor, while definitely a function of the Euler angles, might have to be calculated over and over again. What if we only had to do it once!

The trick is to decide to think about how things feel in the frame of the rigid body. In what we call the “body-frame” it sits there while the world moves around it. Specifically the angular velocity vector is all over the place. But, and this is the important part, the inertia tensor is fixed! You just have to calculate it once in that orientation. So the problem becomes determining what the angular velocity vector (which we calculated above) looks like in that body-frame.

But luckily the Euler angles were built to transform any vector into a new rotated state. So we just need to apply the Euler rules to the angular velocity vector above. Except that we’re trying to describe how things look to the body frame, so what we really need is the reverse of the Euler rules, namely:

  • negative rotation of \phi around the z-axis, followed by. . .
  • a negative rotation of \theta around the y-zaxis, followed by . . .
  • a negative rotation of \psi around the z-axis

I know, that seems weird, but, trust me, it does the trick. All of those rotations are pretty straightforward for something like Mathematica to do so we’re in business: Calculate the inertia tensor once, figure out the lab-based angular velocity vector, and then do the steps above to find the body-frame version of the angular velocity vector. Then you calculate the kinetic energy.

Here’s what that looks like in Mathematica:

After the “KE=…” line is the usual construction of how I do Lagrangian Mechanics. The “Method->…” stuff was suggested by Mathematica after the first solve seemed to fail. That Method fix did the trick.

You’ll note that my inertia tensor (“Ithing”) is really simple. There are no off-diagonal terms in the matrix because I wanted to investigate the same sort of thing I did in this post about flipping handles in space. The section below called “Principal Axes frame” is how I teach my students to find the magical orientiation of an object where its inertia tensor has that simple setup. Note that it’s not necessary for the work I’ve laid out here.

Note that while this makes it pretty do-able, you’re still going to want a Computer Algebra System to help. Here’s just the \theta Euler-Lagrange equation (all of this equals zero) for a generic inertia tensor:

Putting it to use to answer Rhett

Ok, so now I had the tools to answer Rhett’s question. I’ve got an interesting system made up of a hoop and a bead. What I decided to do was to fix \psi for the hoop at zero because otherwise I would have to keep track of the bead’s sliding separate from the rotation of the hoop around it’s (original) vertical axis. The orientation of the hoop when all the Euler angles are zero is horizontal. So I first set \theta to \pi/2 to stand it up and then I set \dot{\phi} to some initial spin rate.

For the bead I just determine it’s contribution to the inertia tensor of the system. Luckily that’s just a simple function of \psi (remember I fix that at zero for the hoop). So I basically do the same thing as before. Here’s the Mathematica:

Really I have to carefully do the hoop and the bead separately when doing the \omega calculations, but otherwise it’s similar to what I did above.

Here’s the results!

You can see how the hoop starts to tilt over. It has to to maintain a constant angular momentum (since there’s no torque in the problem if you don’t fix the axis of the hoop).

Principal Axes frame

I noted above that you don’t have to use the principal axes frame but I wanted to put this in for my future self:

When I teach this I point out how difficult things are when the angular moment is not parallel to the angular velocity. I ask them if they’d like it if we could make that disappear. Usually they say yes.

I tell them that there are some magic scenarios when the two are parallel. Essentially there are three axes that if you rotate around them you get an angular momentum in that same direction. There are a few ways to figure them out, but most are either difficult to understand mathematically (eigenvectors of the inertia tensor) or would take forever (just spin the object along all possible axes and look to see when the angular momentum is along that same axis). But there’s one way I’ve found that my students can do. I have to explain how a diagonalized inertia tensor does the trick. I do that by saying that we’ll just keep orienting the object until we find that all three of the magic axes are along the normal lab axes (x-, y-, and z-axes). In that orientation, you get what we need as long as there are no non-zero off-diagonal components of the inertia tensor.

So I provide my students with a Mathematica document that lets them freely adjust the Euler angles for an object and it constantly recalculates the inertia tensor for that orientation (see the last section above to understand why the tensor changes with orientation). What’s cool is that they can see how certain adjustments to the Euler angles tends to reduce at least one of the 3 off-diagonal components and they usually quickly find a systematic way to get them all to be vanishingly small. When they do that, they’re staring at the orientation that puts the magic axes along the lab axes!

Your thoughts?

I’d love to hear what you think. Here are some starters for you:

  • I love this! What I especially am surprised by is . . .
  • This is dumb. I especially hated . . .
  • So whatever Rhett says you just do? That sounds interesting . . .
  • Once again no python, loser
  • I thought you only had a chromebook at home. How are you doing all this Mathematica?
  • Whoa, it looks like you do the \psi rotation first in your Mathematica code. Didn’t you say it comes last?
  • Wait, I don’t have to diagonalize the inertia tensor before doing these kinds of simulations? That’s so cool!
  • Do you provide tools to help your students unbolt their chairs from the desks?
  • I want to hear more about how students can systematically find the principal axes by just adjusting the Euler angles!
  • Why don’t you have to do the normal accelerated frame adjustments for your kinetic energy? (it’s quite late as I type this and I realize I don’t have a good answer beyond “it just works”)
Posted in mathematica, physics, syllabus creation, teaching, twitter | 1 Comment

Creating bike routes with python

This weekend my goal was to ride 50 miles to and from my house. In my last post I showed four ways to find where I could get to for a certain distance, but I really hate “there and back” rides, so I wanted to find loop-based routes that would have my target distance and not have any doubling back. I used basically the same tools and I’m decently happy with the results.

tl;dr here’s the route

How’d I do it?

Here’s the basic gist of what I did:

  1. Get the lat/long of my house to have a starting point
  2. Do some geo-based math to find points on a regular polygon that includes my house. I set the polygon perimeter length to my target length
  3. Ask openrouteservice to find directions to all of the points in order (beginning and ending at my house)
  4. Adjust the total polygon perimeter until the actual path is my target distance
  5. Construct a google maps URL so that I can navigate while riding

Step one: lat long of my house

This one is easy. Just open google maps and right click anywhere to get the latitude and longitude

After doing that here’s what’s in your clipboard: 44.9240104020827, -93.09348080561426

So just paste that somewhere in your python code as a list and you’re good to go.

Step 2: Geo-based math for other points

I’m pretty sure I could have used some package to do this, but I figured it wouldn’t be that hard to just code up a sloppy version, especially since I knew I’d have to make adjustments to the polygon positions so that the actual travel distance would be what I wanted (see step 4).

If you have the lat/long of a point and want to find a point some known distance away in some particular direction (say \theta degrees counterclockwise from east), you need to realize that while moving north/south follows a great circle (with a radius equal to the earth’s), moving east/west is not moving along a great circle (which is why you have to turn slightly north when driving straight east). Instead you’re moving on a circle whose radius is the earth’s radius multiplied by the cosine (ugh) of the latitude. If you were thinking it should have been sine (like my first thought) it’s likely because you’re used to using polar angles in mathematics (ie, the zero of latitude is the equator, not the north pole). Also, as you’ll note looking at the lat/long I pasted in above, since I’m roughly at 45 degrees latitude, it doesn’t really matter.

So here’s how I coded it:

earthRad=6378.1
def newPoint(start, angle, distance):
  [long,lat]=start
  horizScale=earthRad*math.cos(lat*math.pi/180)
  newLong= long+distance*math.cos(angle)/horizScale*180/math.pi
  newLat=lat+distance*math.sin(angle)/earthRad*180/math.pi
  return [newLong,newLat]

Ok, so now I just need to make a bunch of points on a regular polygon:

def polygon(sides,start,totalDistance,initialAngle=0):
  theList=[start]
  curLoc=start
  for x in range(sides):
    curLoc=newPoint(curLoc,x/sides*2*math.pi+initialAngle,totalDistance/sides)
    theList.append(curLoc)
  return theList

You might wonder about that “initialAngle” part. I found that being able to tilt the initial polygon helped because I have a big river just east of me that the mapping software doesn’t like (only if you try to do bicycle navigating right into the river).

Step 3: Use openrouteservice to map it out

Here’s the code I used:

import openrouteservice as ors
from openrouteservice import convert
import folium
import math
token="GET YOURS FROM OPENROUTESERVICE"
client=ors.Client(key=token)
home=[PASTE THIS IN FROM GOOGLE MAPS]
directions=client.directions(profile='cycling-regular', coordinates=polygon(7,home[::-1],60,3*math.pi/4))
directions['routes'][0]['summary']['distance']

You can get your own token by signing up (for free!) at https://openrouteservice.org/ It has some limitations regarding how often you can use it, but they’re very generous, at least for this kind of thing.

The last line is what you have to do to get the total distance for the journey. You can see that the result of the “client.directions” command is a highly structured object (it took me a while to figure out how to extract that total distance). You get the map I generated above with this:

map=folium.Map(width=500, height=500,location=home)
folium.GeoJson(convert.decode_polyline(directions['routes'][0]['geometry'])).add_to(map)
map

Now you see why I had to import folium up above.

Step 4: repeat and adjust until you get the distance you want

You’ll note in the code above that I submitted points on a polygon that would have a total perimeter of 60 kilometers, but it returned a path that’s just over 80 km, or 50 miles (my original goal). The reason is mostly due to the taxi-cab geometry I talk about in my 2nd approximation in this post. In other words, there’s no way there are dead straight roads among all the polygon vertices, so you’re going to travel further.

Step 5: create a google maps URL

While I like the programmatic capabilities of openrouteservice, I *love* google maps as a navigation aid while riding. So I wanted to figure out how to get my newly found route into google maps. Solution: use their url api! On that page you can learn how to create URLs with starting points, ending points, way points (points along the journey), and type of directions (I wanted bicycle directions). Here’s how I built that:

pgon=polygon(7,home[::-1],60,1*math.pi/4)
baseUrl="https://www.google.com/maps/dir/?api=1&"
start="origin="+str(home[0])+","+str(home[1])
travelmode="travelmode=bicycling"
waypoints="waypoints="+"|".join([str(x[1])+","+str(x[0]) for x in pgon])
destination="destination="+str(home[0])+","+str(home[1])
baseUrl+start+"&"+travelmode+"&"+waypoints+"&"+destination

That gives you a URL that produces this:

And I’m finally ready to ride!

How it went (the actual ride)

It was a beautiful late October day and I really did have a blast. The mapping worked really well but there are some issues I had to deal with.

If you look at the image above, you’ll probably note the little spur on the eastern edge of the route. Essentially google’s mapping algorithm (really both algorithms if you go back and look at the first image on this post) didn’t know of a way to get to and from that particular waypoint without a little bit of doubling back. I didn’t want to do that, so I just figured I could delete that waypoint (google makes that pretty easy to do). However, if you delete it before you begin, google will try to find a faster way from the prior waypoint to the next one, so I just waited until I got close before deleting it. That worked well.

The other waypoints weren’t so far off a normal path, but I did do a lot of weird extra blocks to get to the exact point I asked for. I probably could have saved a few miles not doing that, but that wasn’t really the point.

Overall I think it worked out pretty well. I’m excited to do it again. I’ll probably just tilt the 7-sided polygon over just a little further to get a completely different ride.

Your thoughts?

Got any feedback for me? Here are some starters for you:

  • I love this! Can you do it for . . .?
  • I hate this! It’s obvious you could do this in a much easier fashion by . . .
  • What do you mean you have to turn a little north to drive straight east?
  • I think we should redo latitude so that the zero is at the north pole
  • Why python, I thought you loved Mathematica?
  • Why are the google and openrouteservice directions slightly different?
  • Why only 50 miles?
  • How long did it take you? (8 hours but that was with a really long lunch and a few other stops to read and watch rugby)
  • How did you find rugby to watch?
  • What’s wrong with using cosine (you said “ugh” after doing so)?
Posted in fun, math, programming, technology | Leave a comment

Gabriel’s Horn (guest post)

Last week I was intrigued by this post:

I asked, on Facebook, whether filling it with paint would essentially be painting it on the inside and I had a suspicion that my good friend Art Guetter (Professor of Mathematics at my institution) would help me learn something. I was right, and he’s been kind enough to type up his thoughts for this post. Here he shows that it’s pretty tough to use a finite amount of paint to paint anything infinitely long.

Guest post from Art Guetter

Gabriel’s horn is a solid created by rotating the graph of f(x) = 1/x, defined on the interval (1,\infty), around the x-axis. Two “paradoxes” arise. In the first, the horn has a finite volume, despite being created by rotating a region with infinite area around an axis. The second is that the horn has finite volume but infinite surface area, leading to the apparent paradox that the horn could be filled with a finite volume of paint, yet the paint would not be sufficient to coat (that is, paint) the surface. The resolution of this “painter’s paradox” is that the thickness of the paint would need to decrease to 0 in the limit as x tends to infinity. The assumption being made here is that painting requires a uniform thickness of paint. Note that I can paint the entire plane if I am allowed to decrease the thickness of paint as I move far from the origin.

So could I paint an infinite solid of revolution (to a uniform depth h) if the surface area were finite? As a first example, replace f(x) = 1/x from Gabriel’s horn with a piecewise constant function f(x) = r_n when n < x < n+1 for n = 1,2,3,\ldots, and the constants r_n to be determined later. The surface will consist of an infinite collection of right circular cylinders, and each cylinder will have surface area 2 \pi r_n. If the r_n are chosen so that the sum \Sigma r_n < \infty, can I paint the surface with a finite amount of paint? The answer appears to be “yes”, but this involves the assumption that I roll each cylinder open, so that the amount of paint used is simply the surface area multiplied by the thickness of the paint, say h. (Each cylinder can be rolled open without issue because they have thickness 0.)

What about painting if the cylinders aren’t rolled out? I will assume that painting to a thickness h means that the depth of paint at any point is h measured along the normal, in the outward direction. The amount of paint needed to paint one of the cylinders is then given by \pi [(r_n + h)^2 - r_n^2] = \pi [2 r_n h + h^2] = 2 \pi r_n [h + h^2/(2r_n)] = A_n [h + H h^2], where A_n is the surface area of the cylinder and H = 1/(2r_n) is the mean curvature of the cylinder. Summing this over n will lead to an infinite volume of paint, no matter how fast the r_n tend to 0.

A more general theorem has that the volume of a surface that has been thickened by an amount h in the direction of the normal to the surface (assuming that h is small enough that there is no self-intersection) is given by

V = h \cdot SA + h^2 \cdot \int H \; dA +\frac{1}{3} h^3 \cdot \int K \; dA

where SA is the total surface area, dA is the surface element, H is the mean curvature of the surface, and K is the Gaussian curvature. (These are constant for the cylinder, with values H = 1/(2 r_n) and K = 0.) The amount of paint needed to paint a surface to uniform thickness depends on the curvature of the surface.

For a surface created by revolving the curve y = f(x) around the x-axis, the values of H and K depend only on x and are given by (in general and then for f(x) = 1/x)

H = \frac{1 + f'(x) - f(x)f''(x)}{2 f(x) (1 + f'(x)^2)^{3/2}} = \frac{x^7 - x^3}{2(1 + x^4)^{3/2}}

K = -\frac{f''(x)}{f(x) (1 + f'(x)^2)^2} = -\frac{2 x^6}{(1 + x^4)^2}

Posted in guest post, math | 1 Comment

What’s my 30 mile cycle limit?

UPDATED WITH 4th APPROXIMATION!

Last weekend I went hammock camping by towing all my gear behind my bike. I loved it and now I’m interested in finding other adventures that won’t tax me too much. I really think that, for now, 30 miles in one day towing the trailer is a good limit for me. It leaves me enough energy to make camp and I’m able to relax and enjoy myself.

1st approximation

My first thought was to just look at a map with a 30-mile radius circle centered on my house. I figured if I could find any campgrounds in that circle I’d be good to go.

30 mile radius circle around my house

The problem with this approach is that there aren’t any roads or paths that go straight from my house to the edge of this circle. Any way I’d ride out to that edge would be more than 30 miles of riding.

2nd approximation

I’ve noticed that I spend most of my time riding in the Cardinal directions (North, South, East, and West). What If I could figure out how far I could ride only going in those directions?

I realized that I could do a polar plot around my house if I could find a relationship between r (the radius) and theta (the angle with respect to East). It took a little head scratching, but here’s what I came up with:

r \cos(\theta)+r\sin(\theta)=D

or

r(\theta)=\frac{D}{\cos(\theta)+\sin(\theta)}

Actually that only works in the first quadrant. For all the quadrants you want the absolute value of both the cosine and sine term:

r(\theta)=\frac{D}{\left|\cos(\theta)\right|+\left|\sin(\theta)\right|}

A polar plot of that equation with D=1

Weird and surprising, right? I hope it was for at least a few of you. It surprised me! But of course, after some thought it makes some sense. As you give up height you get exactly that much more width. So a straight line with a slope of 45 degrees it is!

Now with the Cardinal direction limit.

So why doesn’t the blue square touch the red circle, you might be asking? Well, I’m not really sure. I got the locations of those 4 corners from Google Maps using the “Measure Distance” feature, whereas the red circle comes from the Circle command in python’s folium library. It’s weird that they’re so far off from each other, isn’t it?

The obvious problem with this solution is that not all roads/paths run along the Cardinal directions. I think that’s likely even more true for bike paths.

3rd approximation

I realized that Mathematica has the command TravelDirections that lets you put in two locations and a TravelMethod to get decent directions. I used “Biking” for the TravelMethod for all of this work.

I had Mathematica get directions for 36 evenly spaced locations on that red circle above. Then I took a look at the actual travel distance and they were all well over 30 miles. So then, for all 36 of those, I crept in towards my house by 1 kilometer until I got a path that was 30\pm1 miles away. Then I just plotted them all to get a decent sense of my options:

The gray is a 30 mile radius disk and red are all the bike journeys I calculated

So I think this is my best approximation yet. It gives me a real good sense of how far I can get comfortably in one day. I haven’t really investigated the directions, but certainly I can see the exact path I took this past weekend in there (it’s the one that heads just a little west of south)

UPDATED: 4th approximation

Ok my mind was just blown with something called isochrones in the OpenRouteService available to python. Here’s my colab. And here’s the amazing result:

Isochrone approach limited to 30 miles of biking distance

Your thoughts

So what do you think? Here are some starters for you:

  1. This is cool! Have you thought of doing this . . .?
  2. This is dumb. Taxi drivers figured this out years ago
  3. Why do you only sometimes capitalize the carDinAl directions?
  4. Could you please share the colab doc you used for the folium maps?
  5. Could you please share your Mathematica code (saved at work and I’m typing this at home for the moment)?
  6. I live ____ and can ride ___ miles in a day. Could you do this for me?
  7. I know why the blue and red don’t touch, it’s because . . .
  8. That blue square didn’t surprise me at all. You’re dumb.
  9. That blue square totally surprised me! I learned something today!
  10. You’re going to tease us about your bike and not bother showing a pic?
my bike (Priority 600) and my homemade trailer
Posted in fun, math, mathematica, technology | 6 Comments

Classroom photo sharing app

For a long time I’ve wanted an app that could

  • Allow my students to take a picture of their work and share it with the class
  • Certainly my computer should be able to display it, but with thumbnails for all the images
  • Bonus if all the images are on everyone’s device

This post is about my attempt to make just such an app in Google Apps Script. It includes detailed instructions for how you can make your own copy. If you want to see it in action, see this vid.

Before jumping in, I want to say thanks to some great twitter friends who had thoughts about other ways to do this. I haven’t tried them all, but I think they all could work really well:

How to get it working for you

Below I’ll explain a little more about how it works, but one important feature is the use of webSockets. To do that I use a free account on pusher.com. It limits all my apps to 100 total simultaneous connections. For the purposes of this app, that means 100 students all sharing images with each other. That’s fine for me, since I’ve never taught a class that big, but I can’t just let everyone use my copy of the app since then I’ll hit that limit pretty quick. The good news is that if you follow the steps in this section you can have your own version up and running that’ll have its own 100-user limit.

So: Step 1, sign up for a pusher account. I described how in this post.

On to step 2. Make your own copy of this google sheet. Make a google folder to hold the images and make sure its contents are viewable to the world. Then follow the instructions in this vid (:

I made a copy in my consumer gmail account

And you’re done! Have fun!

To start the next class all you have to do is clear out the spreadsheet. You can delete the files in the image folder too if you want, but you wouldn’t have to.

Ways to use it

I think my ideal situation is when I ask students to work on something, either individually or in groups, and then letting them all review everyone’s work. Without this app I’ve done that in the past by gathering their personal whiteboards and showing the class what I’ve found or asking them to walk around and look at other groups’ whiteboards. This way they see everyone’s and they can zoom in all they want without trying to see past people. It’s somewhat anonymous since the system doesn’t capture who uploads the image (actually it does capture the email and the date but I haven’t bothered to use that info in the UI).

Here’s an example: “Ok everyone, write down something you know about the dot product. When you’re done, submit your photo.” It seems to me you’d get some cosine (ugh) thoughts, some vectors saying “hey you, how much of you is parallel to me”, and some component-by-component scribbles. Imagine the students scrolling through all of that and then asking what patterns they’re seeing!

I also think this could be interesting in other situations. Imagine if your students were out and about doing some observations or something. They could stay on top of what everyone is seeing with this app. In other words, they don’t have to be in the same room to use the app.

How it works

What’s stopped me in the past from making this app is the annoying habit of phone manufacturers where they make their cameras have incredible pixel densities. That means that the typical file sizes are multiple megabytes. That’s not really a big deal for your phone, but it means your uploads take a while and the server just gets hammered.

This time around I googled a little and found that you can resize images before uploading them. The local machine (your student’s phone) creates a much lower density copy and uploads that instead of the original. This also makes the final product much more responsive.

Next comes the hassle of uploading and saving files in google apps script. Here’s a chapter of my GAS book explaining how I do that, though I did a little different trickery due to the copy of the image that gets created by the local machine.

Next comes making a thumbnail-based image viewer. I just used what I found here. Of course I had to make it a little more dynamic, essentially recreating the thumbnails when a new image comes in. I don’t bother to make tiny images for the thumbnails, they’re just the full image shown small. But since they’re all pretty small to begin with, that works fine.

So all together then, here’s the life cycle:

  • A student loads the page
  • They’re shown any images that have already been uploaded.
  • If they upload a new one:
    • Their machine makes the low density version
    • and then sends it to google.
    • Google saves the image
    • retrieves the download url
    • saves the url in the spreadsheet (which is why the page loaded all the old ones to begin with)
    • and sends that url to all connected devices through pusher
  • When a pusher message is received, the student’s machine requests the new image from the google server
  • Then that image is added to the array of images at the top of the page.

Some issues:

  • For some crazy reason, this only works on iPhones if they use an incognito window
  • If for sure everyone loads the page at the beginning of class, you could skip the spreadsheet. But if they hit reload they’d lose all the old ones
  • I gather google isn’t crazy about being an image hoster. It’s possible they’ll throttle you. In my (admittedly small scale) testing, I haven’t seen that. Famous last words, maybe
  • I used a cool trick to automatically resize the image and upload it after a student selects the image. There’s no separate “submit” button. I like that
  • On both iPhones and android when you select a new image you get the choice of your camera to take a picture. I love that as it reduces the steps for the students.

Your thoughts?

I’d love to hear what you think. Here are some starters for you:

  • I love this. What I particularly like is . . .
  • I can’t believe you reinvented the wheel, especially after such great suggestions from twitter. Loser
  • I know why it has to use an incognito window on the iPhone . . .
  • What’s wrong with cosine?
  • Google here: cease and desist
  • pusher here: are you ever going to pay for our service (note, I have paid when my school has scaled up some of my projects – they let you go month to month which is great)
  • Here’s another cool way you could use it . . .
  • I don’t think you need to be so worried about the original sized images and here’s why . . .
Posted in Google Apps Script, syllabus creation, teaching, technology | Tagged | Leave a comment

Audience ranking questions

I’m helping to run a workshop this week and we realized that instead of using post-it notes through the day to capture “burning questions” we could use a tool that both collects questions and allows the audience to vote on them.

After the first day of post-it notes I volunteered to get us set up with one of the many systems that exist that do this. I went back to my hotel room and did the usual google searching, finding tools I had used (including the one that the Physics Education Research Conference organizers used earlier this month at their opening session) and others I hadn’t heard of. Unfortunately it seemed that most had a free level that only worked with 7 participants but it got expensive after that. Also most of the tools offered all kinds of extra bells and whistles, like tying in content or finding other ways to foster community, that we didn’t need.

I knew I could use the Q&A feature built into Google Slides (don’t know about that? Check it out!) but I knew there were two major issues: 1) it’s not really built to be open all day. In my experimentation in the past I noticed that the URL changes when you go back to the Q&A window a few hours later. Really that’s because it’s meant to be used in a typical 1-2 hour lecture, rather than a place to keep questions for a longer time. 2) while anyone can submit questions, you have to be logged into your Google account in order to be able to vote on questions. I’m not sure that’s a huge deal, but it always gives me pause when I’m dealing with an audience from diverse places. Certainly since my school has all our email powered by Google, it would be my go-to choice. Really it’s the first issue above that caused me to go down a different road for this workshop.

So as I was taking a shower after that disappointing google searching, it occurred to me that I might have the skills necessary to write my own solution. This post is about that.

Let’s first see it in action (there’s no sound, just trying to show all the features):

What are you seeing?

  • A big box at the top to add your question.
    • It has a placeholder encouraging you to type in your question. I find this much cleaner than a label above the box.
    • It clears the box out after the submission has been accepted (typically ~2 seconds)
  • A ranked list of questions that gets updated in real time
    • If the question is from you, you can’t “upvote” it
    • You can only vote once on a question
    • You can’t downvote or unvote a question
  • A quick description at the bottom indicating the limitations

If you care, the rest of this post explains how you could build this yourself. Why not just use mine? Because of the limitations listed at the bottom of the vid: it can only handle 100 simultaneous users. That’s totally fine for this workshop, but I can’t just provide this tool to everyone. Instead, all you need is:

  1. A google account (you need one, your users don’t)
  2. A pusher.com free-level account (this enables the webSocket technology and has the 100 simultaneous user limitation)
  3. A burning desire to help your fellow human beings

Set up Pusher

  • Sign up for a pusher.com account, specifying the “sandbox” plan.
  • Under “Channels” create an app
You should probably choose a name that’s easier to remember than the default one they’ll prefill that field with
  • Hit “Create app”
  • Navigate to the new app and click on “App Keys”
You’ll need all 4 of these in your Google Apps Script below.
  • You’re done with pusher!
  • Now on to setting up your google apps script

Set up the spreadsheet

You only two tabs:

  • “questions”
    • id
    • question
    • date
  • “votes”
    • id
    • date

Then open the App script window:

You don’t have to call your google sheet “pusher q&a sheet”, of course

In the app script window you should have a code.gs file. Add a main.html file and a pusher.gs file:

You can name your script whatever you want. Also note that after you hit the “+” button it asks what type of file you want so you don’t have to add the extension.

Here’s the code for Code.gs:

var funcs=[];
function doGet()
{
  var ss=SpreadsheetApp.getActive();
  var sheet=ss.getSheetByName("questions");
  var data=sheet.getDataRange().getValues();
  data.shift();
  var vsheet=ss.getSheetByName("votes");
  var vdata=vsheet.getDataRange().getValues();
  vdata.shift();
  var vobj={};
  data.forEach(d=>vobj[d[0]]=0);
  vdata.forEach(v=>
  {
    // if(!vobj[v[0]]) vobj[v[0]]=0;
    vobj[v[0]]++
  })
  var arr=data.map(m=>[m[0],m[1],vobj[m[0]]]);


  
  var t=HtmlService.createTemplateFromFile("main");
  t.funcs=funcs;
  t.funcnames=t.funcs.map(f=>f.name);
  t.globals={data:data, votes:vobj, arr:arr, myvotes:[]};
  return t.evaluate().setTitle("Sorted Questions");
}

const init=()=>{
  var cook=document.cookie;
  var str=cook.match(/votes=(.*)/);
  if (!str)
  {
    var list=[];
  } else {
    var list=str[1].split(",").map(m=>Number(m));
  }
  myvotes=[...myvotes,...list];
  displayQs();
}
funcs.push(init)

const displayQs=()=>
{
  arr.sort((a,b)=>b[2]-a[2]);
  // html=`<h1>Sorted Questions</h1>`;
  // html+=`<div><textarea rows="4" cols="50" id="newQuestion" placeholder="enter your question"></textarea></div><div><button onclick="newQlocal()">submit</button></div><hr/>`;
  var html="";
  html+=arr.map(a=>
  {
    var h=`<p>(${a[2]} votes) ${a[1]}`;
    if (!myvotes.includes(a[0]))
    {
      h+=` <button onclick="vote(this, ${a[0]})">upvote</button>`;
    } else {
      h+=` You upvoted`
    }
    h+="</p>";
    return h;
  }).join(" ");
  document.getElementById("main").innerHTML=html;
}
funcs.push(displayQs);

const newQlocal=()=>
{
  var question=document.getElementById("newQuestion").value;
  google.script.run.withSuccessHandler(newQBack).sendNewQuestion(question);
}
funcs.push(newQlocal);

const newQBack=(id)=>
{
  myvotes.push(id);
  var str=myvotes.join(",");
  document.cookie = "votes="+str+"; SameSite=none; secure";
  document.getElementById("newQuestion").value="";
  displayQs();
}
funcs.push(newQBack);

const vote=(el,id)=>
{
  // myvotes.push(id);
  google.script.run.withSuccessHandler(newVBack).sendVote(id);
}
funcs.push(vote);

const newVBack=(id)=>
{
  myvotes.push(id);
  var str=myvotes.join(",");
  document.cookie = "votes="+str+"; SameSite=none; secure";
  displayQs();
}
funcs.push(newVBack);


function sendNewQuestion(question) 
{
  var sheet=SpreadsheetApp.getActive().getSheetByName("questions");
  var data=sheet.getDataRange().getValues();
  data.shift();
  if (data.length==0)
  {
    var newId=1;
  } else {
    var newId=Math.max(...data.map(m=>m[0]))+1;
  }
  Logger.log(`new id is ${newId}`)
  var d=new Date();
  sheet.appendRow([newId,question,d]);
  sendToPusher("newQ", {row: [newId,question,0]});
  return newId
}

function sendVote(id)
{
  var sheet=SpreadsheetApp.getActive().getSheetByName("votes");
  var d=new Date();
  sheet.appendRow([id,d]);
  sendToPusher("newV", {id:id});
  return id;
}

Here’s the html for main.html:

<!DOCTYPE html>
<html>
  <head>
    <base target="_top">
    <link href="https://cdn.jsdelivr.net/npm/bootstrap@5.0.2/dist/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-EVSTQN3/azprG1Anm3QDgpJLIm9Nao0Yz1ztcQTwFspd3yD65VohhpuuCOmLASjC" crossorigin="anonymous">

    <script src="https://js.pusher.com/5.1/pusher.min.js"></script>

  </head>
  <body onload="init()">
    <div class="container">
      <div>
         <h1>Sorted Questions</h1>
  <div><textarea rows="4" cols="50" id="newQuestion" placeholder="enter your question"></textarea></div><div><button onclick="newQlocal()">submit</button></div><hr/>
      </div>
      <div id="main">
      </div>
      <div>
        <hr/>
        <p>If your vote isn't registered after ~10 seconds, it means
        the system lost your vote. Feel free to try again. The system
        has a limit of 30 simultaneous votes and 100 connected users.</p>
    </div>

<script>
var globals = <?!= JSON.stringify(globals) ?>;
Object.keys(globals).forEach(key=>window[key]=globals[key]);
var funcnames=<?!= JSON.stringify(funcnames) ?>;
var funcs=[<?!= funcs ?>];
funcnames.forEach((fn,i)=>window[fn]=funcs[i]);

var pusher = new Pusher(key, {
    cluster: 'us3',
    forceTLS: true
  });

var channel = pusher.subscribe('my-channel');
channel.bind('newQ', function(data) {
  arr.push(data.row)
  displayQs();
  });
channel.bind('newV', function(data)
{
  var row=arr.find(f=>f[0]==data.id);
  row[2]++;
  displayQs();
})

</script>

  </body>
</html>

And here’s the code for pusher.gs:


var app_id = "YOUR APP_ID HERE";
var key = "YOUR KEY HERE";
var secret = "YOUR SECRET HERE";
var cluster = "YOUR CLUSTER HERE";


function sendToPusher(event,data) {
  var pvals={
    appId: app_id,
    key: key,
    secret: secret,
    cluster: cluster,
    encrypted: true
  };
  
  var url = `https://api-${pvals["cluster"]}.pusher.com/apps/${pvals["appId"]}/events`;
  var body = {"name":event,"channels":["my-channel"],"data":JSON.stringify(data)};
  var bodystring = JSON.stringify(body);
  var now=new Date();
  var d = Math.round(now.getTime() / 1000);
  var auth_timestamp = d;
  var auth_version = '1.0';
  var bodymd5 = byteToString(Utilities.computeDigest(Utilities.DigestAlgorithm.MD5, bodystring));
  var wholething = `POST
/apps/${pvals["appId"]}/events
auth_key=${pvals["key"]}&auth_timestamp=${auth_timestamp}&auth_version=${auth_version}&body_md5=${bodymd5}`;
  var wholethingencrypt = byteToString(Utilities.computeHmacSha256Signature(wholething,pvals["secret"]));
  Logger.log(wholethingencrypt);
  

  var options = {
    'method' : 'post',
    'contentType': 'application/json',
    // Convert the JavaScript object to a JSON string.
    'payload' : bodystring,
    'muteHttpExceptions' : true
  };
  var urltry = UrlFetchApp.fetch(url+`?auth_key=${pvals["key"]}&auth_timestamp=${auth_timestamp}&auth_version=${auth_version}&body_md5=${bodymd5}&auth_signature=${wholethingencrypt}`, options);

  
  }
  
function byteToString(byte) {
  var signature = byte.reduce(function(str,chr){
    chr = (chr < 0 ? chr + 256 : chr).toString(16);
    return str + (chr.length==1?'0':'') + chr;
  },'');
  return signature;
  }

Make sure you make the changes in the top 4 lines before moving forward.

Next you have to create a web app. Make sure to set it to execute as you but be available to everyone.

The first time you should choose “new deployment” and choose “web app”

Now you can go to the url provided and it should be working!

How it works

When someone goes to the url, google sends them all the questions along with the vote tally for each. It also checks the local users cookies to see if they’ve supplied any votes. If they have, it doesn’t let them vote again. This prevents ballot box stuffing, but you should know that it’s pretty easily defeated using incognito windows.

The user can enter questions or “upvote” existing questions. When they do, an AJAX call is made to the google server (any time you see google.script.run… that’s what’s happening) where google either saves the new question (making sure to give it a unique id) or saves the vote. In either case it saves the timestamp.

After updating the spreadsheet, the “sendtopusher” function runs, sending along either the new question (along with its id and an initial vote count of zero) or the id of the new vote. That uses the webSocket that pusher has set up to send that info to all connected devices.

If a device receives a new question from pusher it adds it to the list and re-displays all the questions. If a new vote comes in from pusher it adds to the vote count for that id and re-displays all the questions (this also involves sorting based on the vote count so that the highest voted question is always at the top).

Your thoughts?

I’d love to hear your thoughts. Here are some starters for you:

  • Why do you sometimes capitalize google and other times you don’t capitalize Google?
  • This is cool, can it also . . .
  • This is a rip-off of my cool idea. You can send checks to . . .
  • This is dumb. How does this preserve the time honored tradition of the first person asking not so much a question as a 10 point rebuff of everything they heard?
  • You really can’t draw very straight arrows. It’s almost as if you’re writing this post on a chromebook on your lap at LAX
  • I can tell you made that vid on Loom. Why didn’t you just embed that instead of downloading it from loom, posting it on youtube, and then embedding that?
  • I was at this workshop and I found this very useful.
  • I was at this workshop and this was incredibly distracting
  • Gross, I didn’t need to know that you had this idea in the shower.
  • Here’s another way to protect against ballot box stuffing . . .
Posted in Google Apps Script, syllabus creation | Tagged | 2 Comments

Brachistochrone for rolling things

The Brachistochrone curve is the shape of a wire for beads to slide down (friction free) to get from point A to point B the fastest. Note that since I used the word “down” there I’m implying this happens in gravity. Here’s an old post of mine describing how I go about teaching it. This post is all about scratching an itch I’ve had for a while: What if instead of sliding beads we want to roll balls. Is the shape the same? Spoiler: Nope, not the same.

My first thoughts had to do with how you’d factor rolling into the typical analysis. Normally you determine the integral formula for the time to go from A to B on an arbitrary curve given by y(x):

\text{time}=\int_A^B dt=\int_A^B \frac{ds}{v(s)}=\int_A^B\frac{\sqrt{1+y'^2}dx}{\sqrt{2\text{KE}/m}}

where y’ is the slope of the curve, s is how far along the curve the bead has gone, v is how fast the bead is traveling, and KE is the kinetic energy which is usually a function of y (since you’re cashing in gravitational potential energy). So if it’s rolling without slipping, my first thought was that all I had to do was add in some rotational kinetic energy:

\text{KE}=\frac{1}{2}m v^2+\frac{1}{2}I\omega^2

But then I realized that I had to know exactly where the center of mass was in order to figure out how much potential energy had been cashed in and I went down a rabbit hole.

Does the center of mass follow the same curve? (No!)

You’ll see that at this point I jumped on twitter for the first time in a while (hey, my job is different now and a lot of what I do is untweetable, give me a break).

If you scroll through you’ll see some curves that I no longer stand behind

If you have a curvy road and you know the mathematical formula for one side (let’s say the road is going left to right along what we’ll call the “x-axis” and that it doesn’t turn back on itself so we can call it a function). Do you know the formula for the other side of the road? Is it just the same function with a shift? Nope. It took me a while to convince myself but this is the figure that sold me:

Blue: right side of the road. Orange: right side of the road with a constant added. Green: the left side of the road

The blue curve is a pure sine function (why would I ever use cosine?). The orange curve is something like “sin(x)+0.32”. The green curve is what took me a while to derive but it’s really what enables a 0.32 diameter ball to fit between green and blue everywhere. Note that green and orange have the same amplitude and same frequency. Therefore, since they don’t overlap, the green curve is not a sinusoid.

So how do you derive the green curve? Well, here’s how I did it:

This represents a zoomed in version where locally the curve is flat. Also note that if you’re rolling up the other side (when the slope is positive) you need to make some adjustments to the signs in those equations. But that’s really all it takes. If you have a function for y(x) and you can calculate its slope at every location (y’), then you can figure out where the center of mass of the ball will be when you know the contact point with the curve.

Obviously I coded that in and ran it for a sine curve to get the figure above, but my same code would work with any (differentiable) function. Note that if the curve has a constant slope, the adjustments for the center of mass location are constant and then the other side of the road is truly just a shifted function. But that’s the only case that leads to that simple conclusion.

Alright! Let’s solve a complicated differential equation!

Ok, so now we know how to find the center of mass location when you know the contact point. It seems like we could figure out the potential energy drop (and hence the kinetic energy) since we know the vertical drop of the ball. Seems like we’d be in business! Alas, no, we’re not. The problem is the angular frequency, or \omega in the equation above.

For rolling without slipping on a flat surface, you know that your linear speed and rotational speed are tied together, namely \omega=v/r. Unfortunately, that’s not the case when rolling on a curved surface. This web page helped me understand this a little better. When you have a curved surface that has a local radius of curvature, \rho you get this for \omega:

\omega=\frac{\rho-R}{\rho R}v

where v is the speed of the contact point along the surface.

No big deal, right? it’s just some weird multiplier in front of the speed. That should make solving for the speed from the kinetic energy easy! Well, that’s what I thought, and certainly that’s what led to me erroneous twitter posts (if you scrolled through). Unfortunately, \rho, you know that pesky local radius of curvature, is not easy to deal with. From wikipedia I learned that:

\rho=\frac{\left(1+y'^2\right)^{3/2}}{y''}

Ugh! Do you see that denominator?! Suddenly you need to know not just the slope of the function but its curvature as well. Let me tell you, that makes things gross.

Ok, gross maybe we can handle. We know how to calculate the kinetic energy and it’ll be all in terms of the (unknown) function, its slope, and its curvature. Maybe we can just close our eyes and throw it to Mathematica. Here’s where we’re at:

\text{time}=\int_A^B \frac{\sqrt{1+y'^2}}{\frac{\left(1+y'^2\right)^{3/2}R/y''}{\left(1+y'^2\right)^{3/2}/y''-R}\sqrt{\frac{M g \left(y_0-\left(y-\frac{R/y'}{\sqrt{1+1/y'^2}}\right)\right)}{\frac{1}{2}(I+MR^2)}}}dx

Fun right! Anyways, it’s technically all set to use the calculus of variations, but I’ve tried it, and wasn’t able to make any progress. 😦 I think the biggest problem is the y”s in there because they lead to a third order differential equation, which means I need to supply not only where to start the curve and what direction to head, but also the local curvature right there. Needless to say, I didn’t make much progress. If you have ideas, I’m all ears!

By the way, here’s what it looks like if you’re just doing a bead sliding down a wire:

\text{time}=\int_A^B\frac{\sqrt{1+y'^2}}{2 g (y_0-y)}dx

Muuuccch easier, trust me. (Also note that if you thought I’d be using the word “cycloid” by now, you don’t get there this way. You only do if you swap x and y. You know an “obvious” thing surely your students would think to do.)

When in doubt, check the literature

So I started googling. Here’s an awesome paper from 1946 that helps us put it all together. What they’re saying is that even when rolling on a curved surface, you can use \omega=v/r as long as you’re using the speed of the center of mass, not the speed of the contact point. Alas, even though they’re always moving in parallel, they don’t have the same speed (think about going up and over a hump in a roller coaster, you’re moving faster than the contact point on the track). Note that they’re also saying that the center of mass follows the traditional brachistochrone! So what is this post all about!? Well, we want to know the shape of the track the ball is rolling on, and if you’ve read what I wrote above you’d know that’s different!

How did they prove it was the traditional curve? Because you get the very simple equation above instead of the incredibly ugly one if you use the coordinates of the center of mass and not the contact point. With that same simple equation, you get the same simple result (if you must: a cycloid).

But now we can put it all together. If I have a normal brachistochrone, I can find the curve for the ball to roll on by doing the coordinate shift in the figure above in reverse!

Blue is the shape of the track for the ball to roll on and red is the path of the center of mass

I know, I know, the blue path (the track for the ball to roll on) sure looks like a standard brachistochrone, but it’s not, because of what I was talking about above. Don’t believe me, let me hear it!

Update!

I don’t know why I didn’t do this last night, but here’s that same image with an added brachistochrone from the start to the finish of the track in green. See, I told you the blue curve wasn’t a brachistochrone:

Red: brachistochrone that the center of mass follows for a ball rolling without slipping on the blue track. Blue: The correct shape of a track for a ball with radius 0.1 to roll from the start to the finish the fastest. Green: the correct brachistochrone for a bead to slide on from the beginning of the track to the end of the track the fastest.

Starters

Your thoughts? Here are some starters for you:

  • What do you mean all you had to do was say “down” to imply gravity?
  • Seriously, I have to read a whole other post of yours just to be able to read this one. No way! I’m unclicking. You can’t count my click.
  • What do you have against rabbits? Why does going down their holes feel like an interminable complicated journey?
  • What do you mean about your job being untweetable?
  • What do you have against cosine?
  • Duh, of course you needed to know about the curvature. What are you, an idiot?
  • I know exactly where you made a mistake in that big ugly equation. For $5 I’ll tell you.
  • Of course switching x and y is obvious. What’s you’re point?
  • Hang on, this was solved back in 1946 and I had to read nearly your whole post to get there? Jerk.
  • That blue curve is a brachistochrone and I’ve blogged about this a bunch. Try reading some time.
Posted in math, mathematica, physics, teaching | 3 Comments