Modeling the Spread of COVID-19: Part 3 - Mitigating Strategies

In the previous columns on this topic, I have looked at the problematic criterion for lifting the lockdown in Pennsylvania, and the statistics behind whether you get infected if exposed.  In this column, I will discuss the efficacy of mitigating strategies like mask wearing and social distancing.

The most visible of the public health strategies thus far has been the lockdown.  We have all been hunkered down in our houses, afraid of our groceries, and suddenly, we're all home baking bread because gluten intolerance is no longer a thing (/s).  Additionally, the skies are clearing and due to a lack of traffic, there is verifiably less pollution all around the world.  All of this is in the name of social distancing.  If you aren't exposed to the virus, you can't get it.  Lockdowns work.  Period.  While there is some debate about how long they should last to be effective, there is evidence that a five week lockdown could end the pandemic once and for all (if you can't get it and it runs its course, you can't spread it!).  Others feel that because such lockdowns can never be total (where will you get food? who keeps the lights on? etc.) it will be more like 18 months or about when a viable vaccine is developed.

Social distancing measures can also be used in a less all-or-nothing manner.  The Centers for Disease Control (CDC) recommends, among other measures, when in public spaces keeping a 2-meter distance from others who aren't in your household in order to prevent the spread of COVID-19.  However, if someone who is infected is not wearing a mask, their region of influence, that is, how far they can spread the virus, is not limited to two meters.  Jayaweera et al just published a paper titled: "Transmission of COVID-19 virus by droplets and aerosols: A critical review on the unresolved dichotomy" which calls into question the 2-meter doctrine, especially in closed spaces.

Other supporting work includes that of Qureshi et al (2020) in a paper titled: "What is the evidence to support the 2-metre social distancing rule to reduce COVID-19 transmission?"  Their results can be summarized via their graphic here:

The paper is a fantastic review of prior work by other authors and was conducted by members of the Center for Evidence Based Medicine at Oxford University.    Another work included in the summary is published in the Lancet by Chu et al (2020).  This work looks at the odds ratios (risks) of contracting COVID-19 and considers a number of factors.  We'll be using some of the numbers from this review later in this column.  In the graph below, we can see how social distancing can decrease the risk of infection.

As you might imagine from the Jayaweera et al graphic, mask wearing is one of the other major mitigation strategies.  Masks work mainly by keeping the infection in (you cough into your mask, and the virus is trapped) rather than by keeping the virus out.  How well do they work?  Well that depends on whether you are social distancing, what type of mask you have, and how well it fits.  This is, again, covered in the Chu et al paper in the Lancet, but the bottom line is that N95 masks work much better than surgical masks, but any mask is better than no mask.  Thicker masks that fit well protect you better than thin masks with gaps.  Most masks are relatively effective in stopping YOU from spreading your germs.

OK, so now that we have some background, let's try to model what happens with or without social distancing, and with or without mask wearing.  Let's do some probability calculations to get a sense of how our prior probabilities are modified by either wearing a mask -OR- social distancing (or both!).  In the last column, the probability of contracting COVID-19 for an 18-44 year old was modeled based on the number of contacts one had.  Recall that for this group, the probability of infection at close quarters was 0.022 or 2.2%  On the x-axis is the number of infected people one is exposed to, while on the y-axis is the probability (0.50 = 50%). 

This result was a little bit scary in that it implies that without any precautions, one can very quickly get infected.  It does not model the infectiousness of a given carrier, the time you have spent in their proximity or anything of that nature.  It assumes that all carriers are equal and you've been exposed at close quarters.

The Lancet paper offers some estimates of how masks and social distancing affect infection using what is called an "odds-ratio".  The odds-ratio is the probability that an event happens divided by the probability that an event doesn't happen.  In the case of rolling a 6-sided die and trying for a one, the odds ratio would be: (1/6)/(5/6) = 1/5 = 0.20.  So while your probability of rolling a one is 1/6 = 0.166666..., your odds of rolling a one are 1 in 5.  In the Lancet paper, the odds-ratio for how mask wearing reduces transmission (compared to not wearing a mask at all) is 0.15 with error bounds of 0.07 to 0.34.  The odds-ratio of social distancing at a distance of one meter or more is 0.18 with error bounds of 0.09 to 0.38.   To be clear, there really are a continuum of values here dependent on a number of factors, but starting with the mean values is a way to get started.

Let's consider the social distancing scenario.  If we use the mean odds-ratio of 0.18, we could calculate the probability, but that would be extra work in this case.  We already know that compared to NOT social distancing the ratio probabilities of getting COVID-19 with social distancing to getting COVID-19 without social distancing is 0.18/1.0.    This means that we can just modify any probability we calculate by multiplying it by the odds-ratio. 

• Probability of getting COVID-19 if exposed to an infected individual = 1.0*0.022 = 0.022 or 2.2%.
• Probability of getting COVID-19 if exposed to an infected individual, BUT with at least 1.0 meter of social distancing: 0.18*0.022 = 0.00396 = 0.396%

Likewise, we can determine how wearing a mask affects things. 

• Probability of getting COVID-19 if exposed to an infected individual, BUT with a decent mask on: 0.15*0.022 = 0.0033 = 0.33%

And, we can compute how the two effects work together, presuming they are truly independent:

• Probability of getting COVID-19 if exposed to an infected individual, BUT with a decent mask on AND with social distancing: 0.15*0.18*0.022 = 0.000059 = 0.0059%

OK, this is looking pretty good.  So instead of having a 2.2% chance of getting COVID-19 with every contact with a sick person (this is a bit more than a 1 in 50 chance), we have reduced our risk to something like 59 chances out of a million or a bit less than 3 in 50,000. 

Now, let's put all this in context with the number of contacts we have who have COVID-19.  Here are the initial parameters.

p = 0.022    // probability of getting COVID-19 from close contact
pNot = 1 - p    // probability of NOT get it
pSD = 0.18;    // Odds ratio of social distancing converted to probability
pM = 0.15;    // Odds ratio of mask wearing converted to probability

Then when we calculate the probabilities for the number of contacts (i+1) the code will be as follows:

  pInf[i] = (1 - pNot**(i+1));
  pMask[i] = pM*pInf[i];
  pDist[i] = pSD*pInf[i];
  pMD[i] = pM*pSD*pInf[i];

Recall that the a**b notation in JavaScript means a to the b power.  Other languages like using the "^" character. 

When we put all this together and rerun our simulation with average values, our results look like this:

Look at that!  By using these mitigating measures, especially in combination, we can really substantially reduce our risk.  And this is when we are pretending that everyone we meet is infected.  This sounds great!   So let's go with the social distancing bit and all wear masks, and we're probably pretty good!  Again not so fast!  While these measures make a difference, they don't make you bulletproof.

In the next column, we'll look at the incidence of COVID-19 in the population, and how this affects the model.  Hint - currently the incidence is quite low, so all these numbers will go down.  But at that point, we'll be ready to start considering the ensemble risk due to a group of people, and how risk increases from day to day.

In the meantime, the take home message for today is wear a mask and social distance!  Especially since the number of new cases in Pennsylvania for June 22 is 456.  The epidemic is far from being under control.

References:

M. Jayaweera, H. Perera, B. Gunawardana, J. Manatunge, Transmission of COVID-19 virus by droplets and aerosols: A critical review on the unresolved dichotomy, Environ Res. 188 (2020) 109819. https://doi.org/10.1016/j.envres.2020.109819.

Zeshan Qureshi, Nicholas Jones, Robert Temple, Jessica PJ Larwood, Trisha Greenhalgh, Lydia Bourouiba. What is the evidence to support the 2-metre social distancing rule to reduce COVID-19 transmission? - CEBM, (June 22, 2020) https://www.cebm.net/covid-19/what-is-the-evidence-to-support-the-2-metre-social-distancing-rule-to-reduce-covid-19-transmission/ .

D.K. Chu, E.A. Akl, S. Duda, K. Solo, S. Yaacoub, H.J. Schünemann, D.K. Chu, E.A. Akl, A. El-harakeh, A. Bognanni, T. Lotfi, M. Loeb, A. Hajizadeh, A. Bak, A. Izcovich, C.A. Cuello-Garcia, C. Chen, D.J. Harris, E. Borowiack, F. Chamseddine, F. Schünemann, G.P. Morgano, G.E.U.M. Schünemann, G. Chen, H. Zhao, I. Neumann, J. Chan, J. Khabsa, L. Hneiny, L. Harrison, M. Smith, N. Rizk, P.G. Rossi, P. AbiHanna, R. El-khoury, R. Stalteri, T. Baldeh, T. Piggott, Y. Zhang, Z. Saad, A. Khamis, M. Reinap, S. Duda, K. Solo, S. Yaacoub, H.J. Schünemann, Physical distancing, face masks, and eye protection to prevent person-to-person transmission of SARS-CoV-2 and COVID-19: a systematic review and meta-analysis, The Lancet. 0 (2020). https://doi.org/10.1016/S0140-6736(20)31142-9.

Modeling the Spread of COVID-19: Part 2 - What are my chances?

Many universities and schools are planning to go back into session in person in the fall.  In a later post, we'll look at why human nature says this is a Very Bad Idea, but for now, let's just get familiar with the statistics involved with understanding the spread of infection.  In this post we're not going to be too complicated in our modeling.  For example, we aren't going to get into R values, incubation times, or such, we're going to focus on pure probabilities and then see how that plays out with our chances of getting infected.

OK, if you've read this far, you've probably got a passing interest either in COVID-19 or probability.  Before we delve into the world of disease, let's first think about one of the common examples used to teach and illustrate probability, the roll of a 6-sided die.  Since it has 6-sides, if it is a fair die, then your probability of rolling any number, say 1, is 1 chance out of 6, or 1/6 = 0.166666... To roll a one twice in a row is, we can agree, more difficult.  Its probability should be lower.  One can find the probability making a table of possible rolls and counting up the possibilities that present themselves.  Below [1,1] stands for a roll of 1 followed by a second roll of 1.

• [1,1], [1,2], [1,3], [1,4], [1,5], [1,6]
• [2,1], [2,2], [2,3], [2,4], [2,5], [2,6]
• [3,1], [3,2], [3,3], [3,4], [3,5], [3,6]
• [4,1], [4,2], [4,3], [4,4], [4,5], [4,6]
• [5,1], [5,2], [5,3], [5,4], [5,5], [5,6]
• [6,1], [6,2], [6,3], [6,4], [6,5], [6,6]
  

As you can see, the rolls [1,1] only shows up once out of 36 times, so it can be seen that the probabilities per roll (1/6) just multiply: (1/6)*(1/6) = (1/36).  

So it is much more probably NOT to roll two ones in a row.  How much?  This is 1 - (1/36) = (35/36).  In the same way, the probability of NOT rolling a one the first time is 1 - (1/6) = (5/6).  On the other hand, if we change the question and ask, "what is the probability of rolling a one at least one time?" the odds change.  Per roll, the probability is (1/6), but the ensemble probability is now greater because I would have more tries.   Let's look at that table again, but color in any combination roll which has a one.

• [1,1], [1,2], [1,3], [1,4], [1,5], [1,6]
• [2,1], [2,2], [2,3], [2,4], [2,5], [2,6]
• [3,1], [3,2], [3,3], [3,4], [3,5], [3,6]
• [4,1], [4,2], [4,3], [4,4], [4,5], [4,6]
• [5,1], [5,2], [5,3], [5,4], [5,5], [5,6]
• [6,1], [6,2], [6,3], [6,4], [6,5], [6,6]
  

In this case, I have a probability of (11/36) of rolling a one at least once in two rolls. To come up with a closed form formula to represent this, we have to think a bit backwards. First, "what is the probability of not rolling a one in one roll?". From before we found that this was 1-(1/6) = (5/6). Next, "what is the probability of not rolling a one in two rolls?" From our previous paragraph we expect that that is (5/6)*(5/6).   Now we can ask, "what is the probability of NOT 'not rolling a one in two rolls'?"  This must be 1 - (5/6)*(5/6) = 36/36 - 25/36 = 11/36.  So we have a way to represent what are called independent probabilities.  The formula for such a prediction reduces (in JavaScript) to:

pOne = 1/6;                // Probability of rolling a one

pNotOne = 1 - pOne;    // Probability of NOT rolling a one

nRolls = 2;                  // Number of rolls

pAtLeastOne = (1 - pNotOne**nRolls)

where the ** operator in JavaScript means "raise to a power".   What you see here in the last formula is an expression that can be reused later.

How does this relate to the spread of disease?  When doing a simple model, one might say that you have X% chance of catching it.  A percentage can be represented in terms of a probability.  Fifty percent is 50/100 or 1/2 or 0.50.  Two percent is 0.02.  You get the idea.  If you search the literature, you can get a sense of just how probable it is to catch COVID-19 from an encounter with someone who has it.  Now let's be clear about this.  Such a value is an average from the analysis of a great many encounters of various durations and conditions which involve everyone from the mildly contagious to superspreaders.  So such probabilities are very rough at best and are good for making general inferences. 

One such paper, Modes of contact and risk of transmission in COVID-19 among close contacts is found on medrxiv.org, a repository for preprints of medically oriented scientific papers.  This group of researchers examined 4950 close contacts in Guangzhou, China and estimated risk of infection based on age group among other things.  Close contact settings were broken into groups including "Cruises", "Public Transport", "Healthcare Settings", "Households", and "Multiple Modes".  Multiple modes".  Of these, there were more cases due to multiple modes of contact (12 out of 92 individuals in the category) than in any other group.  This is approximately a 13% rate of infection.  By age group, the breakdown was as follows:

• 0-17 years old: 14/783 = 1.8% infection rate
• 18-44 years old: 51/2338 = 2.2% infection rate
• 45-59 years old: 29/997 =  2.9% infection rate
• 60+ years old: 35/824 = 4.2% infection rate

Given these data, we can use the same technique to estimate risk of infection from coming into contact with someone who is infected.  Let's say that you are between 18-44 years of age and have a 2.2% chance of being infected if you are in close contact with an individual with COVID-19.  If you just meet one such person, it's pretty straight forward.  But what if you meet 2 or more?  While your chance of being infected each time will, on average, remain the same, your overall chance of being infected should increase.  You could phrase it as this: "What are my chances of being infected at least once?"  We can use our formula from above.  Instead of writing it like this:

pAtLeastOne = (1 - pNotOne**nRolls)

 

we can write:

pInfectedOnce = (1 - pNoInfection**nContacts)

In this case, pNoInfection = 1 - pInfection = 1 - 0.022 = 0.978 (or a 97.8% chance of not being infected).   Let's run some numbers to see if this behaves the way we think it should - more contacts equals a greater chance of infection.

• pInfectedOnce(1) = (1 - 0.978**1) = 0.022
• pInfectedOnce(2) = (1 - 0.978**2) = 0.0435
• pInfectedOnce(3) = (1 - 0.978**3) = 0.0645
• pInfectedOnce(4) = (1 - 0.978**4) = 0.0851
• pInfectedOnce(5) = (1 - 0.978**5) = 0.105 ...

Remember that a probability of 0.105 is a 10.5% chance of being infected.  This is what this function looks like.  By the time we get in close contact with about 35 people we have a 50% chance of contracting COVID-19. [The code to produce this is included below].

All this presumes that 1) the paper's estimation of contagion is correct, 2) we're all in relatively close contact, 3) no-one is taking any precautions.  Under what real world conditions are these likely to be true?  My guess: bars and restaurants, frat parties, and the like.  Where have we seen such uncontrolled spread of the virus?  Where large numbers of people congregate.  In Missouri, for example, there were raucous, uncontrolled parties over Memorial Day weekend.  Look at the increase in case-number about June 14, about 3 weeks after Memorial Day.  I suspect that this trend will continue with time.

And in PA, there's this:

"PHILADELPHIA (CBS) — Twelve people in Bucks County who attended Memorial Day parties at the Jersey Shore have tested positive for the coronavirus. The Bucks County Health Department discovered this cluster of COVID-19 cases through contact tracing. One positive case led to the 11 others. (emphasis added)"

So it's clear that throwing caution to the wind is a real problem.

What about classes?  Well in theory, classes are supposed to start social distancing and wearing masks - which could mitigate the potential for contagion.  In theory, universities are going to try to apply sanctions to students (or faculty, one would hope) who do not take precautions.  In practice, we'll see.  Behavior in a class is one thing, but since most time is spent outside of class, it's not clear what steps universities are willing to take to police this without alienating their paying population.  In the next installment, we'll look at how taking such precautions reduces the likelihood of transmission based on the evidence in the wild.

Code for use in: http://www.niiler.com/JSvg2/analysis.html to produce graph of probability of infection versus number of contacts.

p1 = 0.022    // probability of getting COVID-19 from close contact
p2 = 1 - p1    // probability of NOT get it

N = 50;        // Number of contacts

Ncontacts = [];    // Array to store numbers of contacts (for the x-axis)
pInf = [];        // Array to store the resultant probabilities (for the y-axis)

// Loop through calculating the probability with each number of contacts (i+1)
for (var i=0; i< N; i++) {
  pInf[i] = (1 - p2**(i+1));
  Ncontacts[i] = i+1;
}

// Labels for the graph
optsSeries[0].name = "Model"
optsGraphTitles[0].name = "Number of Contacts";
optsGraphTitles[1].name = "Chance of being infected";
optsGraphTitles[2].name = "";
optsGraphTitles[3].name = "";

// Graph it
drawSeriesPlot([Ncontacts], [pInf], 0,0,'')

Modeling the Spread of COVID-19: Part 1 - Why a lockdown is still needed

As the first wave of COVID-19 in the United States starts to recede, at least for those of us in the Northeast, there is talk about reopening and going back to school in the fall (or even before that), because, well, it's been long enough, dammit!   Current plans are that Pennsylvania will reopen on June 26th according to the governor's office.  Their guidelines for "going green" are that:

"A target goal for reopening was initially set at having fewer than 50 new confirmed cases per 100,000 population reported to DOH in the previous 14 days."

This sounds pretty reasonable, right? 

What does this mean?  There are 12.7 million people in the state of Pennsylvania.  So if you take the number of new cases, divide it by this population, and then multiply by 100,000, you will get the number of cases per 100,000 population.  As of June 21st according to Google, there were 464 new cases in PA.  So 464 cases out of 12.7 million (the population of PA) is 3.6 per 100,000 individuals.  This means we're safe, right? 

Not so fast! The guidelines don't call for fewer than 50 new confirmed cases per 100,000 per day, they are stricter than that.  They call for fewer than 50 new confirmed cases per 100,000 per two week period.  Using the data from the state's own website, we can quickly get this total for the last 14 days - and it is 6293 new cases in two weeks.  This divides down to 49.55 cases per 100,000 as of today (June 22).   Therefore we must surely be safe!

I want to point out that in the context of the epidemic thus far, this number (while impressive sounding) is ultimately pretty arbitrary.  While it is true that there are more options to treat and there is more contact tracing, the number of new cases yesterday, 464, is 23 times higher than the number of cases per day on March 16, one day after the lockdown in PA began.

On March 15 in Chester County, Pennsylvania, the schools closed, and a lockdown was initiated. On March 16, according  to the county data dashboard, there were 16 positive COVID-19 tests.  The reason the curve shown above is falling is due to mitigating measures.  Going "green" at a point where we have many more cases when things started removes many such measures and ensures that the curve will rebound.  This is especially true in light of the fact that PA does not currently have sufficient contact tracing.

Unfortunately, the virus doesn't particularly care about our preferences, politics, or anything at all, really.  It's an inanimate bit of genetic material in a protein and lipid capsule which is a little nano-machine.  It doesn't set out to give you pneumonia or to clot your blood, these are side effects of its replication program, and evolution has determined that by behaving in this way, it has the maximum chance of replicating.  (Note: it's not even trying to replicate, replication is just all it can do.)  Furthermore, depending on the R value (which is the number of individuals infected per person infected), the growth rate of this can well go exponential.  Starting with 464 new infections per day is much worse than 16 as in March so far as exponential growth is concerned.  So unless I miss my guess, we're going to be back to square one or worse in 2-3 weeks, and the Northeast, and PA in particular are going to have a second wave.

Teaching Gauss' Law with Equipotentials

How do you teach Gauss' Law?

It's a serious question which depends largely on the audience. Students who are learning this in an algebra-based physics class are often learning Gauss' Law as a "gee-whiz" phenomena rather than as a serious computational tool, and as such, they may ultimately memorize a handful of charge distributions for which it works (if that!). Students in a calculus-based class will often start with flux and then make it to the usual Maxwell equation:

Next they learn to apply this equation to a bunch of geometric situations. However, for many students, this smacks of black magick and arcane knowledge. How, for example, do we know the exact direction of the electric field? How do we choose a Gaussian surface? Although one can discuss surface normals and trying to ensure that the Gaussian surface normals are, in fact, parallel to the electric field, this sometimes seems circular to students. The result is memorization of certain charge configurations.

As I've thought about this over the years, it has occurred to me that starting from equipotential surfaces might be the way to go. Consider a 2D ring of charge. The following image shows what this charge's scalar potential field will look like. The charge has, of course been broken into a number of discreet charges, so it doesn't look like a perfect ring of charge. But it's close, really close.

The units of this graph are in volts. Equipotential curves are drawn in white around the distribution. The equipotential curves for such symmetrical distributions indicate where the Gaussian surface should go. If your students are used to visualizing equipotentials, then they also have a sense of which Gaussian surface to pick. And so long as the charge distribution is fairly symmetrical, the electric field on that surface will be constant. (I won't get into trying to use Gauss' Law with dipoles or the like, since then the electric field has azimuthal components as well, and Coulomb's Law is probably easier to apply.)

Code to generate this sort of charge distribution (in the interest of visualizing equipotentials) using JSvg2 follows below. I'll discuss the generation of equipotential graphs in another post.

First the code to calculate the Coulomb potential at a field point given the position and charge of the source charge, and the position of the field point.

function calcPotential(qS, xS, yS, xF, yF) {
  // qS - source charge
  // xS, yS - source x,y positions
  // xF, yF - field x,y positions
  
  var k=8.99e9; // N*m^2/C^2
  var thresh = 20; // Max abs value of voltage to display
  
  var dx = xF - xS; // Difference in x and y between source and field pts
  var dy = yF - yS;
  var r = Math.sqrt( dx*dx + dy*dy );  // Pythagorean distance 
                                       // between source and field points
  
  var V = k*qS/r; // Potential in units of volts
  
  if (V > thresh) V = thresh; // For display purposes
  if (V < -thresh) V = -thresh;
  
  return V;
}

Next, define some preliminaries.

In this diagram, the boundaries of the scalar potential field are given by (xo,yo) on the upper left, and (x1,y1) on the lower right. The coordinate system for computer graphics is left handed, unfortunately. We define a total charge, Q, which will be spread out into N discreet point charges along a circle of radius R1. The angular increment between point charges is dphi. We will store our scalar potential field in the array V, and our charge distribution from which it is generated in the array chargeArr.

// Here V is a 1D array
var V = []; // A global array in which to store my potential field
var d2r = Math.PI/180; // Degree to radian conversion

// Define the boundaries of our potential field
var xo = -10.5; // m
var yo = -10.5; // m
var x1 = +10.5; // m
var y1 = +10.5; // m

var Q = 1.0e-8; // total charge on the bar
var N = 90; // There will be N charges!
var dphi = 360/N;   // Angular increment to go around the circle

var R1 = 5.0; // m  - Radius of circle
var dq = Q/N; // the amount of charge - on each charge in circle
    
var chargeArr = []; // new array to store charges

Next, we fill the array by starting on the x-axis at a distance R1 from the origin, and then proceeding counter-clockwise around the circle. As we change the phi position, we recalculate the x and y positions of the point charges using trigonometry. Each charge in chargeArr is itself an array of three elements, [q,x,y], that is, charge, x-position, and y-position.

var ya = 0;   // Initial y and x position of charges in array
var xa = R1;

for (var i=0; i< N; i++) {
  
  chargeArr[i] = [dq,xa,ya];
  
  // Update angle around the circle
  var phi = (i+1)*dphi;
  var phiR = phi*d2r; // Then convert it to radians

  xa = R1*Math.cos(phiR); // New x and y positions
  ya = R1*Math.sin(phiR);
  

}

Next, we actually calculate the potential at each point in the grid. We start at xo, yo. Then we step through all the y-values of the grid (changed by adding dy1 to the prior y-value). Then, we set y back to yo, add dx1 to xo, and then step through y again.

// initialize 2 counters - one for the x direction: i
//                       - one for the y direction: j
// The purpose of the counters is to index the V array.

var i=0; // Counter along x direction
var dx1 = 0.5; // scalar field sampling is 1/10 of a meter
var dy1 = 0.5; //  instead of a meter

// Use embedded loops to change field positions
for (var x=xo; x<=x1; x+=dx1) {
  V[i] = []; // Initialize ith value of V array as an embedded 1D array
  j=0;  // Counter along y direction, initialize to zero prior to
        // looping over y.
  for (var y=yo; y<=y1; y+=dy1) {
    
    V[i][j] = 0; // Initialization
    //console.log("V[i][j] = "+V[i][j]+", x = "+x+", y = "+y) 
    
    for (var ii=0; ii< chargeArr.length; ii++) {
      var charge = chargeArr[ii]; // The ii-th charge in chargeArr
      var qa = charge[0]; // the charge of this charge
      var xa = charge[1]; // the x and y positions of this charge
      var ya = charge[2];
      
      // Update the potential using superposition
      V[i][j] += calcPotential(qa, xa, ya, x, y);   
    }
    
    j++; // Increment this each time through the y loop 
  } 
  i++; // Increment this each time through the x loop
}

Finally, we use some plotting code from JSvg2. Note that everything above this is doable in generic JavaScript. What comes below is not. You need to have the JSvg2 library and its associated functions to work. In any event, the inputs are the colortable of choice, the x and y dimensions of the image, and the number of contours you want drawn. This code takes the input array V and performs a linear interpolation across x and y to estimate in-between values before displaying them.

colortable = 0; // What color do you want the output to be
xdim = 500; // How big a grid do you want
ydim = 500; // 
nContours = 15; // How many equipotentials should be drawn

// Draw it - only works on http://www.niiler.com/JSvg2/analysis.html
field_setupColorMap(V,colortable,xdim,ydim,nContours);

The class had previously worked at creating all kinds of equipotential fields from different charge distributions, so when I projected this image on the whiteboard, they knew what they were looking at. I then had a student come up and draw the electric field lines that corresponded to the equipotentials. He nailed it. They were all perpendicular as should be. (His drawing was somewhat similar to what is below - arrows pointing outwards)

Next, I asked the class how many field lines went through each circle. Answer: the same number. Then I asked if the field was getting stronger or weaker as one moved away. There were several answers. Some students remembered the Coulomb field of a point charge decreasing as one over r-squared. Others said that the field lines were spreading out as one moved further from the charges, and that therefore, the field strength was decreasing. Then I told them to imagine this as a 3D scenario. What shape would the equipotentials be? "Spherical" was the answer. Which way would the field lines be pointing compared to the equipotentials? "Right angles." Next: how did the field strength times the surface area compare from an inner surface to an outer surface. After some discussion, the students arrived at the answer: "it didn't." Outer surfaces had larger surface areas and weaker E fields, while inner surfaces had smaller surface areas and stronger E fields. They balanced out.

It was only at this point that we began to talk about the formalities of Gaussian surfaces, flux, and surface normals. By the end of the class, the students were able to apply Gauss' Law to a solid sphere of charge, and a solid cylinder of charge. And rather than memorizing the shapes of the surfaces, they were able to reason out why the chosen Gaussian surfaces had their shapes, and how conservation of flux was relevant.

This is not the be all and end all of teaching Gauss' Law to 2nd semester calculus-based physics students. I'm sure there are some other good ways to present the material. But it certainly is the best way I've presented it in the nearly 8 years I've been teaching this class.

Cheers!

XUL is going away

In the past, JSvg was part of OpenTrack which was a Firefox extension to track video data. The first post in this blog shows how to install OpenTrack and then use JSvg. Sadly (or perhaps, not so much, depending on your viewpoint), Firefox is changing their extension architecture. One change that will make it impossible to continue with the status quo, is the deprecation of XUL which is the widget framework upon which the entire browser is based. Another change is the security model which, for the present, still allows signed extensions, but will only allow easy installs from the Firefox Marketplace.

In the first case, since both OpenTrack and JSvg were written using XUL, these changes break the applications unless users choose to maintain an older version of Firefox. While doable, it is generally not recommended for security reasons unless the user only uses Firefox for this particular extension. This generally means setting up separate user profiles for Firefox, and then invoking the browser using the command line (or changing the icon commands) to something like "firefox -P username --no-remote". Although this is relatively easy for the initiated, to do, it also means that if a user wishes to use Firefox for other things, they will need a completely different and up to date installation as well. Many people (students) will not have the knowledge to pull this off.

With regards to signing, Mozilla is following a tighter security model than previously. Not only must the extension be signed (thereby identifying the developer responsible), but also the extension must be distributed by Mozilla after a review. All good, right? Mostly. The issue for instructors such as myself who are sole maintainers is that occasionally students will find bugs that must be addressed immediately. In the past, I have been able to roll out a fix overnight, and when students start up OpenTrack or JSvg the next time, the patch was automagically downloaded and applied via the Firefox extension system. At worst, I would have to tell students to restart their browsers. Now, if I wish to hold to this new distribution model, I need to both completely rewrite these apps (which still basically needs to happen), and then submit them for review. Review of updates to extensions can sometimes take a couple of weeks. Therefore, fixing an issue that might be critical for students can't happen overnight. And for our purposes, this new model won't work very well. Please, don't get me wrong, in general, Mozilla's new way of dealing with security is good.

So what to do? At this point, I have done a minor rewrite of the JSvg graphics and coding portion of the library and hosted it at niiler.com so that an online version is available for students to use. This is now more of a cloud based solution, which seems to be working for my current physics class. The downside at present is that I haven't made it pretty. All the functionality is there (more or less), but the random person stumbling across the page will have no idea what to do. I will deal with this as there is time. In the meantime, go to the above link to continue to use the analysis portion of this library.

Cheers!

Some coming updates

This will be just a short post to show a couple of the new graph types that will be available in the next release of OT/JSvg. Both of them are the result of a project I've been working on for the orthopedics guys at duPont. The first type, a bar graph with error bars, is rather standard, but was not possible using OT/JSvg without all your own extra coding.

It's nothing fancy, but the syntax is rather straight forward:

      drawBarSeriesErrorPlot(xnames, yarr,yerr,fn);

In this case, xnames is a 1d array of labels for the x-axis, yarr is a 2d array of y-values for plotting in series form, yerr is the corresponding 2d array of error values to be added or subtracted from yarr, and fn is the filename to save.

But moving along, we have something, which is to my mind quite a lot more impressive: mapping! Using shapefiles converted from the ogr2ogr utility (used by GRASS and other GIS systems), demographics files from the US census, and then event files (for example number of crimes or accidents), one can then paint a graph as follows.

This map doesn't include events at the current time as that data is not yet public. If you know Delaware, you'll also note that parts of New Castle are missing. That's because the MULTIPOLY parsing isn't working exactly right yet. Other features are an ability to overlay other layers for visualization. So you can draw the map, color the map dependent on the feature, insert points for events, and insert overlay data (gradients, heat maps, etc) to show additional trends.

     drawMap(xCoords, yCoords, refLongs, refLats, overlayData, overlayCT, xPoints, yPoints, coordType, fn, legendStuff) {
    // xCoords - x values or longitude data to be converted (2D)
    // yCoords - y values or latitude data to be converted (2D)
    //      For each series above, a line color and fill type are needed
    // refLong - a reference longitude for drawing map.  Greenwich is far away and 
    //      causes distortion of map over great distances
    // overlayData - 3D array containing data for coloration
    //      This will always be a rectangular array whose upper left position
    //      corresponds to the minimum xCoord and minimum yCoord, and whose
    //      lower right position corresponds to the maximum xCoord and yCoord.
    //      This routine will do the registration of overlayData
    // //overlayCorners - 3D array containing coord data for overlayData
    //      =[ [ [xTop1,xBot1],[yTop1,yBot1] ], ...];
    //      These should be within the xCoords/yCoords range
    // overlayCT - colortable for overlayData 
    // xPoints - x values of points to overlay on top of everything else (2D)
    // yPoints - y values of points to overlay on top of everything else (2D)
    // coordType - 'grid' values or 'longlat' data  if the latter, then convert
    //      For xCoords, yCoords, and xPoints, yPoints, all should be of same type
    // fn - filename for svg file.

For you GIS types, this might not be that impressive, but what comes along with this is a regional ANOVA analysis, Getis-Ord* statistics by point and region, and some frequency analyses by region. It's not quite ready for prime time, at this point, but as I settle on the API, this function and its relatives will make it into the main trunk.

More shortly as there is time.

Cheers!

Snowday data modeling

It's a snowday here in South-Eastern Pennsylvania, so I'm off from work where I teach physics, astronomy and Taijiquan. As the weather is so unpredictable in these parts in the spring semester, we've come up with some computationally based labs to help students better understand the various physics phenomena they've been studying in the event that labs are canceled. Before I go any further, let me say that we attempt to teach the students the basics of JavaScript and give them a number of exercises before going any further.

Today's exercise is a bit more complicated. We're going to generate electric field lines for a distribution of charges. There are a number of ways to do this, but the way we will pick is based on the kinematics of a test charge. In principle, test charges move directly along field lines, and so if one maps out their paths, one also gets a sense of the field line direction. Sounds pretty easy, right? Basically yes, but the devil is in the details. Here is the general algorithm:

1. Make arrays to contain the positions of the charge distribution. These charges are point charges that are fixed in place so we will only need to specify a list of x positions and y positions. Therefore, 1D arrays will due. Call these variables: xpos and ypos, respectively.
2. Make arrays to contain the positions of the test charges. These charges are point charges that move. If we want to draw them, we need to have their positions at a number of different times. Therefore, these will be 2D arrays. The first dimension will refer to the test charge, and the second dimension will refer to the time. Call these variables: x and y.
3. Pick a reasonable time increment dti, and test charge mass masst.
4. Start with initial velocities v0x = 0; and v0y = 0;.
5. Calculate the x and y components of force on a given test charge due to ALL the charges in the charge distribution. Call these values Fx and Fy. These are determined by the Coulomb force.
6. Calculate the accelerations of each of the given test charges based on their masses (all masst. The accelerations will be ax = Fx/masst and ay = Fy/masst.
7. Update the positions of each of the test charges using the regular equations of kinematics. The x and y positions are given by
8. x[charge][time] = x[charge][time - 1] + v0x*dti + 0.5*ax*dti*dti
9. y[charge][time] = y[charge][time - 1] + v0y*dti + 0.5*ay*dti*dti
Remember: we only need to store the x and y positions to draw the field lines.
10. Update the x and y velocities of each test charge:
11. v0x = v0x + ax*dti
12. v0y = v0y + ay*dti
13. Go back to step 5. Rinse, wash, repeat.

In code, we will start by layout out the initial conditions.

var q = 1.0e-6;     // C (magnitude of each charge in charge dist)
var qt = 1.0e-10;   // C (test charge magnitude)
var masst = 1.0e-10;// kg (test charge mass)
var k = 8.99e9;     // N*m^2/C^2  (Coulomb's constant)

// Positions of charge array
var ypos = new Array();
var xpos = new Array();

var N = 3;         // number of charges in charge distrib
var scale = 100;   // used in setting test charge positioning

var extra = 4      // how many MORE test charges than charges in distrib
var Nt = 20;       // Number of time increments
var dt = 5.0e-6;   // seconds (time jump per increment) 

var Ntest = N+2*extra;  // how many test charges 
                        //   2*extra on top AND bottom

Next we set the charge distribution positions. In this case, it's going to be a vertical line at the origin, not quite centered in y.

// These set the charge positions
for (var i=0; i< N; i++) {
    ypos[i] = (i-N/2)/scale;
    xpos[i] = 0;           
}

Having set these up, we then create an array full of our test charges. We'll space them vertically along side our charge distribution, and let them be offset just a bit in x. Why the offset? Because the Coulomb force is infinite at a distance of zero. But before we go on, let's set them according to the scale of the fixed charge distribution. dy is the difference in height between any two charges in this distribution.

var dy = Math.abs(ypos[1]-ypos[2]);   // Difference in charge positions

Now, create the arrays for the test charges, and set an initial position for the lowest test charge. We'll set initial positions for the rest of them shortly.

// These arrays will be 2D and will store the positions of the test charges
// which are a proxy for the field lines.
var x = new Array();    
var y = new Array();

var x0 = 0.005;  // m
var y0 = ypos[0]-extra*dy; // m   

We finally come to the looping which will iterate over test charges and over the number of time increments we chose. First a note on Coulomb's Law. While it is typically written as a 1/r² law, if one applies vector math properly, it is possible to get components without using trigonometry explicitly as indicated in this sketch.

In that case, there is an r³ term in the denominator. The distance between a charge from the distribution and a test charge is given by

            // Distance (x and y) between field point and charge in array
            var deltaX = x[i][j-1] - xpos[p];
            var deltaY = y[i][j-1] - ypos[p];

The distance r³ is then given by this expression

            var r3 = Math.pow(deltaX*deltaX + deltaY*deltaY,1.5);

Then, following the rules of vector math, we (finally) get a form for the Coulomb force as a function of distance:

            //  Fx = k*q1*q2*x/r3 and
            //  Fy = k*q1*q2*y/r3
            Fx += k*qt*q*deltaX/r3;
            Fy += k*qt*q*deltaY/r3;

At this point, you're ready to see the full code for the positioning of the test charges at any time during the simulation. Note that there is also some code to draw arrows.

for (var i=0; i< Ntest; i++) {
    // Loop through Nt test charges, each starting from a point x0 to right of
    // the charge array but at the same height as the charge to the left.
    
    x[i] = new Array(); // Create 2nd dimension of the array
    y[i] = new Array();
    
    y0 += dy;   // m
    
    // Initial positions are known
    x[i][0] = x0;
    y[i][0] = y0;
    
    dti = dt;   //(Math.abs(Ntest/2-i)+0.01)*dt;    
    
    // The test charges start from rest
    var v0x = 0;    // m/s
    var v0y = 0;    // m/s

    // Start loop at j=1 since initial positions are known
    // This is a loop on time increments
    for (var j=1; j< Nt; j++) {
        
        // Now loop through all the charges to determine the Force in x and y
        var Fx = 0;
        var Fy = 0;
            
        // This is a loop on position for a given test charge
        for (var p=0; p< N; p++) {
            // Distance (x and y) between field point and charge in array
            var deltaX = x[i][j-1] - xpos[p];
            var deltaY = y[i][j-1] - ypos[p];
            // r3 = r cubed from Pythagorean theorem ( 3/2 power = 1.5)
            var r3 = Math.pow(deltaX*deltaX + deltaY*deltaY,1.5);
            // These have the form 
            //  Fx = k*q1*q2*x/r3 and
            //  Fy = k*q1*q2*y/r3
            Fx += k*qt*q*deltaX/r3;
            Fy += k*qt*q*deltaY/r3;
        }
        
        var ax = Fx/masst; // x-acceleration of test charge -changes
        var ay = Fy/masst; // y-acceleration of test charge -also changes
                
        x[i][j] = x[i][j-1] + v0x*dti + 0.5*ax*dti*dti;
        y[i][j] = y[i][j-1] + v0y*dti + 0.5*ay*dti*dti;
        
        v0x = v0x + ax*dti;
        v0y = v0y + ay*dti;
 
    }

    if (showArrows) {
        // drawArrowAtEnd(xp,yp,vx,vy,size)
        // x[i][j-1], y[i][j-1], ax, ay are last calculated values in preceding loop.    
        var tmp = drawArrowAtEnd(x[i][j-1],y[i][j-1],v0x,v0y,dy/20);
        var ind = 0;
        for (var jj = j; jj< j+3; jj++) {
            x[i][jj] = tmp[0][ind];
            y[i][jj] = tmp[1][ind];
            ind++;
        }
    }
    
    optsSeries[i].itype=1;        // Be sure it's a line rather than a point.
    optsSeries[i].thick=1;        // Make this 1px thick
    optsSeries[i].stroke="red";   // Make this path red
}

The last bit is code to figure out how to draw arrows at the end of the test charge paths. The inputs are the last calculated point for the test charge, the last velocities, and a length (dy/20) in this case. The output of drawArrowAtEnd() is a 2D array [xpts,ypts] which contains three x and three y points: the x and y coordinates of each side of the arrow, and the point. There is also code here for coloring the paths of the test charges. [Note that drawArrowAtEnd() is not part of OT/JSvg, but is included at the end of this article.

Having done all the calculations, we're now ready to plot. We're going to use our well known drawSeriesPlot() command, so we need to use a 2D array for both x and y. Our test charges are already configured like this, so let's see what happens if we just call our plot command.

Ugh! We certainly don't want a legend, and we can't see the charge distribution. Also, there are no labels. Let's do the easy stuff first and deal with these issues.

optsGraphTitles[0].name = 'x (m)';
optsGraphTitles[1].name = 'y (m)';
optsGraphTitles[2].name = 'Field lines with N = '+N+' charges';

window.optsLegend = false;

OK, this is better, but we also want to add the charge distribution. To do this, we'll just add them to the x and y arrays. Yes, these were originally for the test charges, but now that the calculations are done, we can just add in one more series.

// Also plot the charges

xl = x.length;
x[xl] = xpos;
y[xl] = ypos;

optsSeries[xl].itype=0;            // Make these points
optsSeries[xl].thick=2;            // Make them 2 wide
optsSeries[xl].stroke="blue";      // Color them blue (outline)
optsSeries[xl].color="blue";       // Make the fill blue also

Excellent! Now we are ready to play. Notice how the field diverges? What happens if you add more charges. Let's go with 20 charges, instead of 3. Set N = 3 at the top of the script.

Notice that the field lines in the middle have become more parallel. Want to make the arrows bigger? Change the final parameter of drawArrowAtEnd(). How about dy/2?

Here is the final script used to produce the final graph in its entirety:

/************ Only change things in this section ***********************/

var showArrows = true;

var q = 1.0e-6; // C
var qt = 1.0e-10;   // C (test charge magnitude)
var masst = 1.0e-10;    // kg (test charge mass)
var k = 8.99e9;     // N*m^2/C^2

// Position of charge array
var ypos = new Array();
var xpos = new Array();

var N = 20; // number of charges
var scale = 100;

var extra = 4
var Nt = 20;    // Number of time increments
var dt = 5.0e-6;    // seconds

var Ntest = N+2*extra;

/***********************************************************************/


function drawArrowAtEnd(xp,yp,vx,vy,size) {
    // xp and yp are the end point coordinates
    // vx and vy are vector components - used to get the angle
    // size is the size in x-y space of the arrow end points
    
    // The arrow is drawn by drawing points starting at the end point
    // and doing the lower "fletching", moving back to the end point and
    // then doing the upper "fletching".
    
    var r2d = 180/Math.PI;
    var d2r = 1/r2d;
    
    // calculate the slope of the vector
    if (xp != 0) {
        if (yp/xp < 1000) {
            // To ensure sanity
            var theta0 = r2d*Math.atan(yp/xp);
            var theta1 = (theta0-135)*d2r;
            var theta2 = (theta0+135)*d2r;
        } else { // do the same thing as if xp = 0
            if (yp < 0) {   // It points down
                theta1 = 135*d2r;
                theta2 = 45*d2r;
            } else {    // It points up
                theta1 = 45*d2r;
                theta2 = 135*d2r;
            }
        }

    } else {
        if (yp < 0) {   // It points down
            theta1 = 135*d2r;
            theta2 = 45*d2r;
        } else {    // It points up
            theta1 = 45*d2r;
            theta2 = 135*d2r;
        }    
    }
    console.log(theta1,theta2)
    var dx1 = size*Math.cos(theta1);
    var dy1 = size*Math.sin(theta1);
    var dx2 = size*Math.cos(theta2);
    var dy2 = size*Math.sin(theta2);

    var xout = [xp+dx1,xp,xp+dx2];
    var yout = [yp+dy1,yp,yp+dy2];
    return [xout,yout]; // 2d array containing arrow
}


// These set the charge positions
for (var i=0; i< N; i++) {
    ypos[i] = (i-N/2)/scale;
    xpos[i] = 0;           
}

var dy = Math.abs(ypos[1]-ypos[2]);   // Difference in charge positions

// These arrays will be 2D and will store the positions of the test charges
// which are a proxy for the field lines.
var x = new Array();    
var y = new Array();

var x0 = 0.005;  // m
var y0 = ypos[0]-extra*dy; // m   


// Now start a test charge at a distance of x0 from the array of charges and 
// see where it goes.  This will indicate the direction of the field line.  To
// do this, one must first calculate the net field at the point in order to get
// the force on the charge and thereby its acceleration.
// Use equation x(i+1) = x(i) + v0x*ttot + 0.5*ay(i)*ttot^2 (pseudo code)
// Ditto with y(i+1)
// In this case the particle is likely to go VERY fast.  Choose a very small dt:

for (var i=0; i< Ntest; i++) {
    // Loop through Nt test charges, each starting from a point x0 to right of
    // the charge array but at the same height as the charge to the left.
    
    x[i] = new Array(); // Create 2nd dimension of the array
    y[i] = new Array();
    
    y0 += dy;   // m
    
    // Initial positions are known
    x[i][0] = x0;
    y[i][0] = y0;
    
    dti = dt;   //(Math.abs(Ntest/2-i)+0.01)*dt;    
    
    // The test charges start from rest
    var v0x = 0;    // m/s
    var v0y = 0;    // m/s

    // Start loop at j=1 since initial positions are known
    // This is a loop on time increments
    for (var j=1; j< Nt; j++) {
        
        // Now loop through all the charges to determine the Force in x and y
        var Fx = 0;
        var Fy = 0;
            
        // This is a loop on position for a given test charge
        for (var p=0; p< N; p++) {
            // Distance (x and y) between field point and charge in array
            var deltaX = x[i][j-1] - xpos[p];
            var deltaY = y[i][j-1] - ypos[p];
            // r3 = r cubed from Pythagorean theorem ( 3/2 power = 1.5)
            var r3 = Math.pow(deltaX*deltaX + deltaY*deltaY,1.5);
            // These have the form 
            //  Fx = k*q1*q2*x/r3 and
            //  Fy = k*q1*q2*y/r3
            Fx += k*qt*q*deltaX/r3;
            Fy += k*qt*q*deltaY/r3;
        }
        
        var ax = Fx/masst; // x-acceleration of test charge -changes
        var ay = Fy/masst; // y-acceleration of test charge -also changes
                
        x[i][j] = x[i][j-1] + v0x*dti + 0.5*ax*dti*dti;
        y[i][j] = y[i][j-1] + v0y*dti + 0.5*ay*dti*dti;
        
        v0x = v0x + ax*dti;
        v0y = v0y + ay*dti;
 
    }

    if (showArrows) {
        // drawArrowAtEnd(xp,yp,vx,vy,size)
        // x[i][j-1], y[i][j-1], ax, ay are last calculated values in preceding loop.    
        var tmp = drawArrowAtEnd(x[i][j-1],y[i][j-1],v0x,v0y,dy/2);
        var ind = 0;
        for (var jj = j; jj< j+3; jj++) {
            x[i][jj] = tmp[0][ind];
            y[i][jj] = tmp[1][ind];
            ind++;
        }
    }
    
    optsSeries[i].itype=1;
    optsSeries[i].thick=1;
    optsSeries[i].stroke="red";
}

// Also plot the charges

xl = x.length;
x[xl] = xpos;
y[xl] = ypos;

optsSeries[xl].itype=0;
optsSeries[xl].thick=2;
optsSeries[xl].stroke="blue";
optsSeries[xl].color="blue";

optsGraphTitles[0].name = 'x (m)';
optsGraphTitles[1].name = 'y (m)';
optsGraphTitles[2].name = 'Field lines with N = '+N+' charges';

window.optsLegend = false;
/*
window.xprec = 4;
*/
// OK, let's plot this stuff

drawSeriesPlot(x,y,0,0,'fieldLines.svg');

Until next time: cheers!