emfieldcalc2.mws

Program 1- Maple Computer Program to Calculate the Classical Fields of a Moving Charge.

(A Theoretical Retarded Potential Method Using Lienard-Wiechert Potentials)

>

BRIEF THEORY OF THEORETICAL METHOD:-

The formula for the field from a relativistic charge and its derivation is given in many textbooks on electromagnetism. The exact form in which the equation is expressed can vary although the calculated field will of course be the same. The technique of retarded potentials is particularly well explained in Schwartz, Principles of Electrodynamics (ref 7) which covers this topic in considerable detail. The electric field E(rmeas,t) at the measuring point (rmeas) and time (t) can be separated into a radiation field and an induction field The induction field is independent of the charge acceleration and contains a term proportional to whereas the radiation field is proportional to the charge acceleration and only has a term so decays more slowly. The field equations use charge parameters at the retarded position (rret) and at the retarded time (tret). The retarded position being the earlier position of the charge where the field which left it at the speed of light and at the retarded time has just reached the measuring point at time t. The distance between the retarded position and the measuring point is and therefore is the time light takes to travel this distance. The retarded time, at which light left the charge is therefore given by the equation:- The formula for the radiation electric field in Schwartz is in cgs units but expressing this in rationalised MKSA units gives:- Where q is the charge magnitude in Coulombs.

vret is the charge velocity vector at the retarded position in metres/sec.

aret is the charge acceleration vector at the retarded position in metres/sec/sec

c is the velocity of light in metres/sec - a scalar quantity ie just a magnitude. is a unit vector from the retarded position (rret) at time tret to the measuring point (rmeas) at time t.

. is a vector dot product.

x is a vector cross product.

e0 is the permittivity of free space.

To find the unit vector (epsilon) we divide the vector from the retarded position to the measuring point by its length ie:- The formula for the induction electric field Eind(rmeas,t), in volts/metre at the measuring point rmeas and time t in MKS units is:- Where is not a vector but just the scalar magnitude of the charge velocity in metres/sec at the retarded position.

The H field in ampere-turns/metre can be found using a similar formula but the simplest method is to derive it from the E fields. The required formula is difficult to find expressly stated in MKS units but can be seen in Stratton (ref 8) to be equivalent to:- These equations should give the correct fields for a charge travelling at any velocity from stationary up to almost the speed of light. A program was written using Maple 6 as this provides a relatively painless method of calculating the fields from the above formulae as shown below. As frequently happens with symbolic mathematics programs some of the expressions resulting from an intermediate calculation are quite long and complicated and to save space these have been suppressed in the following display as the main interest is in the final values provided for the field magnitudes. For anyone able to run the program the intermediate formulae can easily be examined.

For this program to work you must have Maple software installed on your computer and first manually create a folder c:/field on the c: drive of your pc. If an alternative path is to be used then all references to this path must be changed in the program. All user selectable values are indicated by starting the side comment with a triple asterisk (***). The program as shown has initially been set up to calculate the field from an electron rotating in a circular orbit around the origin with a radius of 0.7*10^(-15) metres and a velocity (gamma) of 137. By varying the parameters any radius and any charge velocity can be used instead. It can easily model a number of charges, either stationary, moving in a straight line (linear motion), or oscillating backwards and forwards (harmonic oscillators) or combinations of all these. Although the program correctly calculates the classical field values from such arrangements it does not check that such a charge motion is possible nor does it take into account any 'quantum' effects or electron spin which may be important at small distances. Any other charge path could also be modelled provided an equation of motion is provided. In order to provide this flexibility there are a number of different options which must be selected correctly for the program to work. For this example all of sections 1.1, 1.2, 1.3, 1.4B, 1.5, 1.6, 1.7,1.9 and 1.10 should be entered. These options are probably rather difficult to follow initially so please read the comments carefully.

>

1.1 SET THE TEST NUMBER AND NUMBER OF DIGITS TO USE:-

> restart;

> testnum:=1; *** This is used as the main identification number for stored files. Data for each test will be stored in a separate file called "c:/field/params"||testnum||"-1.m". Increment testnum for each different simulation. If testnum is kept the same previous results will be partly overwritten and all data could be ruined. Files for archive can be periodically moved from the c:/field folder manually and stored in a separate folder. A program for checking details of these files stored in the field folder is included at the end of this page. > Digits:=60; *** Digits needs to be very high when calculating the fields near a very highly relativistic charge for the Maple fsolve command to find a solution in some situations. > digits:=Digits; This enables the value assigned to Digits to be stored later in the program . 1.2 DEFINE THE VALUE OF COMMON CONSTANTS USED BY THE PROGRAM:-

> u0:=4*Pi*10^(-7); Permeability of a vacuum. > c:=299792458; Velocity of light in a vacuum (metres/second). > e0:=evalf(1/(c^2*u0)); Permittivity of a vacuum. > dl:=10^(-15); ***Enter the scaling factor dl. Dimensions will be dl times the number entered. Suggest a multiple or sub multiple of ten selected. To avoid the complications of a scaling factor just make dl=1 and it can then be ignored. If so note that if any value is less than 10^(-15) then Maple 6 will not plot it - but another scaling factor can then be used just for plotting if needed.. >

1.3 DETERMINE THE EQUATIONS OF MOTION OF THE CHARGE:-

> N:=1; N is the number of charges. Assume just one charge initially. This can be altered later if more charges are required. > t:='t'; a:='a'; f:='f'; gamm:='gamm'; This ensures these variables are not initially assigned any value.    Enter the charge position equations for x, y and z coordinates in metres. Typical entries would be as follows where dist, dist1 and dist2 (in metres/dl) and v (in metres/sec) must be numeric values . Other entries should be entered as symbols as the program will evaluate them:-

dist*dl for a stationary charge. Where dist*dl is the distance along the co-ordinate axis in metres.

dist1*dl+v*t for a charge moving with linear motion where v is the charge velocity (metres/second) along the co-ordinate axis and dist*dl is the charge position in metres at t=0.

dist1*dl+a*dl*sin(2*Pi*f*t) in one dimension for a harmonic oscillator. Where dist*dl gives the numerical distance of the centre of oscillation along the co-ordinate axis and a*dl will be the peak displacement (in metres).

dist1*dl+a*dl*sin(2*Pi*f*t) in one dimension and dist2*dl+a*dl*cos(2*Pi*f*t) in the other dimension for an orbiting charge. Where a*dl is the radius of the orbit and dist1*dl and dist2*dl are distances along the co-ordinate axis to the centre of orbit.(dist2=-a for z coordinate to position the charge at x=0,z=0 coordinates at t=0).

> Xqpos:=dl*a*sin(2*Pi*f*t); ***In this case the equation of motion for a charge in a circular orbit around the origin. > Yqpos:=0; *** Enter the y axis position equation > Zqpos:=dl*a*cos(2*Pi*f*t); *** Enter the z axis position equation Find the charge velocity by differentiating the charge position equations with respect to time:-

> VXqpos:=diff(Xqpos,t); The charge velocity in metres/sec. > VYqpos:=diff(Yqpos,t); > VZqpos:=diff(Zqpos,t); To find an expression for the maximum charge velocity, enter the time at which the velocity is a maximum:-

> Tmaxv:=0; ***In this case velocity is constant so enter Tmaxv:=0. A different value may be needed for another set of equations of motion. > t:=Tmaxv; > Vmax:=sqrt(VXqpos^2+VYqpos^2+VZqpos^2); This is the maximum velocity. The formula for gamma corresponding to this maximum velocity is:-

> eqn1:=gamm=1/sqrt(1-Vmax^2/c^2); If the value of gamma has been numerically evaluated, which it will for a stationary charge or one in linear motion, check it looks correct (gamm should be a real number between 1 and just less than infinity). >

>

>

>

>

1.4 THE FOLLOWING OPTIONS ALLOW VARIOUS METHODS OF FINDING PARAMETERS a AND f IN THE EQUATIONS OF MOTION AND IS REQUIRED FOR OSCILLATING OR ORBITING CHARGES ONLY. For stationary charges or ones in linear motion jump to section 1.7:-

NB. Remove the hash symbol (#) at the beginning of a line if it is required to run it.

SELECT EITHER A or B or C BELOW:-

1.4A) To select a desired frequency and a maximum value of gamma and then find the corresponding displacement:-

> #f:=1; #438.5*10^6; *** Enter a desired charge oscillation frequency.

> #gamm:=10; #137; *** Enter a desired value for gamma.

The corresponding value for the maximum displacement is now obtained by solving eqn1 for a:-

> #res:=solve(eqn1,a);

> #a:=abs(res);

>

1.4B) OR ALTERNATIVELY to select a displacement and maximum value of gamma and solve for the frequency:-

> a:=0.7; *** Enter the displacement in metres/dl (remember the number entered is automatically multiplied by the scaling factor dl in the equation of motion). > gamm:=137; *** Enter a desired value for gamma. The corresponding value for the frequency is now obtained by solving eqn1 for f :-

> res:=solve(eqn1,f);  > f:=abs(res); Select the positive root from the above expression. (% in Maple just calls the last result) >

1.4C) OR ALTERNATIVELY to select a displacement and maximum velocity in metres/second and solve for the frequency:-

> #a:=0.7; *** Enter the displacement in metres (remember the number entered is automatically multiplied by the scaling factor dl).

> #v:=2.185541354*10^6; ***Enter the desired maximum velocity in metres/second.

To find the corresponding value of gamm:-

> #gamm:=1/sqrt(1-v^2/c^2);

The corresponding value for the frequency is now obtained by solving eqn1 for f :-

> #res:=solve(eqn1,f);

> #f:=abs(res); Select the positive root from above.

>

>

1.5 CHECK THE VALUES THAT HAVE BEEN ASSIGNED TO f AND a TO CONFIRM THEY ARE CORRECT (only applies to oscillating and orbiting charges):-

We should now have the values for f and a. The method of obtaining them is a little laboured but has the advantage that the peak velocity can not exceed the velocity of light. If instead a desired frequency and displacement are selected then it would be necessary to check that the maximum velocity is less than c otherwise strange field values will be obtained....as might be expected!! Next is a final check on the values selected:-

> f; a;  1.6 TO PLOT THE CHARGE MOTION:-

NB This is not always required and there are several different plotting methods but this one often works OK for oscillatory motion. (This plot method is unsuitable for a stationary charge or one with linear motion)

> with (plots):

```Warning, the name changecoords has been redefined

```
` `

> plotsetup(inline);

> t:='t'; > plot([Xqpos/dl,(Zqpos/dl),t=0..1/f],scaling=constrained,colour=red,labels=["X","Z"]); Note that lengths are divided by dl. 1.7 ENTER THE CHARGED PARTICLE MAGNITUDE AND MASS.

> q:=-1.601864*10^(-19); *** Enter charge values. For an electron q:=-1.601864*10^(-19)

> m:=9.1*10^(-31); *** Enter charge masses.For an electron m:=9.1*10^(-31). This is not needed to calculate field magnitudes but may be useful in future programs.  >

>

>

1.8 FOR A SECOND CHARGE DIAMETRICALLY OPPOSITE (Skip this section if there is only one charge);-

For N charges sections 1.3, 1.4, 1.5, 1.6 and 1.7 above could be copied N times and position and velocity equations Xqpos[N],Yqpos[N],Zqpos[N] and VXqpos[N], VYqpos[N], VZqpos[N] found using variables a[N], f[N], gamm[N],Tmaxv[N],Vmax[N],q[N],m[N] etc. A common time t should be used.

NB. This section (1.8) is a simplified version for a second charge which is ideal for simple simulations such as two orbiting charges and uses the same values of a and f used above.

> N:=2; *** Enter the total number of charges.

> t:='t'; Ensures t is unassigned.

> Xqpos:=-dl*a*sin(2*Pi*f*t); ***In this case the equation of motion for the position of a charge orbiting diametrically opposite the first one

> Yqpos:=0; *** Enter the y axis position equation

> Zqpos:=-dl*a*cos(2*Pi*f*t); *** Enter the z axis position equation

Find the charge velocity by differentiating the charge position equations with respect to time:-

> VXqpos:=diff(Xqpos,t); The charge velocity in metres/sec.

> VYqpos:=diff(Yqpos,t);

> VZqpos:=diff(Zqpos,t);

>

>

>

In this case, for second charge and assuming q is a positron:-

> q:=1.601864*10^(-19); *** Enter charge values. For a positron q:=1.601864*10^(-19)

> m:=9.1*10^(-31); *** Enter charge masses. For a positron m:=9.1*10^(-31)

>

>

>

>

1.9 TO CALCULATE THE RELATIVISTIC FIELD EQUATIONS:-

This section is the heart of the program which calculates the field equations. NB. It would be possible to replace all the above programming with just a list to enter the value of constants q[nq],e0, u0 & c. Set Digits & N and enter the correct charge equations of motion Xqpos[nq], Yqpos[nq] & Zqpos[nq] which are differentiable functions of t. This would be sufficient to evaluate the field values. To save the test data details, testnum and digits are really required to be assigned values as well.

> with(linalg):

```Warning, the protected names norm and trace have been redefined and unprotected

```
` `

> t:='t'; tret:='tret'; tact:='tact'; This ensures no value is initially assigned to t or tret or tact.   Finally check the position equations of the charges is correct. In metres these are:-

> Xqpos; Yqpos; Zqpos;     > Xqpos; Yqpos; Zqpos;   > N; >

>

> for nq from 1 to N do The whole of the next field equation generation block will be run N times so that a complete set of equations is produced for the field from each charge.

>

The retarded charge position vector (rret) is just the vector defining the charge position, but with the time now being the retarded time (tret) and distances still in metres ie:-

> t:=tret[nq]:

> rret[nq]:=convert([Xqpos[nq],Yqpos[nq],Zqpos[nq]],vector):

>

> t:='t'; This unassigns t.

Let the vector defining the measuring position be:-

> rmeas:=convert([xmeas,ymeas,zmeas],vector):

The x distance from the charge retarded position to the measuring point is obtained using:-

> evalm(rmeas-rret[nq]):

The unit vector from the charge retarded position to the measuring point is:-

> epsilon[nq]:=evalm((rmeas-rret[nq])/(sqrt(evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2))):

The charge velocity at the retarded position is obtained by differentiating the charge retarded position vector (rret) with respect to retarded time tret:-

> vret[nq]:=convert([diff(rret[nq],tret[nq]),diff(rret[nq],tret[nq]),diff(rret[nq],tret[nq])],vector):

The charge acceleration at the retarded position is similarly obtained by differentiating the retarded velocity (vret) with respect to the retarded time once more:-

> aret[nq]:=convert([diff(vret[nq],tret[nq]),diff(vret[nq],tret[nq]),diff(vret[nq],tret[nq])],vector):

The denominator of the radiation field is:-

> denominator1[nq]:=evalm(c^2*(sqrt(evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2))*(1-dotprod(vret[nq],epsilon[nq]/c,'orthogonal'))^3):

The complete formula for the radiation field is:-

The induction field has components:-

> denominator2[nq]:=evalm((evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2)*(1-dotprod(vret[nq],epsilon[nq]/c))^3):

The complete expression for the induction field is:-

> Eind[nq]:=evalm((q[nq]/(4*Pi*e0)*(epsilon[nq]-vret[nq]/c)*(1-(vret[nq]^2+vret[nq]^2+vret[nq]^2)/c^2))/denominator2[nq]):

The total combined radiation and induction fields for a single charge (numbered nq) are:-

The magnetic field is given by:-

For the x,y and z total H field components:-

> Hxt[nq]:=(Ht[nq]):

> Hyt[nq]:=(Ht[nq]):

> Hzt[nq]:=(Ht[nq]):

>

>

To find phi, the scaler potential:-

> St[nq]:=q[nq]/(4*Pi*e0*sqrt(evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2)*(1-dotprod(vret[nq],epsilon[nq]/c,'orthogonal'))):

To find A, the vector potentials:-

> Axt[nq]:=St[nq]*vret[nq]/c^2:

> Ayt[nq]:=St[nq]*vret[nq]/c^2:

> Azt[nq]:=St[nq]*vret[nq]/c^2:

>

>

The equation relating the retarded time and the actual time is:-

> eqn2[nq]:=tret[nq]=t-(sqrt(evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2+evalm(rmeas-rret[nq])^2))/c;

rret[nq] in eqn2 also depends on tret[nq] so solve the above eqn2[nq] for tret[nq]:-

> tretsol[nq]:=solve(eqn2[nq],tret[nq]); In this line Maple attempts to find a theoretical solution for tret. Although often a complicated expression this is usually possible for a stationary charge, linear motion and an harmonic oscillator. In this case it would be possible to obtain solved equations for field values directly. For orbital motion Maple 6 solver fails. Because of this in the following field calculation programs this theoretical solution is not used and instead a numerical solution is found with the Maple fsolve command. This can usually find a solution whatever the charge path provided Digits has been set high enough

> print("Maple Equation for tret is tret["||nq||"]= ",tretsol[nq]); This will print the formula for the solution if Maple can find one. As stated above, however, this solution is not normally used.

>

> end do: >

1.10 CHECK RESULTS, ENTER DESCRIPTION OF TEST AND STORE THIS CHARGE SIMULATION:-

Check the equation for tret (should be a function of tret1, t, xmeas, ymeas and zmeas):-

> tret:='tret'; > eqn2;    If a second charge has been entered eqn2 should be a function of tret2, t, xmeas, ymeas and zmeas (and similarly for the Nth charge) :-

> eqn2; Check the value of gamm:-

> gamm; Enter test description:-

> details:="One negative charge in a circular orbit in x/z plane, gamma=137, dl = 10^(-15).Starting at x=0 z=0.7*10^(-15) at t=0.Digits=60.To find field. Derived using emfieldcalc2.mws"; ***Alter as appropriate. >

Store the above equations and data for future use in plotting programs etc.:-

>

>

To Obtain Field Values

Having stored the above equations they can be used to obtain field values at any measuring point defined by coordinates xmeas, ymeas and zmeas and with the charge positions defined by the value assigned to time t. A variety of plotting programs can be written to enter these parameters and obtain and plot the field values automatically. About the simplest program possible for calculating the resultant field from either one or many charges is listed below and uses manual entry of the input data. The corresponding field contributions from each charge are calculated and summed and the total field shown.

> restart;

> testnum:=1; ***Check testnum is correct and alter if necessary. This will be used to select the charge data file produced above. > read "c:/field/params"||testnum||"-1.m": This reads the stored field equations and test parameters

> Digits:=digits; Alter digits to the stored value. Check the test details are as expected:-

> details; To obtain field values just enter the required co-ordinates of the field measuring point, in metres, remembering to use the scale factor (dl):-

> xmeas:=0; ***Field measuring point x co-ordinate. > ymeas:=0; ***Field measuring point y co-ordinate. > zmeas:=-a*dl; ***Field measuring point z co-ordinate. In this case in the path of the orbiting charge and 180 degrees from it. Enter the time at which the field is to be calculated (this will determine the charge location(s)):-

> t:=0; ***Could be t=0 or other possible times to enter might be a fraction of 1/f. >

Run the field calculation routine:-

First set the combined field stores to zero

> Extc:=0: Hxtc:=0: Eradxc:=0: Eindxc:=0: Axtc:=0: Stc:=0:

> Eytc:=0: Hytc:=0: Eradyc:=0: Eindyc:=0: Aytc:=0:

> Eztc:=0: Hztc:=0: Eradzc:=0: Eindzc:=0: Aztc:=0:

> for nq from 1 to N do This loop selects each charge in turn.

>

> tret:='tret'; Ensure no values assigned to tret, tret etc

> tretans[nq]:=fsolve(eqn2[nq],tret[nq]); Solve for the retarded time (tret) for this charge

> if type(tretans[nq],numeric) then Check that Maple has found a numerical solution for tret.

> tret[nq]:=tretans[nq];

> print("Retarded Time tret["||nq||"] = ",tret[nq]);

>

> else If no solution found print error message.

>

> ERROR("SORRY BUT MAPLE FSOLVE IS UNABLE TO FIND A SOLUTION FOR TRET IN EQUATION ",tretans[nq],"TRY A HIGHER VALUE FOR DIGITS"); NB. If a higher value for Digits is needed all data on this page must be re-entered and all equations re-calculated and saved.

> end if;

>

> Extc:=Re(evalf(Ext[nq]+Extc)); Sum the electric fields from each charge.

> Eytc:=Re(evalf(Eyt[nq]+Eytc));

> Eztc:=Re(evalf(Ezt[nq]+Eztc));

> Hxtc:=Re(evalf(Hxt[nq]+Hxtc)); Sum the magnetic fields from each charge.

> Hytc:=Re(evalf(Hyt[nq]+Hytc));

> Hztc:=Re(evalf(Hzt[nq]+Hztc));

> Eindxc:=Re(evalf(Eind[nq]+Eindxc)); Sum the induction fields.

> Eindyc:=Re(evalf(Eind[nq]+Eindyc));

> Eindzc:=Re(evalf(Eind[nq]+Eindzc));

> Axtc:=Re(evalf(Axt[nq]+Axtc)); Sum the vector potentials.

> Aytc:=Re(evalf(Ayt[nq]+Aytc));

> Aztc:=Re(evalf(Azt[nq]+Aztc));

> Stc:=Re(evalf(St[nq]+Stc)); Sum the scalar potential.

>

> end do:

>

Print the results

>

> print("Total Ex Field = ",Extc," Volts/metre");

> print("Total Ey Field = ",Eytc," Volts/metre");

> print("Total Ez Field = ",Eztc," Volts/metre");

> print("Total E Field = ",sqrt(Extc^2+Eytc^2+Eztc^2)," Volts/metre");

> print("Total Hx Field = ",Hxtc," Ampere-turns/metre");

> print("Total Hy Field = ",Hytc," Ampere-turns/metre");

> print("Total Hz Field = ",Hztc," Ampere-turns/metre");

> print("Total H Field = ",sqrt(Hxtc^2+Hytc^2+Hztc^2)," Ampere-turns/metre");

> print("Total Ex Induction Field = ",Eindxc," Volts/metre");

> print("Total Ey Induction Field = ",Eindyc," Volts/metre");

> print("Total Ez Induction Field = ",Eindzc," Volts/metre");

> print("Total Ax Vector Potential = ",Axtc," Weber/metre");

> print("Total Ay Vector Potential = ",Aytc," Weber/metre");

> print("Total Az Vector Potential = ",Aztc," Weber/metre");

> print("Total Scaler Potential = ",Stc," Volts");

>

>

>

> # print("xmeas=",xmeas,"ymeas=",ymeas,"zmeas=",zmeas);

>

>                   >

>

NB If you just get Float(undefined) for all values on first running this program you have probably entered details for the second charge and so are measuring the field at the location of the second charge where the field will be infinite. In this case re-enter all the program again but ensure that charge 2 is omitted.

>

Interestingly the field of a stationary charge at the same location is :-

> evalf(-q/(4*Pi*e0*(a*dl*2)^2)); The minus sign is required as this formula does not take into account the field direction, which will be the positive z direction. This is very similar to the total field of the moving charge, although the direction is approx 47.6 degrees clockwise:-

> sqrt(Extc^2+Eztc^2); Total field of moving charge at measuring point. > evalf(360/(2*Pi)*arctan(Extc/Eztc)); Field angle at measuring point, in degrees clockwise from charge direction (ie from z axis). >

TO CHECK THE DETAILS OF THE STORED TESTS:-

NB If the description is rubbish enter better details next time!!

> testrange:=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,60,61,70,71,80,81,100,101,102,103,104,105,110,111,115,116,117,118,119,120; ***Enter the testnum number of file details required (eg testrange:=1,2,3,4,5,6; )  The following code block prints the test details. NB If a test involves a lot of charges (ie several thousand) then reading the file can take some time so please wait.

> for tests in testrange do

> try

>

This code allows the program to continue if a file is missing:-

> catch "could not open":

> next:

> end try;

For a file which is present print the test details:-

> print("TEST "||tests||" DETAILS ARE:- "||details);

> end do:

>

>

> 