Concerning Contagion

A measles outbreak has been in the news lately, and since the disease was previously considered eradicated in the US, I wanted to try a simple model of disease spread with and without vaccination. To represent a population, I used Python's NetworkX package, which gives us nodes (people) connected by edges (contact that can spread disease). In addition to vaccination status, I wanted to see how different levels of contact affected spread, so I generated the layouts using the Watts-Strogatz Graph, which takes a given number of nodes n, a community size k, and a rewiring parameter pwire. The nodes are arranged in a ring, with each connected to its k nearest neighbors. Then we randomly switch some of those connections according to pwire to randomize the interactions a bit.

Turning now to the nodes of these networks, the model for people I chose tracked a few qualities: Current infection, vaccination/immunity status, susceptibility to infection, and infectious period. Based on some stats I found from the CDC, I decided after being infected a person would be sick for 22 days, but only contagious for days 10-18. While contagious, each day a person's neighbors can become infected according to their susceptibility (0.5% for immune, 50% otherwise). I started the model with 10% of people infected, and ran until all cases were gone or a maximum of 300 steps (days). I ran several different values for k and pwire, as well as different vaccination levels.

After running the simulations, I tried to come up with some summary statistics. The first thing that came to mind was the fraction of people that get sick:

The x-axis shows the k parameter for each run, which represents the number of people each person has contact with. For the lower values of k, this shows clear separation in the different levels of vaccination – Low vaccination leads to high infection rates. With enough connection though, everyone encounters an unvaccinated person and gets sick.

That's just the number of people who ever get sick, but we can also look at how long the disease spreads before being eradicated:

Again we see a big uptick for high community size, to the point where disease is still spreading after 300 steps. For the smallest size, we again see separation in the vaccination rates, but there's also an effect from the pwire parameter. This is easier to understand if we look at a map of one of these communities:

Blue dots are people who have never been ill, red are currently sick, and green are vaccinated, or recovered. The twisting shape shows how the pwire parameter has connected a few people on opposite sides of the initial ring. The black dot in the plot above shows that for low pwire, the disease has to make its way around the whole ring, increasing the time to eradication. We can also look at a case with larger communities and more interconnection:

In this case, almost everyone immediately gets sick, which you might imagine would reduce the time before the disease is eradicated, but the increased connection also results in many more breakthrough infections.

I'm really impressed with the complex behavior this simple model was able to show! That said, I fully recognize it's just a toy, and should not be used for setting policy. I'm not headed to that farm upstate just yet...

Teach Your Physicist to Suck Eggs

Recently I was reminded of a neat science demo where burning a match inside a bottle can suck in a boiled egg that otherwise wouldn't fit:

I saw similar demos growing up, and I believe I was given the incorrect explanation of the match using up the oxygen to lower the pressure. As the video above points out, it's actually the heat from the fire that makes the air expand and escape from the bottle. After the air cools, the pressure is lower and the atmosphere pushes the egg into the bottle. I was curious whether I could determine the pressures involved in this demo.

Burning the match releases heat, which raises the temperature of the air in the bottle. According to the Energy Information Administration, burning a match releases 1 BTU of energy. Using air's heat capacity, density, and a volume of 1 L, we can calculate that this will raise the temperature by an impressive 1220 K. We can use the ideal gas law to see what that does to the air:

Increasing the temperature, T, will cause the pressure, P, to to go up. This pushes up on the egg allowing some air to escape. If we suppose that the pressure required is very small, we can instead assume all the energy from heating the air goes into expanding the volume, V. We could use this to find the final volume of the air, and thus how much escapes from the bottle, but we don't really care about the air that escapes. Instead, we can just suppose our 1 L bottle of air at a room temperature of 300 K is now at its new temperature of 1520 K, with the pressure remaining constant.

Now we can flip the relation, and suppose the volume and number of particles, N, stays constant. We can write

where the subscripts specify the initial and final states. Using the temperature we calculated and 1 atm of pressure, returning to room temperature would bring the pressure to approximately 0.25 atm. That implies the atmosphere is pushing down with 0.75 atm of net pressure. We can get the area of the egg's cross-section and multiply by the pressure for a force of 110 N, or about 25 lbs. Not enormous, but it's easy to see how applying that force unevenly by pushing with your fingers would just mush the egg.

Turning to the latter part of the video, how about getting the egg out again? The maximum pressure a person can blow is about 0.1 atm, 40% of the pressure that pushed it in. If you watch the video above though, you can see the egg goes into the bottle with a great deal of speed, while getting it out isn't nearly as violent, so these results appear to be consistent. Cool experiment to try if you've never seen it before!

Window's Vista

This week, Futility Closet had a post about the Vista Paradox. The paradox is connected to an effect I discussed before, where the angular nature of our vision causes objects to move in unexpected ways. In this case, we're viewing a distant tower through a window. As we approach the window, the window takes up proportionately more of our vision that the tower, causing it to appear to shrink. I tried making an animation of this to test my understanding:

As expected, the tower (red) changes very little, while the closer window (blue) expands rapidly. Initially, I was concerned with how the window's proportions appear to change, but I checked with my mother Sally, and she gets the same expression I did. There still may be an error in my implementation though.

After reading the Futility Closet post, I was reminded of an effect often used in filmmaking called a Dolly Zoom: The camera moves forward (backward) while zooming out (in), causing the edges of the frame to stay where they are, but the center to contract (expand). The scene that always comes to my mind for this is one from Fellowship of the Ring:

makeagif

Using the same framework as above, I plotted a series of rings at different distances, but the same size. Then I moved the camera toward them, while fixing the edges of the plot to the nearest ring:

This effect makes the center appear to shrink, the opposite of the motion in the film, but exactly the same technique. It's neat to see this artifact of human perception appear in both architecture and filmmaking!

Pilfered Power

Some time ago I saw an article about the concept of harvesting power from high-voltage transmission lines. The article reported disturbingly large voltages obtained from the equipment they set up, which led them to stop the experiment for the sake of safety (and legal protection). I was curious if I could figure out a way to predict the voltage we could expect from such a situation.

The electric system transmits power with alternating current (AC) at high voltages, which gets stepped-down by transformers when it is delivered to our houses. To see why the transmission is done at high voltage, we can imagine sending out a power P through lines with resistance R:

The power lost is inversely proportional to the square of the voltage, meaning raising the voltage can quickly decrease the power we lose to resistance.

When current travels through a wire, it produces a magnetic field in circles around the wire. The equation for this field is

where μ is the permeability of free space, I is the current in the wire, and ρ is the radial distance from the wire. For alternating current we can write

where I0 is the amplitude of the current and f is the oscillation frequency, 60 Hz in the US. If we set up our fence directly below the transmission line with the z-axis pointing along the fence and the y-axis pointing up toward the wire, then the B-field at this point will be

where h is the height of the transmission wire, and we've replaced I0 with the power & voltage carried by the line. Now to get the voltage induced in the fence, we can use Ampere's Law:

Plugging in our B and integrating to get E gives

Electric field is measured in volts/meter, so this will tell us the relation between the length of our fence and the voltage we can expect between the two ends. The example case discussed in this Wikipedia article uses P = 1000 MW and V = 765 kV. This page gives the minimum height for wires carrying that voltage as 45 ft. Plugging in those values, along with the 60 Hz frequency I mentioned, gives 330 megavolts/meter, which is pretty insane! Now I did assume best (or maybe worst) case numbers in this calculation, but if you watch the video in the article that started this, the presenter seems a bit concerned with the size of the spark he gets at the end. Good (grounded) fences make good neighbors!

Duck and Cover

[I considered "Don't Look Up," but I have a soft spot for Red Scare propaganda.]

This week there was news of an asteroid with a 1-in-82 chance of hitting the Earth in 2032. Most reports seemed to spin this as a reassuringly low chance, but as someone well-versed in low probabilities, I found it to be uncomfortably high, so I wanted to dig into it a bit. NASA's Small Body Database has an entry for this asteroid, which gives its orbital parameters and a nice interactive plot of its path through space. In order to do some calculations with this info though, I found a tool developed by MIT called poliastro, which can connect to the database and read the orbit definitions.

As a first pass, I tried duplicating the type of plot NASA included on its page, showing the orbit of the asteroid compared with a few of the planets around the time of the possible collision:

This shows the orbits of Earth in blue, Venus in cyan, Mars in red, and the asteroid in magenta. The planets' orbits are all fairly circular, while the asteroid's is sharply elliptical, and crosses the paths of the planets. The predicted flyby is around December 2032, and if you watch the animation you can see the asteroid skirts behind the Earth with a fair gap.

All predictions include errors though, so what happens if this orbit isn't quite the right one? Poliastro offers an object called a maneuver that can perturb a given orbit. Starting with the orbit above, I applied a random change in velocity with magnitude up to 10% of the initial value. Then I measured the distance between the Earth and the asteroid in time for each variation in the orbit:

The x-axis is the Modified Julian Date, a measure of time often used in astronomy. The black line is the predicted orbit, which only gets down to around 200 million km separation. For reference, that's about 1,000 times the distance from Earth to the Moon. However, some of those variations get a lot closer. To see how many of them are getting close, we can look at a histogram comparing the closest approach for each of the variants to the original:

Values below zero indicate the variant gets closer than the base prediction, and the red line marks the distance that would result in a collision with the Earth. Several of those variants are getting awfully close, so much like my pressure cooker post, I'm not sure I feel better at the end of this!

Water Wedge

The streets around our house have many potholes, which fill up with ice:

Water expands when it freezes, which means when a hole fills with water that freezes, the ice can potentially make the hole bigger. To do this, the ice puts pressure on the asphalt, and vice versa. Adding pressure can make it more difficult for ice to freeze, so I looked up some info on asphalt strength. This paper gave a plot comparing the bond strength for different aging times:

Figure 8

These are all around 2-4 MPa. We can look at the phase diagram for water to see how this would affect freezing:

Wikipedia

This is well within the range where the freezing point is still around 0°C, so we don't need to worry about the pressure we're applying. Water expands by about 9% when it freezes, so we can imagine water filling a crack, then pushing out in all directions:

That's the general idea I had in mind, but getting into the details, the steps I used are

1. Begin with a flat road. Put a dent in the center with depth y0 (*parameter).
2. Find all points below the road, and use them to define a polygon.
3. Find the centroid (center of mass) of the polygon and push all points out by 9%.
4. Push the polygon down according to its weight (*parameter).
5. Add points anywhere the separation has passed a threshold.
6. Adjust points along line for equal spacing.
7. Repeat from 2.

I added steps 5 and 6 since otherwise the crack just scaled up from the initial dent. We can look at the type of crack these steps give as we cycle through freezes/melts:

This seems a little extreme, but we can take a look at how this changes with the two parameters I marked above. The two qualities I was mainly interested in were the width and depth of the crack. We can see how those change in time with different y0 and weight values:

It makes sense these curves look exponential, since the increase depends on the current size. What really surprised me though was how well-separated these two qualities are according to the input weight and y0. These plots suggest that the width is almost wholly determined by the weight, and the depth by the initial depth. For the depth, this can easily be tied to the exponential behavior – Starting deeper increases the rate the depth increases. A little less obvious is the weight, but if we look at the initial diagram of the setup, the wedge shape of the ice means pushing down also pushes out. A quality of exponential systems is that they can very quickly get away from their starting point if left unchecked – A reminder to keep those roads smooth!

Beet the Traffic

We watch the local news every morning, and this week a story caught my interest: Alternatives to using road salt to avoid harmful risks. The key point of the story is that runoff from salting the roads can damage the surrounding environment, but adding beet juice to the solution can make it stick to the roads better. Reading up on the idea, I've found it's still debated whether this idea is really better for the environment, but the bit I was curious about is the ability to stay on the road.

Looking into ways I could model this, I found a paper discussing how droplets spread over time based on their surface tension. Their model was a bit more involved than I wanted to get, so I made some simplifications: The droplet takes the form of an ellipsoid with constant volume and circular base. This constrains the relationship between the height and the radius. The paper defines a value they call h* where gravity and surface tension balance. After my simplifications, it takes the form

where ρ is the density of the fluid, g is the acceleration due to gravity, and γ is the surface tension. For a given fluid, we can look up the surface tension and density. I decided to try a salt-water solution (γ, ρ), a sugar-water solution (γ, ρ), and molasses (γ, ρ).

The paper gives a t^1/5 form for the height of the drop, so we can start each of these fluids as a hemisphere and see how they spread as they approach their respective h* values:

This shows the sugar and salt spreading at roughly the same rate, contrary to the idea given in the report, so I expect my model is not capturing all their qualities. What I find really interesting though is that the thicker molasses actually spreads faster, because it's more dense that the other two, so gravity exerts a stronger force. Naturally, this brought to mind a bit of history from my home state, when a flood of molasses from a ruptured tank cut a swath of destruction through Boston!

Make Like a Tree and Get Out of Here

Swarthmore's Scott Arboretum traditionally gives incoming students a plant to care for in their dorms. Shockingly, the Hawaiian Schefflera I received in 2007 is still going 17 years later! Over the years as I've moved from place to place and kept the plant in different environments, I've been impressed by its ability to track the sunlight, frequently growing lopsided as it reaches toward the nearest window until I think to turn it. This has resulted in some twisting, gnarled branches:

If you've been reading this blog, you can probably guess where this is going – I was curious if I could make a simulation of my plant's heliotropic tendencies. I decided to model the plant as a collection of connected branch segments, each a fixed length and pointing at an angle relative to the vertical. At each step, we iterate over all the segments and pick an action:

• If the segment has no children, i.e. it's at the tip of a branch, we add a new segment on the end with a probability p_grow/size, where size is the number of existing segments.
• If the segment does have children, we add a new one with probability p_sprout/size.
• If neither of those occur, we adjust the angle of the branch to point closer to the sun's current position. The adjustment is proportional to how far off the angle is, how many branches are on the end of this one, and a constant stiffness for the plant.

I tried a bunch of values for the different parameters until I landed on a range that gave plants looking reasonably similar to the real thing (click to enlarge):

The numbers along the top give the stiffnesses, and the ones on the left give the sprout probability. The sun moves back and forth sinusoidally, which you can see in the snaking of the plants. I wasn't able to get my digital plants to spread as much as the analog one, possibly because I'm not accounting for the plant casting shadow on itself, but I'm still pleased with the results – The top center one seems particularly good. If you'd like to try for yourself, the code is here.