Real life ghost mode

This post describes how I figured out how to ride to work while racing earlier versions of myself. Think of it like Mario Kart’s Ghost mode, only without cool 3D virtual reality and with a phone screen that shows your position and the position of 3 random rides along the same route.

1. What does it look like?
2. Don’t care about the details, just tell me how to do it
3. What I learned (or: “I care more about your journey of learning than just taking what you’ve done and implementing it for myself”)
   1. .tcx files
   2. Maps
   3. Drawing on html canvas elements
   4. Animation on a web page
4. Gratuitous gifs
5. Your thoughts?

What does it look like?

Shows poor phone GPS (actually on the bridge)

These are all taken from one of my rides. The first one shows some poor phone GPS because I was actually on that big bridge. The rest were are different points on a ride I was pretty proud of, since I won! (note that I also came in second, third, and fourth!). The gray rectangle at the bottom is a mismatch between my html (see below) and my phone screen size that I haven’t bothered to fix yet.

The screen updates every time the phone gets a new GPS location. For my phone (Google Pixel 7 with Verizon) it’s about every 8 seconds. The map is set to keep all four markers in the frame, so it constantly adjusts the pan and zoom to do that.

Don’t care about the details, just tell me how to do it

I’ll put the details of what I learned below. Here’s how to do it yourself:

1. Get yourself a device that can track your workout that allows you to export .tcx files. That’s what fitbit uses, but I gather it was invented by Garmin so I think there’s lots of ways to do that.
2. Create a folder in Google Drive that will hold all the tcx files (make a different folder for different types of rides).
3. Create a new Google Apps Script file.
4. In code.gs, put in this:

var funcs=[];
var allData=[];
function doGet(e) {
  var t=HtmlService.createTemplateFromFile("main");
  t.funcs=funcs;
  t.funcnames=t.funcs.map(f=>f.name);
  prepData(e.parameter.dir);
  var timeTotal=Math.max(...allData.map(m=>m.times).flat());
  t.globals={allData:allData,
             currentTime:0,
             timeTotal:timeTotal,
             colors:["red","blue","green"],
             initTime:0,
             map:null,
             markers:[],
             locationMarker:null,
             locationCircle:null,
             group:null,
            };
  return t.evaluate();
}

const getLocations=(time)=>
{
  var locations=[];
  allData.forEach(a=>
  {
    var targetIndex=a.times.findIndex(f=>f>=time);
    if(targetIndex==-1) targetIndex=a.times.length-1;
    locations.push([a.lats[targetIndex], a.longs[targetIndex]]);
  })
  return locations;
}
funcs.push(getLocations);

function prepData(dir="to")
{
  var folder=DriveApp.getFolderById(dir=="to"?"[FOLDERID FOR THE 'to' RIDES]":"[FOLDERID FOR THE 'from' RIDES]");
  var files=folder.getFiles();
  var fileIds=[];
  while(files.hasNext())
  {
    var file=files.next();
    fileIds.push(file.getId());
  }
  var upToThree=[];
  while(upToThree.length<3)
  {
    var id=fileIds[Math.floor(Math.random()*fileIds.length)];
    if(!upToThree.includes(id)) upToThree.push(id);
  }
  allData=upToThree.map(id=>readTcxFile(id));
}

function readTcxFile(id)
{
  var file=DriveApp.getFileById(id);
  var text=file.getBlob().getDataAsString();
  // Logger.log(text.slice(0,10));
  var times=[...text.matchAll(/<Time>(.*)<\/Time>/g)].map(m=>new Date(m[1]).getTime());
  var fT=times[0];
  times=times.map(m=>m-fT);
  var lats=[...text.matchAll(/<LatitudeDegrees>(.*)<\/LatitudeDegrees>/g)].map(m=>m[1]);
  var longs=[...text.matchAll(/<LongitudeDegrees>(.*)<\/LongitudeDegrees>/g)].map(m=>m[1]);
  var alts=[...text.matchAll(/<AltitudeMeters>(.*)<\/AltitudeMeters>/g)].map(m=>m[1]);
  return {times:times, lats:lats, longs:longs, alts:alts};
  
}

5. Change the “[FOLDERID FOR THE ‘to’ RIDES]” to the folder id for your “to” rides and similar for your “from” rides.
   • note that you could change the logic on line 38 to do any number of different folders. I mostly use this for riding to and from work so I use this logic.
6. Make a new file in the script and call it main.html

7. Put this code in the “main.html” file:

<!DOCTYPE html>
<html lang="en">
<head>
	<base target="_top">
	<meta charset="utf-8">
	<meta name="viewport" content="width=device-width, initial-scale=1">
	
	<title>Mobile tutorial - Leaflet</title>
	
	<link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css" integrity="sha256-p4NxAoJBhIIN+hmNHrzRCf9tD/miZyoHS5obTRR9BMY=" crossorigin=""/>
    <script src="https://unpkg.com/leaflet@1.9.4/dist/leaflet.js" integrity="sha256-20nQCchB9co0qIjJZRGuk2/Z9VM+kNiyxNV1lvTlZBo=" crossorigin=""></script>

	<style>
		html, body {
			height: 100%;
			margin: 0;
		}
		.leaflet-container {
			height: 400px;
			width: 600px;
			max-width: 100%;
			max-height: 100%;
		}
	</style>

	<style>body { padding: 0; margin: 0; } #map { height: 100%; width: 100vw; }</style>
</head>
<body onload="init()">

<div id='map'></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]);
	function init()
  {
    map = L.map('map');
  map.on('load', ()=>
  {
    initTime=new Date().getTime();
    console.log(`initTime onload=${initTime}`);
  })
  map.setView([45,-90],11);
	const tiles = L.tileLayer('https://tile.openstreetmap.org/{z}/{x}/{y}.png', {
		maxZoom: 19,
		attribution: '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>'
	}).addTo(map);

  allData.forEach((d,j)=>
	{
		var latlngs=d.lats.map((m,i)=>[m,d.longs[i]]);
		L.polyline(latlngs, {color: colors[j]}).addTo(map);
		var marker=L.marker([d.lats[0], d.longs[0]]);
		marker.addTo(map);
		markers.push(marker);
	});

  locationMarker=L.marker([0,0]);
  locationMarker.addTo(map);
  locationCircle=L.circle([0,0],10);
  locationCircle.addTo(map);
  group = new L.featureGroup([locationMarker, ...markers]);

	

	map.on('locationfound', onLocationFound);
	map.on('locationerror', onLocationError);

	// map.locate({setView: true, maxZoom: 16});
  window.requestAnimationFrame(draw)
  }
  function draw()
  {
    map.locate();
    window.requestAnimationFrame(draw)
  }

  function onLocationFound(e) {
		const radius = e.accuracy / 2;
    currentTime=new Date().getTime();
    // console.log(`initTime=${initTime}, currentTime=${currentTime}, diff=${currentTime-initTime}`);
    var newlocs=getLocations(currentTime-initTime);
    console.log(newlocs);
    newlocs.forEach((locs,i)=>
    {
      markers[i].setLatLng(locs);
    })
		// const locationMarker = L.marker(e.latlng).addTo(map)
		// 	.bindPopup(`You are within ${radius} meters from this point`).openPopup();
    locationMarker.setLatLng(e.latlng);
    locationCircle.setRadius(radius);
    locationCircle.setLatLng(e.latlng);
    

    map.fitBounds(group.getBounds());

		// const locationCircle = L.circle(e.latlng, radius).addTo(map);
	}

	function onLocationError(e) {
		alert(e.message);
	}
</script>



</body>
</html>

8. Set up a new web app deployment with this script and you’re good to go! (note that you need to run the doGet function once just to get all the permissions set).
   • To use the logic on line 38 of the code.gs, add “?dir=to” or “?dir=from” to the url you get from the deployment.

What I learned (or: “I care more about your journey of learning than just taking what you’ve done and implementing it for myself”)

Here’s a quick vid showing what I was originally excited to build. It shows what it does and walks through how to do the code:

Below I’ll talk about all the interesting technical things I learned how to do:

.tcx files

With my new Pixel Watch (technically the Pixel Watch 2 if you must know), you get fitbit for free (though they strongly encourage you to move up to premium which I haven’t done). Tracking a ride is really easy. Just wake up the watch, swipe once, then hit “bike” and it’s tracking. At the end of the day you just have to export the tcx file wherever you want it:

Here’s the top few lines of one of my rides:

<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<TrainingCenterDatabase xmlns="http://www.garmin.com/xmlschemas/TrainingCenterDatabase/v2">
    <Activities>
        <Activity Sport="Biking">
            <Id>2024-07-15T16:37:48.000-05:00</Id>
            <Lap StartTime="2024-07-15T16:37:48.000-05:00">
                <TotalTimeSeconds>1876.0</TotalTimeSeconds>
                <DistanceMeters>10604.77</DistanceMeters>
                <Calories>321</Calories>
                <Intensity>Active</Intensity>
                <TriggerMethod>Manual</TriggerMethod>
                <Track>
                    <Trackpoint>
                        <Time>2024-07-15T16:37:48.000-05:00</Time>
                        <Position>
                            <LatitudeDegrees>44.96538329124451</LatitudeDegrees>
                            <LongitudeDegrees>-93.16418159008026</LongitudeDegrees>
                        </Position>
                        <AltitudeMeters>286.8</AltitudeMeters>
                        <DistanceMeters>0.0</DistanceMeters>
                        <HeartRateBpm>
                            <Value>83</Value>
                        </HeartRateBpm>
                    </Trackpoint>
                    <Trackpoint>
                        <Time>2024-07-15T16:37:53.000-05:00</Time>
                        <Position>
                            <LatitudeDegrees>44.96538329124451</LatitudeDegrees>
                            <LongitudeDegrees>-93.16418159008026</LongitudeDegrees>
                        </Position>
                        <AltitudeMeters>286.8</AltitudeMeters>
                        <DistanceMeters>0.0</DistanceMeters>
                        <HeartRateBpm>
                            <Value>85</Value>
                        </HeartRateBpm>
                    </Trackpoint>

Typically it seems I get a new “Trackpoint” every second. So most of what I did had to do with scraping all of the Times, Latitudes, Longitudes, Altitudes, Distances, and HeartRates. Typically I do that using some regex. In Google Apps Script you can see an example on lines 63-65 in the code up above. I just grab them all and store them in javascript arrays so that they’re all indexed the same (if I’m looking at the 51st time, then the 51st element of the “lats” array is the latitude at that time).

Maps

As you can likely tell, I ultimately went with Leaflet.js for controlling the maps on the web pages and with openstreetmaps for what are called the tiles that Leaflet needs to show. I started doing a ton with the Google Maps Static API but they’re static and not supposed to be used with Leaflet. Also I get to make 1000 Google static maps per day using Google Apps Script but to connect it with Leaflet would mean using their normal API and that starts costing money much quicker. OpenStreetMaps, at least for small projects like this, is a perfect choice, it seems to me.

Brief aside about how to choose the appropriate initial zoom when you know all your latitudes and longitudes: It turns out that all the mapping things I’ve played with for this project use the same calculation for zoom level. I was interested in making sure that I picked the integer zoom level that for sure contained all my position data. So I’d make sure to know the extremes (most northerly, southernly, westerly, and easterly) and then I needed to figure out what zoom level for sure contained that. At first I was constantly making 600×400 maps so I just needed to make sure that everything would fit. It turns out that the key is to know that at zero zoom level, the full equator should fit in 256 pixels. Every zoom after that leads to an additional factor of two. I knew I wanted all my horizontal stuff to fit in 600 pixels and all my vertical stuff to fit in 400 pixels. I’d figure out the decimal zoom that would achieve those (two different answers often) and take the most zoomed out one as the one that would capture everything. Then I would round that down to the nearest integer and use that for the initial zoom level. You can see that in the video above.

Once you tell Leaflet to use openstreetmaps, you can then load the map with additional markers and lines. What you see in all the things above are the lines of all the routes shown in the race and markers for the current location. Here’s some pseudo code that makes a map and adds some lines and markers:

// put this all in the <head> of your html:
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css" integrity="sha256-p4NxAoJBhIIN+hmNHrzRCf9tD/miZyoHS5obTRR9BMY=" crossorigin=""/>
    <script src="https://unpkg.com/leaflet@1.9.4/dist/leaflet.js" integrity="sha256-20nQCchB9co0qIjJZRGuk2/Z9VM+kNiyxNV1lvTlZBo=" crossorigin=""></script>

// then put this in your script:
var map=L.map('map').setView([45,-90],11) // initial center and initial zoom level. 
// 'map' is the id of the <div> you need to have in your <body> to place the map
const tiles = L.tileLayer('https://tile.openstreetmap.org/{z}/{x}/{y}.png', {
		maxZoom: 19,
		attribution: '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>'
	}).addTo(map); // this actually gets the "tiles" to draw
var marker=L.marker[45,-90])
marker.addTo(map) // actually puts it on the map but it's also a variable you can adjust

var listOfLocations=[[45,-90], [45.001, -90], [45.002,-90]];
L.polyline(listOfLocations).addTo(map) // draws a boring line

Drawing on html canvas elements

At first I thought I had to do a bunch of external drawing (basically an html canvas element on top of the map) but Leaflet let’s me do all of that just like my pseudo-code above. However, for the speed, heart rate, altitude, and distance plots you see in the youtube vid above I still had to figure out how to make those (animated!) plots.

Because they all have all the data constantly and a moving set of markers, I actually use two canvas elements on top of each other. The lower one has the full plots and the top one is constantly redrawn (see the animation section below) with the moving markers.

To have stacked canvas elements, just put them in the same div and make sure their css uses absolute location:

<!-- put this in the head somewhere -->
<style>
    .wrapper {
    position: relative;
    width: 600px;
    height: 400px;
}

.wrapper canvas {
    position: absolute;
    top: 0;
    left: 0;
}
</style>

<!-- put this in the body somewhere -->

<div id="speedDiv" class="wrapper col">
		<canvas id="speedsCanvas" width="300" height="200"></canvas>
		<canvas id="speedsCanvasDot" width="300" height="200"></canvas>
</div>

Here’s the code to draw a simple line to the bottom canvas and a marker on the top canvas:

var bottom = document.getElementById("speedsCanvas");
var top = document.getElementById("speedsCanvasDot");
var context = bottom.getContext('2d');
context.clearRect(0,0,canvas.width,canvas.height);
context.moveTo(100,50); // where to start the line
context.lineTo(110,60); // draw a line from the start to here
context.lineTo(120, 70); // draw a line from where you are now to here
context.stroke(); // actually draw the line

context=top.getContext('2d'); // now grab the other canvas and draw a marker
context.fillStyle="red";
context.fillRect(100,50,10); // at 100 from the left, 50 from the top, and 10 pixels wide

Here’s the part that sucks: the origin (of the pixels) is at the upper left of the canvas. So your horizontal instincts are all correct, but your vertical ones are all backwards.

Animation on a web page

For any tcx files that I’m trying to animate, I reset all of their times to start at zero. I also convert all time stamps to the number of milliseconds since 1/1/1970 using new Date(datestamp).getTime(). Then, when I’m animating things, I just have javascript constantly recheck the time, determine how much time has gone by since the animation started, and then choose the appropriate latitudes and longitudes to display on the map. All of this is done with the magical window.requestAnimationFrame().

Actually, for the ghost mode app, I don’t have it constantly recheck the time. Instead, I wait until there’s a new GPS measurement, and then I check the time. On my phone that’s about every eight seconds. There is vanilla javascript to check the gps location, but leaflet has it built in with the command map.locate(). It’s that command that I repeat over and over again, waiting for the results before updating the map. Here’s the pseudocode for that:

// globals
var map, locationMarker, group

map = L.map('map');

// this next part makes sure to set the initial time once the map is loaded
map.on('load', ()=>
     {
       initTime=new Date().getTime();
     })
map.setView([45,-90],11);
const tiles = L.tileLayer('https://tile.openstreetmap.org/{z}/{x}/{y}.png', {
		maxZoom: 19,
		attribution: '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>'
	}).addTo(map);

locationMarker=L.marker([0,0]);
locationMarker.addTo(map); // do this in two steps so it's a variable you can adjust
group=[locationMarker]; // sets up an easy way to pan/zoom to center you on the map

map.on('locationfound', onLocationFound); // set up what to do when a location is found
map.on('locationerror', onLocationError); // set up what to do if the location doesn't work

window.requestAnimationFrame(draw) // start the animation!

function draw() // function to call with every new animation frame request
  {
    map.locate(); // grabs the gps
    window.requestAnimationFrame(draw) // go to next frame (but map.locate will take a while)
  }

function onLocationFound(e) {
    locationMarker.setLatLng(e.latlng); // moves the marker to your location
    

    map.fitBounds(group.getBounds()); // cool way to center you
	}
function onLocationError(e) {
		alert(e.message);
	}

Gratuitous gifs

Here’s a few fun gifs showing the races I’ve had so far. Here’s all the rides to work I’ve done since I got my watch:

Here’s all the rides home:

And here’s just the 2 fastest so you can see the interesting effects of stop lights:

Your thoughts?

Oof, that’s a lot! Hopefully you found it interesting. I’d love to hear from you. Here are some starter questions for you to consider:

• This is cool! What happens when you get so old that you can’t beat your old times any more?
• This is dumb! _______ does the same thing a lot better! They’re going to sue you
• This is cool! I don’t understand what you mean when you say leaflet and google static maps don’t work together
• This is dumb! I tried it and it doesn’t work. You suck
• This is cool! You should have the bearing always be up (apparently you can’t rotate the tiles in leaflet, sorry)
• This is dumb! I don’t always start my tracking at exactly the same point on my driveway so I can’t be sure I’ve won!
• This is cool! Could I use it with my car to find which way to something much more fun than work is best?
• This is dumb! I don’t think you should be tracking any of this, let alone letting the government get their hands on it!
• This is cool! Could you add _____? It would be so awesome!
• This is dumb! Until you remove this feature _______ it really sucks!
• I know why you have a gray box on the bottom of your phone. Here’s the quick fix.

AppSheet ideas

AppSheet is a google product that enables you to use Google Sheets as a backend database for a mobile app that works in both the Android and Apple ecosystems. I encourage you to check it out, as it’s become a go-to tool for me to think about interacting with students and my colleagues at my school. Quick note: it’s really designed to be used by people on the same domain (yourcoolschool.edu, for example), not really to produce a true mobile app for the world. This post is about some brainstorming I’ve been doing to think about a campus app that could be useful.

What I’ve done so far

Here’s a quick list of the apps I’ve built and some of the cool features they make use of:

• Hamline Go
  • This is an app I built to try to help build community on campus. Students can see faculty and staff avatars on a map and collect them by clicking on them. They get one point for every collection unless its a faculty member who teaches in their major, then it’s 5 points. If instead they actually talk to a faculty or staff member and get their hourly code, they get 10 points. Certain offices on campus also have hourly codes. I also built in flash (1 hour), daily, and weekly challenges like “write a haiku about your major” or “prove you found this spot on campus”. Point leaders after the end of each month get declining balance (money) prizes.
  • Tools used:
    • Google Maps
    • Login features (checks user email and crosschecks the faculty in their declared major)
    • Editable images (you can draw on images before submitting them for the challenges)
    • Google Apps Script running independently to update both the avatar locations (randomized on our campus) and the hourly codes
• Majors t-shirts
  • Every major program was encouraged to design a tshirt that we would then order for any declared major that wanted one. I built an app to show students the designs available to them (only the designs for their declared majors) and let them submit the size that they wanted.
  • Tools used:
    • Login (check their major and crosscheck against the submitted designs)
    • Selfie collection (students are encouraged to submit a selfie wearing their tshirts so we can make a collage)
• Scavenger hunt
  • I made an app that encourages first year students to learn about the various offices on campus. Each campus had a QR code displayed that, if the students scanned it using the app, indicated that the student had visited that office
  • Tools used
    • QR code reader
    • Login

I’ve made it quite a bit up the learning curve, including some of these really useful tips and tricks:

• Make a slice of the users data called “me”. Have it be a filter of all the users that match the useremail() (a built in function). Essentially it’s a data source with just one row in it. Then I use index(me[id], 1) to get the current user’s id when doing a lot of the crosschecking mentioned above. I do the same thing on any of the related collections. In Hamline Go, for example, I make a “my avatars” slice to show the avatars that have been collected by that user. Then I can send that slice as the core data to a “view” that shows avatars.
• Use an automation to send a notification whenever there’s a new flash challenge. They only last an hour, so it’s helpful to let the students know there’s a new one. These notifications come as push notifications to the students phones.
• Pre-fill forms using LinkToForm that can take a collection of values to pre-fill elements. You can also then not show those elements in the form view and so it looks magical!
• Send email notifications to people that include live forms right in gmail. This works really great for approvals.

What I’m hoping to do

I’d love to create an app that I’m currently calling “Hamline All Around”. Here are some of the features I hope to build in:

• Somehow provide a navigation aid through our various systems. We get a lot of feedback from students who say they just don’t understand how to accomplish various things, especially the ones that are supposed to be providing valuable resources to them. A great example is our emergency grant program.
  • One thing I’m thinking of is a dynamic FAQ that might ask what sort of thing they’d like to accomplish and then quickly get them to the office they should start with.
• “Find a study buddy!” I gather other schools do this and so I built a quick test to see if this would work. The student logs in and it shows them the classes they’re taking. For each they can “raise their hand.” If they do, they see all the other students who’ve done the same, all of whom are indicating that they’re interested in finding study buddies. They can see each other’s emails, or potentially use the app’s push notification system to communicate. I also made a form that lets them say “I plan to study tomorrow night in the library from 7-8” and the others who’ve raised their hands can indicate if they plan to come. This is displayed in a nice calendar format that looks a lot like google calendar.
• Show the student group and athletic team calendar events in one place (right now you have to separately subscribe to both ical streams, but appsheet can do the subscription and display them in the app)
• Event check in: All employees and students have bar codes on the back of their ids and appsheet can read them (it’s a setting for one of the built-in form elements). It would be very easy to build in an event check-in system where the event runners (multiple of them) could have the app running on their phones and they could scan in the students/staff/faculty who attend. Once they do, the app on the students/staff/faculty phones could become much more functional (ie the app, upon being updated, would detect that the person is checked in and would then show all the relevant functionality for the event).
• Directions: Often it’s not enough to tell people what building things are in. Directions for how to actually get to the office can come in handy. I’ve been wondering about vids showing someone literally walking from a common place to an office, all the while describing what useful things the office can do for students.

Help!

I’m excited to work on this, but I’d love some more brainstorming, including poking some of these ideas a little. Here are some starters for you:

• I love this. What I especially love is . . .
• This is dumb. What I especially dislike is . . .
• Wait, you’re tracking the location of all faculty and staff in real time?! (note: it took me a long time to convince people I wasn’t doing that. I call the things avatars and try to make it clear that I just randomly place them on the campus map. However, several faculty and staff were deeply suspicious of this app and asked to be taken out of the app.)
• I think the “find a study buddy” app is really interesting. Have you thought about . . .
• I think the “find a study buddy” app is really disturbing. Instead you should . . .
• So after you figure out all the views and the data connections you have to learn how to write it for both the Android and iOS systems?! (no! Appsheet does all that for me. All the students have to install is Appsheet which already exists in both ecosystems)
• We’ve got a student app on our campus and I love it. Here are my favorite parts . . .
• We’ve got a student app on our campus and I hate it. What would make it so much better is . . .
• Here’s some additional functionality you should think of . . .
• Can you please do some tutorials on how to do all this in Appsheet?
• What’s an example of a haiku about a major?
  • Minnesota Law
  • Ins and Outs of Policing
  • Criminal Justice

Constraint forces for path in 3D

This is a post mostly for my future self. What got me thinking about it was this comment about my planetary tunnels post:

Straight tunnels? What about ones using the Lagrangian drooping rope optimized as the line?

— Andrew Jagger (@ajagger01) February 27, 2023

which is actually referring to a really old post of mine about how to dig a well.

In order to add in the effects of a spinning planet, there are two or three major routes I could take:

1. Set the constraints on the path to be for a spinning path.
2. Just add the coriolis and centrifugal forces to my other calculations
3. Do either (1) or (2) with some Lagrange Multipliers so that we could interrogate the constraint forces keeping the ball on in the tunnel

This post is really about (3), though I’ve done (2) and it’s pretty cool. I just put in the fictitious forces, but only those components that are along the tunnel (using a dot product). Here’s the result:

What constraints do we use for the Lagrange Multiplier approach?

When you do Lagrange Multipliers, which you usually do because you want to learn something about the forces that enforce the constraint, you need to have a formula of something that (typically) stays constant for all the dynamics.

I first thought I could do this with a single constraint. I figured if we have points a and b on the line and we want to find a measure of how far a third point, p, is form the line, we’d get a distance and the constraint would be that the distance should always be zero. I can get that distance by asking for the length of the vector given by crossing a normalized vector on the line with the vector from the point of interest, p, to either end point. As I tell my students, when you use a cross product, one vector is saying to the other: “hey you, how much of you is perpendicular to me?” and that perpendicular distance is exactly what I’m looking for. So, I figured this would work for a constraint (it should be zero for any point on the line):

So I tried it! Holy crow, it’s slow in the Wolfram Engine (Mathematica for free). Actually, if I start the simulation on the line, the WE just gives up:

r[t_]:={x[t],y[t],z[t]};
constraint=Cross[r[t]-p1,p2-p1].Cross[r[t]-p1,p2-p1];
KE=1/2 r'[t].r'[t];
PE=4/3 Pi r[t].r[t];
L=KE-PE;
el[a_]:=D[L,a[t]]-D[L,a'[t],t]+lm[t] D[constraint,a[t]]==0;
sol=First[NDSolve[{el/@{x,y,z},
                   x[0]==p1[[1]],
                   y[0]==p1[[2]],
                   z[0]==p1[[3]],
                   x'[0]==0,
                   y'[0]==0,
                   z'[0]==0,
                   D[constraint,t,t]==0},{x,y,z,lm},{t,0,1}]]

NDSolve::ndcf: Repeated convergence test failure at t == 0.; unable to continue.

If I start it just off the line, it meanders up and down the line, maintaining its distance. But, wow, it’s a really slow integration (note: all I did was change the x[0] value to p1[[1]]+0.01 on line 8)

Here’s a plot of the constraint function (note that it stays pretty constant):

Since it was so slow, I wanted to find another way. That’s when I realized that a straight line is given by the intersection of two planes. So maybe I could just say “hey, it’s on this plane and on that plane at all times.” That’s equivalent to saying that we have two constraints, not just the one like above.

So how do you find the equations for a plane (or two for that matter) that the original line is on? Well, the formula for a plane is so we just need to find two sets of a, b, c, and d that do the trick. What I did was:

Solve[{a x1+b y1+c z1+d==0, a x2+b y2+c z2+d==0}, {a,b,c,d}]
{{c -> -((a*(x1 - x2))/(z1 - z2)) - (b*(y1 - y2))/(z1 - z2), 
d -> -((a*(x2*z1 - x1*z2))/(z1 - z2)) - (b*(y2*z1 - y1*z2))/(z1 - z2)}}

and the WE told me the relationships you have among them. Then I just arbitrarily chose two sets of a and b (<0,1> and <1,0>) and calculated two sets of c and d. Then I just used those two equations in the Solve command as my (2!) constraints and it worked. And it was really fast.

c[a_,b_]:=0. + 0.3998911137852679*a - 0.2628716901520667*b;
d[a_,b_]:=0. + 0.9517484637723629*a - 0.5145497461933767*b;
r[t_]:={x[t],y[t],z[t]};
c1={0,1,c[0,1]}.r[t]+d[0,1];
c2={1,0,c[1,0]}.r[t]+d[1,0];
KE=1/2 r'[t].r'[t];
PE=4/3 Pi r[t].r[t];
L=KE-PE;
el[a_]:=D[L,a[t]]-D[L,a'[t],t]+lm1[t] D[c1,a[t]]+lm2[t] D[c2,a[t]]==0;
sol=First[NDSolve[{el/@{x,y,z},
                   x[0]==p1[[1]],
                   y[0]==p1[[2]],
                   z[0]==p1[[3]],
                   x'[0]==0,
                   y'[0]==0,
                   z'[0]==0,
                   D[c1,t,t]==0,
                   D[c2,t,t]==0},{x,y,z,lm1,lm2},{t,0,1}]]

What about non-straight line?

A general path in 3D can be given by three parametric equations where you use a fourth variable to connect them all:

• x=f(t)
• y=g(t)
• z=h(t)

That looks like three constraints but when using Lagrange Multipliers you don’t really want to have that extra variable (t in this case). You want to only use x, y, and z and have extra unknown force terms that the ODE solver solves for.

So really you have to invert one of those 3 equations. Often one of them is amazingly simple, like f(t)=t. If that’s true, then you have two constraints:

• y-g(x)
• z-h(x)

and you’re off and running!

Here’s an example where g(x) is sin(x) and h(x) is cos(x). This is a helical path:

c1=y[t]-Sin[x[t]];
c2=z[t]-Cos[x[t]];
r[t_]:={x[t],y[t],z[t]};
KE=1/2 r'[t].r'[t];
PE=-9.8 x[t];
L=KE-PE;
el[a_]:=D[L,a[t]]-D[L,a'[t],t]+lm1[t] D[c1,a[t]]+lm2[t] D[c2,a[t]]==0;
sol=First[NDSolve[{el/@{x,y,z},
                   x[0]==0,
                   y[0]==0,
                   z[0]==1,
                   x'[0]==0,
                   y'[0]==0,
                   z'[0]==0,
                   D[c1,t,t]==0,
                   D[c2,t,t]==0},{x,y,z,lm1,lm2},{t,0,3}]]

Producing this path:

Note that I put in a potential energy in all of these examples. I could have just set that to zero, but then to get things to move I’d have to set the initial velocity. For the straight lines, that’s relatively easy, but for the helix I’d have to think carefully about what direction it starts moving in. So much easier to start things at rest and let a (conservative) force take over.

Your thoughts?

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

• This is super helpful. What I especially like is . . .
• This is dumb. I hope future you realizes that.
• This is future you. Thanks for this, it’s really come in handy.
• This is future you, watch out for that thing tomorrow.
• Really? You’re going to start using WE instead of “the mathematica kernel you can use for free”? That seems too informal.
• I’m still not sure I get why you insisted on adding a potential energy to all these examples.
• Why do you use for the potential energy in a planet?
• What do you mean by a conservative force? What’s wrong with a liberal force?

Does the hoop hop?

I greatly enjoyed this recent video from StandUpMaths:

The set up is a hoop with a mass attached at one point. It’s rough so it rolls without slipping. It’s released with some angular momentum with the extra mass starting at the top. The question is whether the hoop loses contact with the ground as it rolls around. The original analysis from decades ago focused on a mass-less hoop and suggested that the hoop has to hop when the mass has rolled around just 90 degrees.

I wanted to see what it would take to model this. What I thought would be useful would be to use a Lagrange Multiplier technique so that I could easily check the normal force on the ground, watching to see if it ever goes to zero, as that’s when it would hop. If the normal force goes negative, that means the ground is pulling the hoop down. And the ground doesn’t usually do that.

Math/physics setup

When you use Lagrange Multipliers, you tend to have more degrees of freedom in the problem. In fact you have the normal number of degrees of freedom (one in this case: $latex \theta) plus however many constraints you model with Lagrange Multipliers. I want to model both the normal force and the rolling-without-slipping force. For the normal force, my constraint is:

where r is the radius of the hoop and y is the y-component of the center of the hoop.

For the rolling-without-slipping constraint I use:

where x is the x-component of the center of the hoop and is the rotation angle. When the mass is at the top.

The kinetic energy is given by:

Where we treat the center of mass of the hoop as one free particle and the mass as another. Note, however, that we stipulate that the location of the mass is tied to the location of the center:

Then we just do the usual Euler Lagrange tricks and we’re off and running.

m=1;
mhoop=0;
r=1;
g=9.8;
omega=0.1;
xm[t_]:=x[t]+r Sin[th[t]];
ym[t_]:=y[t]+r Cos[th[t]];
KE=1/2 m (xm'[t]^2+ym'[t]^2)+1/2 mhoop (x'[t]^2+y'[t]^2)+1/2 (mhoop r^2) th'[t]^2;
PE=m g ym[t]+mhoop g y[t];
rollingConstraint=x[t]-r th[t];
normalConstraint=y[t]-r;
L=KE-PE;
el[a_]:=D[L,a[t]]-D[L,a'[t],t]+lmRolling[t] D[rollingConstraint,a[t]]+lmNormal[t] D[normalConstraint,a[t]]==0;
sol=First[NDSolve[{el/@{x,y,th},
          x[0]==0,
          y[0]==r,
          th[0]==0,
          x'[0]==r omega,
          y'[0]==0,
          th'[0]==omega,
          D[rollingConstraint,t,t]==0,
          D[normalConstraint,t,t]==0,
          WhenEvent[lmNormal[t]==0,"StopIntegration"]},{x,y,th,lmNormal,lmRolling},{t,0,10}]];
tmax=sol[[1,2,1,1,2]]

The “WhenEvent” trick stops the integration right at the point of when the hoop would hop. I don’t always use that, you’ll see, in my animations below.

Results

Here’s a plot of the normal force for the case of a truly massless hoop:

You see that it starts at 9.8 because the mass is 1kg and the ground is really just holding the hoop+mass up (remember that the hoop is massless). It’s not actually exactly 9.8, though, since right at the beginning there’s a non-zero angular speed, , and so the some of the gravity force needs to be dedicated to centripetal acceleration.

Here’s an animation that stops right at the point of hopping:

Now lets give the hoop just a little bit of mass: 0.01kg (so 100x less than the attached mass). Here’s a normal force plot showing that I allowed the simulation to go past the hopping point. Note that I don’t let it hop. Rather, I allow the normal force to go negative.

Here’s an animation of that setup:

Finally let’s consider a situation where the mass of the hoop equals the mass of the attached point. Here’s a normal force plot:

Note that it never goes to zero. Here’s an animation:

Cool, huh?

Your thoughts?

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

• This is cool! What I liked the most was . . .
• This is dumb! What you totally got wrong was . . .
• Wait, I thought you saw this on a math(s) youtube channel. This seems like physics.
• I don’t think that mass is undergoing circular motion, so saying that some of the gravity is providing a centripetal force doesn’t seem right.
• I suppose you’re not going to bother posting your code, are you?
• I can’t believe you insist on using a super expensive software package like Mathematica . . . wait, what did you say 3 posts ago?
• That youtube video talks about looking at when the hoop starts to slip. Couldn’t you do that by looking at your rolling Lagrange Multiplier and seeing if it gets larger than the expected friction?

Planetary tunnel oscillators

Pick a point on the earth and start digging. It doesn’t have to be straight down. Keep it straight (careful! it’s harder than you might think) and keep going until you come back to the surface. Ok, now drop a ball into your hole. It’ll come back to you eventually because while gravity will pull it down the hole at the beginning, eventually it’ll start getting further from the center of the earth and start slowing down. Assuming the earth’s density is spherically symmetric, it’ll come to a stop right at the edge of the other end, turn around, and come back.

If in addition to being spherically symmetric the earth’s density was also uniform, the travel time of the ball will be independent of the direction you dug the hole. In other words connect any two points on the surface of this idealized earth and the time it takes the ball to tunnel from one point to the other and back will be the same! Pretty cool, right?

I wanted to see if I could 1) make cool animations of that, and 2) see if my intuition would be right for non-uniform spherically symmetric densities. Before you scroll any further, make a prediction: If the planet is more dense towards the center, will a tunnel that gets close to the center take more or less time than a shallow tunnel?

The model

I modeled this using the Wolfram Engine (for free!) which is basically mathematica but using a Jupyter interface. I used a Lagrangian approach which just means I needed to identify the degrees of freedom (the distance down the tunnel) and figure out how to express both the kinetic and potential energy for a particle using that variable. If p1 and p2 are the two points on the surface of the earth, I parametrized my problem with the variable s such that if s is zero we’re at p1 and if s is one we’re at p2. When s is between zero and one it’s in the tunnel:

Then the kinetic energy is:

which will be a function of s.

The potential energy is a little trickier but basically it’ll always be:

where the total Mass, M, that is closer to the center that the point of interest is given by:

assuming r is less than the radius of the planet.

Simulations

Here’s a plot of s as a function of time for a uniform density planet in a universe where G=1, m=1, and the radius of the planet = 1:

Here’s an animation of 10 different random tunnels all dug from the north pole. Note how they stay in sync:

Ok, time to check your intuition. Consider a planet whose density varies as 1/r. That means it’s technically infinitely dense right at the center and it progressively gets less dense as you get further out. Which tunnels will be faster? The ones that are shallow, and so only ever experience the lower densities, which means they always feel nearly the full gravitational pull of the planet? Or the ones that dive deep and explore portions of the planet with high densities?

Here’s what I was thinking: The shallow ones never penetrate where the bulk of the mass is. That means their force stays high. Therefore I think they should be faster.

Here’s a plot of 10 random tunnels, note how they don’t stay in sync:

Here’s the animation:

For me it’s easiest to see the effect if you stare at the upper left purple one and keep your peripheral attention on the central green one. You can tell that the shorter path takes longer! I was wrong!

Ok, what about the opposite: a density that grows linearly with distance from the center. That means there’s no density at the center and the maximum density is at the outer edge.

First a plot of the trajectories:

and here’s the animation:

See how the lower right red one really falls behind? Cool, huh?

Your thoughts?

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

• This is really cool! What about a planet that . . .?
• This is really dumb! What you forgot about was . . .?
• My intuition was right! Here’s what I was thinking . . .
• My intuition was wrong! Therefore all of physics is wrong.
• I just assume the potential energy is always GMm/r and you just have to figure out the M. (that’s what I naively thought too as I started but realized that you have to do the integrations to get it right)
• Did you set all the planets to have the same total mass? (yes)
• Why do you insist on using an expensive piece of software that normal people can’t afford . . . wait, what did you do two posts ago?
• What happens if the density is not spherically symmetric? What are you, lazy?
• You clearly didn’t model the effects of a spinning planet, even though you forced us to click through to a whole other blog post about digging holes. Be consistent!
• Why didn’t you dig identical tunnels for all three planets? It’s super tough for me to directly compare them. My guess is you got all excited about using RandomPoint[Sphere[], n] in the Wolfram Engine.

Dropping ladders

My friend Rhett Allain has really got me interested in this ladder drop posted by Veritasium:

Here’s Rhett’s awesome explanation:

Of course I wanted to see if I could model it with Mathematica, and, after finding I could run Mathematica for free on my windows laptop, I thought I’d revisit this problem, since I think I found something interesting.

Physics involved

Here’s what I did:

• I modeled the rungs of the ladder as rigid sticks with a mass, a length, and moment of inertia
• I modeled the ground as a one-way spring (it only pushes up if either end of the rung goes below the ground, it doesn’t pull back down when either end is above the ground)
• I modeled the strings holding the rungs together as one-way spring as well. If they stretch, they pull. If they’re compressed, they do nothing.
• I modeled frictional energy loss by treating the ground as a viscous fluid. Any time either end is under the ground, there’s a frictional force proportional to the speed of that end that is in the opposite direction of the motion. Here’s a youtube vid I made about that approach.

I did all of that in a Lagrangian formalism. The one-way springs are piecewise potential energy functions, for example.

I had to play around with the strength of the ground one-way spring to make sure the bounces seemed realistic and also with the amount of viscous friction so that the rungs came to a stop relatively quickly.

What it looks like

Here’s an animated gif showing 4 very similar ladders. From left to right:

• Moment of inertia of each rung is
• Moment of inertia of each rung is
• Moment of inertia of each rung is
• Free fall (doesn’t hit the ground so the moment of inertia doesn’t matter)

Here’s a plot (along with a zoomed-in inset) of the height of the middle of the top rung for each of those ladders:

You can see that the smaller the moment of inertia of the rung, the faster it falls.

Using Rhett’s explanation (which I really think is right), when a tilted (that’s important!) rung hits the ground, it slaps down the other end, which stretches the string tied to that end, which then pulls down, at least slightly, on the rung above it. This continues in a domino like fashion. Ultimately all that “pulling down” is what accelerates the top rung faster than free-fall.

The differences for the moment of inertia is what I find interesting. The one is a rung that has a uniform cross-section (ie, the mass is uniformly spread through the length). The one is what you get if all the mass of the rung is at the ends. Really this is the highest possible moment of inertia for a rung. The lowest possible is zero if all the mass is in a point mass at the center, but for fun I did to see a more obvious effect. With a higher moment of inertia, the “slaps down” effect I talk about above is lessened, because the rung is harder to rotate. That means the “pulls down” effect is lessened as well and so you see that there’s less additional acceleration for the top rung. In the plot above you can see that the maximum moment of inertia is indistinguishable from free-fall.

Your thoughts

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

• I love this! I’m going to build a minimal moment of inertia ladder right now.
• This is dumb. What you totally missed in your model was . . .
• I don’t know why you insist on using such an expensive tool like Mathematica . . . wait, what did you say in your last post?
• I don’t like your friction model. Instead I think you should . . .
• I think your rungs are bouncing too high, you should raise your friction amount. (I did that but found that it reduces the “slap down” too as it doesn’t just cause rotation, it also moves the end a little bit. Something to think about, I guess)
• I think all ladder rungs are uniform. This is stupid.
• I think Veritasium doctored that first video. There’s no way there’s two ladders with every rung angle being exactly the same.
• Hang on a sec, it looks like you’ve pinned the center of the rungs on a vertical line. That’s crap and you know it.
• Are you saying that if all the rungs are perfectly horizontal nothing happens? How boring.

Mathematica for free

In this post I’m going to try to capture the steps I took today to get a Jupyter Notebook to run Mathematica commands. I did it on a Windows laptop so if you’re on a Mac of Linux you’ll have to make appropriate adjustments.

Why free?

Mathematica is not free. A single license can be over $1000. Student versions are more like $200 but there’s no question that it’s pricey, especially when compared with python-based computing. Of course, it definitely has value! The amount of development that’s been poured into it, including just in the last few years, is amazing. Want easy to use neural networks? Check. Want lightning fast prototyping? Check. Want access to tons of curated data sets that all play well together? Check.

But, ugh, what a cost. And it’s hard to convince students to embrace it if they know they’ll lose it as soon as they graduate. I’ve spent lots of time looking for a replacement for the types of things I tend to work on, and it’s been a mixed bag.

But ever since I noticed that Wolfram (the company that makes Mathematica) put a free mathematica kernel (what they now call a Wolfram Engine) on the Raspberry Pis, I started keeping an eye out for ways to use the Wolfram Engine for free.

If you go here you’ll see that you can download the engine for developers and use it for free as long as you’re not trying to make money with it. When you download it and install it you can then run wolframscript in a cmd window and you get old-school command-line Mathematica. Graphics don’t work but general calculations do.

To get a more useful interface, you can pair it with a Jupyter Notebook. These are free interfaces that people mostly use to do Python programming. But you can swap out the kernel from Python to Wolfram Engine with the steps that I outline below.

So the kernel (Wolfram Engine) is free (assuming you’re not trying to make money) and Jupyter Notebooks are free. Off we go!

Installation steps on Windows

1. Download and install the Wolfram Engine
   • You’ll have to register on their site (for free) to get it to work.
   • Run wolframscript in a cmd window and it’ll ask for your login credentials for Wolfram.
2. Install python
   • Note, I couldn’t get this to work with Anaconda, so I uninstalled Anaconda and separately installed python, pip, and jupyter
3. Install pip
   • Go to “Method 2” on this page
   • Make sure to install it in the same folder where python is.
     • To figure that out, run py on the cmd line and run this code:

>>> import os
>>> import sys
>>> os.path.dirname(sys.executable)

4. Install Jupyter
   • in that same directory, run this code: pip install jupyter
5. Download the paclet
   • Go to this website
   • Then go down to “Assets” and download the file called WolframLanguageForJupyter-x.y.z.paclet (for me x=0, y=9, z=3)
   • Make sure to save it to the same folder where python and pip are now installed
6. Run a few commands in command-line mathematica
   • In that directory, run the wolfram engine with this command: wolframscript
   • In Mathematica run these commands one at a time:
     • PacletInstall["WolframLanguageForJupyter-x.y.z.paclet"] (make sure to use the right x, y, and z)
     • Needs["WolframLanguageForJupyter`"] (note the extra back tick toward the end there)
     • ConfigureJupyter["Add"]

And you’re done! Now you just need to run jupyter by typing jupyter notebook on the command line and it should launch a new tab in your browser (it takes about a minute on my machine, it seems). Then you can select a new file using the Mathematica kernel:

Finally, here’s a youtube vid of me using it (showing how to model a 10-pendula)

Your thoughts?

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

• This is cool! My favorite part is . . .
• This is dumb. The dumbest part is . . .
• Why don’t you italicize Mathematica?
• Can it do Manipulate commands? (sort of)
• Can it access the curated data? (haven’t tried yet)
• Why didn’t you do an 11-pendulum? Jerk
• I can do all this a lot easier with scipy, loser.
• Can you write this up for Mac users?

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!

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:

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 . 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:

where 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 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:

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 (the x-location of the contact point). Assume .
2. Find expressions for the coordinates and velocities of the center of the ball based on
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
4. Do the usual Euler-Lagrange trick to find the equation of motion for
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 and 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:

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 .

Then step four gives us the following differential equation for :

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 and 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 and as the generalized coordinates. Assume
2. Find expressions for the coordinates and velocities of the center of the ball based on and
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 and .
4. Do the usual Euler-Lagrange trick to find the equations of motion for and
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 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:

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:

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 but what is the magnitude of ? Well, we know what that is in terms of the Euler angles (and their time derivatives)!

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

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 can be found from our simple rolling-without-slipping condition. Putting it all together yields:

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:

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 . 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?

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:

along with the length of the string below the hole:

My mistake was assuming that was equivalent to:

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:

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:

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.