- 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.
- 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.
- Pick a reasonable time increment dti, and test charge mass masst.
- Start with initial velocities v0x = 0; and v0y = 0;.
- 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.
- 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.
- Update the positions of each of the test charges using the regular equations of kinematics. The x and y positions are given by
- x[charge][time] = x[charge][time - 1] + v0x*dti + 0.5*ax*dti*dti
- 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.
- Update the x and y velocities of each test charge:
- v0x = v0x + ax*dti
- v0y = v0y + ay*dti
- 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/r2 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 r3 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 r3 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!
No comments:
Post a Comment