1. This book discusses how to compute the fields of a point charge, moving in any known fashion.
  2. It is clear that Maxwell originally conceived of electric charge as a continuous fluid. In Maxwellian theory the amount of charge that can occur in a finite volume ranges from zero upward. With the discovery of the electron and other charged particles, however, it became clear that electric charge always occurs in quanta of +1.6E-19 coulombs. More generally, wherever charge is present it seems always to be an attribute of a "fundamental" particle. (And such a particle may also have mass, angular momentum, etc.) With the discovery of the neutron it appeared that some particles are uncharged (although the modern view is that they have internal charges ... quarks ... but zero net charge). In any case it proved convenient to conceive of electric charge as being concentrated at points in space. The charge density would be infinite at such points, and zero elsewhere. Maxwell’s paradigm of continuously varying charge density would then amount to the average density over regions of space containing many point charges.

    1. Definition of a Point Charge
    2. Strictly speaking a point charge should be defined to invariably consist of +1.6E-19 coulombs at a point in space. However, common usage is to imagine that a point charge can consist of any quantity of charge (e.g. 1 coulomb), all existing at a point in space and surrounded by charge-free space. This is the concept that will be used in this document.

    3. Field Sources
    4. A source charge by definition engenders an electromagnetic field that can cause other charges to experience external electromagnetic (Lorentz) forces. For example, if source charge q1 engenders field vectors E1 and B1 at the location of charge q2, then the Lorentz force experienced by q2 is

      . (1.2_1)

    5. Field Objects
    6. Field objects are defined herein to be point charges that experience Lorentz forces as a consequence of E and B fields at their positions. The field vectors may be external in nature, as in Eq. 1.2_1. But in addition, certain motions of a field object may result in its own time-varying fields inducing nonzero fields right at its location. When this occurs the object may also experience an internal or self force. The total Lorentz force experienced by an object charge is generally the sum of the external and the internal Lorentz forces.

      1. An Example of a Self Force

      A point charge traveling in a circle experiences a tangential force that acts opposite to the charge’s velocity vector. If this tangential force is not counteracted by an appropriate external "agent" force, the charge will not travel in a circle at constant speed. In the absence of a counteracting agent force, the charge will tangentially decelerate and spiral in towards the circle’s center. If the tangential force is counteracted by a tangential agent force, on the other hand, the charge will persist in traveling in a circle at constant speed. Note that in this monograph "agent" forces will be considered to be non-electromagnetic contact forces. Although non-Lorentzian, they are implicit in Maxwellian theory as, for example, when a constraining agent holds two point charges at a constant distance from each other.

    7. Difficulties with the Point Charge Concept

    The electric field energy of a resting, point charge is infinite. (This fact begs the question of how a point charge could actually exist in nature.) It is also true that the electromagnetic mass (defined below) of a true point charge would be infinite. (It would require infinite forces to accelerate such a charge.) Notwithstanding such difficulties, the concept of point charges has enjoyed a prominent place in electromagnetic theory.

  3. Spherical Shells
  4. Perhaps in order to avoid some of the theoretical difficulties mentioned above, a point charge is often approximated as a spherical shell of minuscule radius, R. At distances of a few R from the shell’s center, the fields are practically indistinguishable from those of a point charge. It is often pointed out that a point charge is the limit of a spherical shell as R shrinks to zero.

    1. Field Energy
    2. In free space the electric field energy density is

      =(2.1_1)

      Thus integration suggests that the field energy of a spherical shell is

      (2.1_2)

    3. Field Momentum
    4. The momentum density in the electromagnetic field is

      g=(2.2_1)

      Integration produces the result that the total momentum in the electromagnetic field of a spherical shell, moving at constant, non-relativistic velocity v, is

      (2.2_2)

      In light of this, the electromagnetic (rest) mass of such a shell is defined to be

      =(2.2_3)

    5. Poincare Stresses
    6. According to the mass/energy equivalence principle, E=mc2. Thus the field energy in Eq. 2.1_2 is less than memc2. In order to understand this inequality, the mechanism whereby a spherical shell might be formed must be considered. The idea is that increments of charge are brought in from infinity and the charged shell is built up, one increment at a time. However, an increment of charge must not only be brought in to a distance R from the shell’s center, but the charge that is already there must be squeezed aside to make room for the newly arrived increment. And in order to do this, stresses in the surface of a shell must be overcome. These stresses are referred to as Poincare stresses. Theoretically the energy that must be expended to overcome these stresses, plus the work that must be done to bring charge increments in from infinity, sum to memc2. In general Poincare stresses are inherent attributes of any Maxwellian, continuous distribution of charge. Their energy exists in addition to the field energy implied by Eq. 2.1_2.

    7. Kinetic Energy

    A moving charge’s kinetic energy also resides in its electromagnetic field. More specifically it resides in the charge’s magnetic field. If one integrates the magnetic field energy density of a spherical shell over all of space, one obtains memv2/2. In effect, as the spherical shell is accelerated, energy flows from the accelerating agent out into the magnetic field. And as it is decelerated, energy flows back in from the magnetic field to the decelerating agent. In both cases the power emitted/absorbed from the field theoretically equals the rate at which the accelerating/decelerating agent does work to accelerate/decelerate the charge.

    The stress energy in the Poincare stresses also increases as a spherical shell is accelerated. This increase can be attributed to the manner whereby the shell length contracts with increasing speed. The total energy of the moving shell (memc2) is the sum of the magnetic and stress energies.

  5. Point Charge Field Solutions
  6. Given knowledge of a point charge’s past motion, the E and B fields at any point in space (other than that occupied by the charge itself) can be computed by combining certain retarded quantities in the proper equations. These equations are relativistically rigorous; they provide the correct values for the field vectors at all source charge speeds.

    1. Retarded Time
    2. The retarded time (tr) for a given point charge varies from field evaluation point to field evaluation point. The retarded time is defined to be the time (in the past) when a light-speed signal would have had to leave the source charge in order to arrive at field evaluation point (x,y,z) at the present time, t. In the case of point source charges, the retarded time is unique, as are the retarded position, velocity, and acceleration. (The retarded position was the source charge’s position at time tr. Similarly the retarded velocity is the source charge’s velocity at time tr, etc.)

    3. General Solutions

    Let D be the vector from a source charge’s retarded position to a field evaluation point. The magnitude of D is c(t-tr). Define utility vector, u, to be

    (3.2_1)

    Then E and B at the field evaluation point are

    (3.2_2)

    (3.2_3)

    Note the use of retarded quantities in these equations. A simple computer routine can be used to closely approximate these quantities. The algorithm below suffices for source charge motion confined to the x-axis and for a field evaluation point lying on the x-axis. We define dt to be t-tr, the time between the present and the retarded time. dtmin and dtmax are defined to be minimum and maximum trial values of dt. dtmin is initialized to zero, and dtmax to twice an impossibly large value for dt. The following algorithm produces the retarded quantities:

    do

    dt=(dtmax+dtmin)/2

    tr=t-dt

    xr=x(tr)

    Dr=x-xr

    if abs(c*dt-Dr)<2^(-30) then exit do

    if c*dt-Dr>0 then

    dtmax=dt

    else

    dtmin=dt

    endif

    loop

    vr=v(tr)

    ar=a(tr)

  7. Example
  8. In this section the fields of a point charge with non-relativistic periodic motion are computed. (As a matter of definition, in the case of non-relativistic periodic motion, the maximum speed of the charge, vmax is <<c.) The charge moves along the x-axis, and the fields are computed at a few points on the y-axis. Of particular interest is the Poynting vector, S:

    (4_1)

    1. An Oscillating Point Charge
    2. Let point charge q=1 coul have the motion

      (4.1_1)

      We shall assume that A=1 meter and w=.01c/A. The wavelength of any radiation emitted by this particle is defined to be

      λ=2πc/ω (4.1_2)

      Program 1 (Section 10) computes Sy at various points on the y axis. Fig. 4.1_1 shows Sy(t) in the "near fields region," at y-axis points .05l above the x-axis. Note how energy appears to flux equally inward and outward periodically. Evidently, at this "near fields" distance, energy is pumped into the fields in one part of a cycle, and is pulled back out of the fields in the ensuing half cycle.

      Figure 4.1_1

      Sy, y=.05l

      In Fig. 4.1_2 Sy is plotted at a distance of .15l on the y-axis. Note how there is again energy flow out and in from the fields. But the flow in this case is predominantly outward. Also, the difference between the maxima and minima is considerably less than in Fig. 4.1_1.

      Figure 4.1_2

      Sy, y=.15l

      Finally, Fig. 4.1_3 depicts things at a full wavelength out on the y-axis. Practically all of the energy appears to pulse outward at this distance. And again, the difference between maxima and minima is less than in Fig. 4.1_2. The distances at which energy pulses outward is referred to as the "far fields" region. It extends out to infinity.

      Figure 4.1_3

      Sy, y=l

       

       

       

      1. Conservative Energy Flow
      2. Let us define energy that flows equally in and out of the fields as "conservative" energy flow. A typical case is illustrated in Fig. 4.1_1. Presumably such energy flow can be traced back to work expended by the external agent that drives the source charge. The energy is dubbed "conservative" because energy poured into the fields in one half-cycle is withdrawn from the fields in the ensuing half-cycle. For reasons to be discussed, we shall also refer to such conservative energy as "inertial" energy.

      3. Non-conservative Energy Flow

    Fig. 4.1_2 depicts things in a fringe zone, where the energy flow is a mixture of conservative and non-conservative. As in Fig. 4.1_1, the conservative energy flows outward and inward in equal amounts. But there is a net pulsing outward of energy into infinite space. In light of the fact that this energy is not recouped by the driving agent, it is dubbed "non-conservative." Here again, for reasons to be discussed, such non-conservative flow is also referred to as "radiant" energy flow.

    Fig. 4.1_3 depicts things in the "far fields" zone (which extends out to infinity). Here the energy flow is outward only, and is pulsed. Each cycle a certain amount of energy is lost from the system. Presumably the driving agent must do a net amount of work per cycle to account for this energy loss.

  9. Reaction Forces
  10. Let us assume that a point charge cannot itself be a source or sink of energy. Any positive work done by a driving agent presumably results in a flow of energy out into the fields. And any negative work results in a flow of energy in from the fields. In this document it is assumed that driving agents non-electromagnetically exert contact forces on point charges, and the same point charges may non-electromagnetically exert contact forces on the driving agents. The point charges can be thought of as "hooks" whereby the driving agents interact with the fields, and vice versa.

    1. Inertial Reaction Force
    2. Since Newton’s day it has been appreciated that mass has inertia. That is, in order to alter the velocity of a mass, an external force must be exerted. If the force has a component in the direction of the driven mass’s velocity, then the work equates to an increase in the mass’s kinetic energy. This energy is conservative in the sense that it can be recouped by the agent if the mass is decelerated.

      Within the present context, Newton’s 3rd law stipulates that whenever a driving agent exerts a force on a charge (and hence on its field), the charge/field exerts an equal and oppositely directed reaction force upon the agent. This reaction force, somewhat mysterious in Newton’s day, can be shown to be a self electric force when the mass is electromagnetic. That is, owing to the way an accelerated charge’s fields vary in time, electric fields are induced right at the charge. The resulting electric force ultimately acts upon the driving agent, who "responds" with an equal and oppositely directed force.

      The agent force does "conservative" work, and hence might be dubbed an inertial force. The self force would then be dubbed the inertial reaction force. In general whenever an inertial force is exerted on an electromagnetic mass, the fields will produce an equal and oppositely directed inertial reaction force right at the charge.

    3. Radiation Reaction Force

    As one might guess, part of the work done by a driving agent might also result in the creation of non-conservative (or radiant) energy. Here too, the fields induce an electric field right at the charge. The resulting electric force is equal and oppositely directed to the non-conservative part of the total agent force. It is accordingly usually referred to as the radiation reaction force

    Abraham and Lorentz were first to work out a (non-relativistic) formula for the radiation reaction force. According to them, in the case of motion along the x-axis,

    (5.2_1)

  11. Energy Conservation
  12. In this section we shall consider a periodically moving charge in order to emphasize the following result:

    The net work per cycle, expended by a driving agent to drive a charge periodically, equals the energy flow per cycle to/from the fields through a surrounding surface.

    Knowing the point charge general field solutions, it is a simple matter to compute the Poynting vector at points on and normal to a surrounding surface. Integration of the Poynting vector’s normal component, over the surface and over a cycle time, then provides the outward/inward energy flow per cycle time.

    The work done per cycle time by the inertial force (ma) is zero. The negative of the radiation reaction force, multiplied by the charge’s velocity (i.e., the power delivered by the agent in counteracting the radiation reaction force) can then be integrated over a cycle time to determine the work per cycle expended by the driving agent.

     

    1. Oscillating Charge
    2. Let a point charge have the motion

      (6.1_1)

      Then

      (6.1_2)

      The radiation reaction force is therefore

      (6.1_3)

      The part of the agent force that counteracts this is –FRadReact. The power expended by the driving agent is

      (6.1_4)

      And an integral of this power, over a cycle time, is the work done by the driving agent per cycle, in order to counteract the radiation reaction force:

      (6.1_5)

      As mentioned, the work done per cycle by the "conservative" inertial reaction force (and/or by the agent’s reaction to it) is zero. Hence the net work done per cycle expended to drive the charge is specified by Eq. 6.1_5.

      By way of an example, let us say that q=1 coul, A=1 meter, and w=.01c/A. Then W in Eq. 6.1_5 equals 18832 joules. Let an enclosing surface be a spherical surface of radius 5A, centered on the origin and with its poles lying on the x-axis. Program 2 (Section 11) computes the surface integral of the Poynting vector over 1 cycle time. The result of that exercise is 18833 joules.

    3. Circularly Orbiting Charge

    The exercise in the previous section can be repeated for a charge going in a circle in the xz plane at constant speed wR. In this case the acceleration points toward the circle’s center (at the origin) and hence is changing its direction (but not its amplitude) with time. da/dt thus has a constant amplitude and points opposite to the velocity vector at all times. Similarly for the radiation reaction force, whose amplitude is q2w3R/6pe0c3. The driving agent must counteract this radiation reaction force if the motion is to remain circular. The work done per cycle is therefore 2pR times the counteraction to the radiation reaction force:

    (6.2_2)

    Note that this is twice the work per cycle to drive the charge sinusoidally (Eq.. 6.1_5).

    Program 3 (Section 12) computes Ex and Ez for points on the y-axis at y=5A. Figs. 6.2_1 and 6.2_2 show the results. Evidently the vector E rotates parallel to the xz-plane. From the perspective of points on the y-axis (normal to the source charge’s orbital plane) the radiation is circularly polarized.

    Figure 6.2_1

    Ex

    Figure 6.2_2

    Ez

  13. Interactive Forces
  14. When two (or more) proximal point charges move in a known way, then each may experience not only its own self inertial and radiation reaction forces, but also an interactive force from the other charge. Perhaps the simplest example is two charges held at rest. Neither charge experiences a self force in this case, but each does experience an interactive (coulomb) force from the other charge. In order to maintain the state of rest, an external agent must counteract each of these interactive forces, which are equal and oppositely directed in the resting case.

    It might seem that whenever the distance between the two charges remains constant in time (e.g. when they share a common, periodic motion), then the net work expended per cycle to counteract the interactive forces would be zero. However, even at low frequencies and small amplitudes of motion this turns out not to be quite the case. If the charges move in prescribed ways, the field solutions (which are relativistically rigorous) and the Lorentz force law indicate that these forces are generally not simply Coulombic.

    1. Oscillating Charges
    2. Let us say that the x positions of two equal point charges (say q1=q2=1 coul) are

      (7.1_1)

      + A (7.1_2)

      y and z equal zero. d (the charge separation) and A (the amplitude of oscillation) are both 1 meter. The angular oscillation frequency is w=.01c/A. Each charge experiences an inertial reaction self force:

      (7.1_3)

      And each charge experiences a radiation reaction self force:

      (7.1_4)

      In order to maintain the specified motion, a driving agent must at the very least counteract the inertial and radiation reaction forces. (But the net work per cycle expended to counteract the inertial reaction force is zero.)

      Each charge also experiences an interactive Lorentz force in the other charge’s electric field. At any given moment, E can be computed using the general field solution algorithm. In order for the specified motion to be maintained, the driving agent must also counteract the interactive forces. The total work per cycle expended by the driving agent is the work done to counteract the radiation reaction forces, plus the work done to counteract the interactive forces.

      1. Work per Cycle to Drive The Two Charges.
      2. The work per cycle that must be expended in order to counteract the radiation reaction force is 37664 joules. (This is twice the work expended to counteract the radiation reaction force when only one, isolated charge is forced to oscillate.). Program 4 (Section 13) shows that the work per cycle expended to counteract the interactive forces is 37652 joules. Thus the total work per cycle that must be expended to counteract both the radiation reaction and interactive forces is 75317 joules.

        Program 2 is modified in Program 4 to compute the energy flux through a spherical surface enclosing the site of the two oscillating charges. Here care must be exercised to add the E and B field vectors before computing the Poynting vector. The energy flux per cycle is 75334 joules ... practically the same as the total work per cycle expended to counteract the radiation reaction and the interactive forces.

         

      3. Mass Coupling

      The electric fields responsible for self forces do not exist solely at the charges. They extend into the surrounding space and may affect other charges in the vicinity. For example, two tiny spherical shells of charge, with equal electromagnetic masses, interact with one another when they move (e.g. when they oscillate) in tandem. The motion-induced electric field right at one of the charges, responsible for the inertial reaction force, may also exist somewhat attenuated at the other charge’s position. Indeed if the charges are virtually on top of one another, then the inertial reaction field at each charge may be twice that of a single, isolated charge. This being the case, the inertial reaction force on the pair is four times that on one charge that is alone and isolated. In brief, unlike mechanical masses, electromagnetic masses do not add. As already indicated, the electromagnetic mass of a charge is proportional to the square of q. Doubling q increases the electromagnetic mass by a factor of 4!

    3. Circling Charges

    Let us suppose that, viewed from the positive z-axis, two equal charges travel in a CCW circle of radius A in the xy-plane. The (centrifugal) inertial reaction forces (and the driving agent’s centripetal counteractions) point normal to the velocity vectors at all time, and thus do no work. (Any magnetic field will have only a z-component. Any magnetic force on a subject charge will also point normal to the charge’s velocity.)

    The radiation reaction forces point opposite to the velocities. The work per cycle, expended to counteract each radiation reaction force, is q2w3A2/(3e0c3). The work per cycle expended to counteract both radiation reaction forces is twice this amount.

    Each charge also experiences an interactive force attributable to the other charge. One of the components of this force acts normal to the subject charge’s velocity, and thus only serves to modify the inertial reaction force and/or the magnetic force. Again, the work per cycle expended by the driving agent to counteract this component of the interactive force is zero. The other component of the interactive force acts tangentially to the subject charge’s circular orbit. It therefore modifies the radiation reaction force.

    Program 5 (Section 14) computes the tangential components of the interactive forces and the work per orbit expended to counteract these components. Perhaps surprisingly, the interactive force acting upon a charge points in the same direction as the charge’s velocity. This decreases the work per cycle that must be expended to counteract the radiation reaction force. Indeed for values used in the program, the work per orbit expended to counteract the radiation reaction forces is 75329 joules, whereas the work to counteract the tangential components of the interactive forces if -75308 joules. The net work per cycle, required to drive both charges in tandem around the circular orbit, is thus only 21 joules!

    It might be expected that by adding more charges, evenly distributed around a circular orbit, the attenuating effect of the interactive forces would be even more pronounced. Indeed if an infinite number of infinitesimal point charges are evenly distributed around a circular orbit, then we have what amounts to a circular current loop. And it takes no power at all to drive a continuous, circular line charge at constant speed around a circle. Of course a corollary to this result is that a constant, circular current emits no radiant energy.

  15. Relativistic Motions
  16. Although the general field solutions are relativistically rigorous, the formulas for the inertial and radiation reaction forces are not. Thus the formulas for the agent’s counteraction to these forces must be modified.

     

    1. The Inertial Reaction Force
    2. The correction to the inertial reaction force amounts to a modification of the electromagnetic mass:

      (8.1_1)

      In words, the mass of a particle increases as its speed relative to one’s chosen inertial frame increases. For values of v<<c, the mass differs little from m0 (the rest mass). But as v approaches c, the inertial mass of the particle approaches infinity. Eq. 8.1_1 applies to both "mechanical" mass (e.g. the mass of an uncharged particle) and to electromagnetic mass.

      It turns out that Newton’s 2nd law, in the form F=d(mv)/dt, is still valid provided the dependence of mass upon speed is taken into account. For example, in 1 dimension,

      (8.1_2)

      Or, in view of eq. 8.1_1,

      (8.1_3)

      1. An Oscillating Charge

      Let an oscillating charge’s motion be described by

      (8.1.1_1)

      A (8.1.1_2)

      Program 6 (Section 15) computes Fx(t) and the work per cycle expended to drive the charge. Fig. 8.1.1_1 plots Fx for the non-relativistic case where w=.01c/A. As expected, the agent force is essentially sinusoidal. The computed work per cycle is 7e-3 joules (practically zero).

      Figure 8.1.1_1

      Fx, Non-Relativistic Oscillator

      Fig. 8.1.1_2 plots Fx for the relativistic case where wA=.95c. The computed work per cycle is -100 joules (again, practically zero considering the extrema).

      Figure 8.1.1_2

      Fx, Relativistic Oscillator

      Since the computed work per cycle is zero in both the non-relativistic and relativistic cases, it is clear that Fx is conservative in both cases. The amount of energy that fluxes out to and in from the fields is relatively enormous in the relativistic case. Yet Fx is symmetric in time, and .

       

    3. The Radiation Reaction Force

    The relativistic version of the Abraham-Lorentz formula for the radiation reaction force is, in one dimension,

    ) (8.2_1)

    In order to maintain the periodic motion, the agent must again counteract FRR. The energy per cycle, expended in doing do, is invariably >0.

    Program 7 (Section 16) computes –FRR(t)vx(t) and the work per cycle expended by this part of the total agent force. Fig. 8.2_1 plots Fxvx for the non-relativistic case where w=.01c/A. The computed work per cycle is 18333 joules. As a check, the energy flux per cycle time through an enclosing surface is computed and is found to be 18333 joules.

    Figure 8.2_1

    Power Expended to Drive Non-relativistic Oscillator.

    Fig. 8.2_2 plots Fxvx for the non-relativistic case where wA=.95c. The computed work per cycle is 1.7E11 joules. Noteworthy in the figure are the 4 short blips of negative radiation!

    Figure 8.2_2

    Power Expended to Drive Relativistic Oscillator

    We can again test the essential correctness of Eq. 8.2_1 by computing the energy flux per cycle through a surrounding surface. That flux is computed to be 1.7E11 joules.

  17. Interactions Over a Range of Distances
  18. In a previous section we computed WInt, the work per cycle expended to drive two oscillating charges that are held at a constant distance d apart. The line of oscillation coincided with the line separating the charges. Specifically, for a separation of d=A=1 meter and an angular frequency of w=.01c/A, it was found that 37652 joules had to be expended each cycle in order for the driving agent to counteract the interactive forces. (Work also had to be expended to counteract the radiation reaction forces.) It is perhaps worth re-emphasizing that this is not the 0 joules which might be expected if the charges were at rest (i.e., if w equaled zero). In the resting case the interactive forces are Coulombic and, among other things, are equal and oppositely directed. In the oscillating case ... even in non-relativistic cases ... the interactive forces are not equal and oppositely directed, and work must be done each cycle in order to counteract those forces and maintain the sinusoidal motion.

    In this section we investigate the dependence of the work per cycle on d, all other things being kept constant. Given a range of d from, say, .5l to 3l (where l=2pc/w), how does WInt change with increasing d? Program 8 (Section 17) provides an answer. Fig. 9_1 plots WInt vs. d when the line of oscillation coincides with the line of separation. Note that WInt does at first drop off with increasing d, but then rises again to a secondary maximum at approximately d=1.5l, etc.

    Figure 9_1

    Work per Cycle to Counteract Int. Forces, Varying Separations

    It is also interesting to determine how WInt depends on d when the line of motion is perpendicular to the line separating the charges. Fig. 9_2 plots the total work expended to counteract the interactive forces vs. d. (Here again, a set amount of work is also expended each cycle to counteract the radiation reaction forces.)

    Figure 9_2

    Work per Cycle to Counteract Int. Forces, Varying Perpendicular Separations

    It is clear in Figs. 9_1 and 9_2 that the work per cycle, expended to counteract interactive forces, oscillates with increasing d. These results may be of interest in antenna design. Also of interest is the fact that, at selected values of d, WInt may actually be negative. A ramification is that, at such values of d, it will take less energy per cycle to drive the two proximal charges than twice what it would take to drive either charge alone. That is, at such selected values of d the interactive forces point opposite to the radiation reaction forces. Hence the interactive forces assist the agent in counteracting the radiation reaction forces. A simple trap in the program, however, indicates that |Wint| is always less than WRR.

  19. Program 1
  20. Private Sub cmdProgram1_Click()

    'PROGRAM 1

    '*******************

    'Compute Sy at points on the y axis, given a charge

    'that oscillates along the x axis.

    '*******************

    Const c As Double = 300000000# 'speed of light

    Const eps0 As Double = 0.00000000000885 'permittivity constant

    Const pi As Double = 3.141592654

    Const q As Double = 1 'charge equals 1 coulomb

    Const A As Double = 1 'amplitude of oscillation equals 1 meter

    Const omega As Double = 0.01 * c / A 'angular frequency (non-relativistic)

    Const steps As Long = 500 'number of iterations

    Const freq As Double = omega / (2 * pi)

    Const tau As Double = 1 / freq

    Const deltat As Double = tau / steps 'time interval between epochs

    Const lambda As Double = c / freq

    'Const Py As Double = 0.05 * lambda

    'Const Py As Double = 0.15 * lambda

    Const Py As Double = lambda

    Dim i As Long 'loop index

    Dim t(steps) As Double 'present time

    Dim Ex(steps), Ey(steps) As Double 'electric field components

    Dim Bz(steps) As Double 'magnetic field component

    Dim Sy(steps) As Double 'Poynting vector

    Dim dtmin, dtmax, dt As Double

    Dim ux, uy As Double

    Dim drx, dry, dr As Double

    Dim tr As Double 'retarded time

    Dim xr As Double

    Dim vr As Double

    Dim ar As Double

    For i = 0 To steps - 1

    Debug.Print i

    t(i) = i * deltat

    dtmin = 0

    dtmax = 5 * (A + Py)

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    xr = A * Sin(omega * tr)

    drx = -xr

    dry = Py

    dr = Sqr(drx ^ 2 + dry ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    vr = omega * A * Cos(omega * tr)

    ar = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - vr

    uy = c * dry / dr

    Ex(i) = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (ux * (c ^ 2 - vr ^ 2) + dry * -uy * ar)

    Ey(i) = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (uy * (c ^ 2 - vr ^ 2) - drx * (-uy * ar))

    Bz(i) = 1 / (c * dr) * (drx * Ey(i) - dry * Ex(i))

    Sy(i) = eps0 * c ^ 2 * -Ex(i) * Bz(i)

    Next i

    Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

    For i = 0 To steps - 1

    Write #1, t(i), Sy(i)

    Next i

    Close

    MsgBox ("Ready for plotting")

    Stop

    End Sub

  21. Program 2
  22. Private Sub cmdProgram2_Click()

    'PROGRAM 2

    '*******************

    'Compute the energy flux per cycle time through a surface containing

    'an oscillating charge.

    '*******************

    Const c As Double = 300000000# 'speed of light

    Const eps0 As Double = 0.00000000000885 'permittivity constant

    Const pi As Double = 3.141592654

    Const q As Double = 1 'charge equals 1 coulomb

    Const A As Double = 1 'amplitude of oscillation equals 1 meter

    Const R As Double = 5 'radius of surrounding surface

    Const omega As Double = 0.01 * c / A 'angular frequency (non-relativistic)

    Const steps As Long = 500 'number of iterations

    Const dtheta As Double = pi / steps

    Const freq As Double = omega / (2 * pi)

    Const tau As Double = 1 / freq

    Const deltat As Double = tau / steps 'time interval between epochs

    Dim i As Long 'loop index

    Dim j As Long

    Dim theta As Double

    Dim dArea As Double

    Dim Px, Py As Double

    Dim t(steps) As Double 'present time

    Dim Ex, Ey As Double 'electric field components

    Dim Bz As Double 'magnetic field component

    Dim Sx, Sy As Double 'Poynting vector

    Dim Snormal As Double

    Dim nx, ny As Double

    Dim dtmin, dtmax, dt As Double

    Dim ux, uy As Double

    Dim drx, dry, dr As Double

    Dim tr As Double 'retarded time

    Dim xr As Double

    Dim vr As Double

    Dim ar As Double

    Dim WorkRR As Double

    Dim EFluxRate(steps) As Double

    Dim Eflux(steps) As Double

    Dim TotalFlux As Double

    WorkRR = q ^ 2 * omega ^ 3 * A ^ 2 / (6 * eps0 * c ^ 3)

    MsgBox ("Work done to counteract Rad React force = " & WorkRR)

    TotalFlux = 0

    For i = 0 To steps - 1

    Debug.Print i

    t(i) = i * deltat

    Eflux(i) = 0

    For j = 0 To steps - 1

    EFluxRate(j) = 0

    Next j

    For j = 0 To steps - 1

    theta = (j + 0.5) * dtheta

    nx = Cos(theta)

    ny = Sin(theta)

    Px = R * Cos(theta)

    Py = R * Sin(theta)

    dArea = 2 * pi * Py * R * dtheta

    dtmin = 0

    dtmax = 5 * (R + A)

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    xr = A * Sin(omega * tr)

    drx = Px - xr

    dry = Py

    dr = Sqr(drx ^ 2 + dry ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    vr = omega * A * Cos(omega * tr)

    ar = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - vr

    uy = c * dry / dr

    Ex = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (ux * (c ^ 2 - vr ^ 2) + dry * (-uy * ar))

    Ey = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (uy * (c ^ 2 - vr ^ 2) - drx * (-uy * ar))

    Bz = 1 / (c * dr) * (drx * Ey - dry * Ex)

    Sx = eps0 * c ^ 2 * (Ey * Bz)

    Sy = eps0 * c ^ 2 * (-Ex * Bz)

    Snormal = Sx * nx + Sy * ny

    EFluxRate(j) = EFluxRate(j) + Snormal * dArea

    Next j

    For j = 0 To steps - 1

    Eflux(i) = Eflux(i) + EFluxRate(j) * deltat

    Next j

    Next i

    For i = 0 To steps - 1

    TotalFlux = TotalFlux + Eflux(i)

    Next i

    MsgBox ("Total Flux = " & TotalFlux)

    Stop

    End Sub

  23. Program 3
  24. Private Sub cmdProgram3_Click()

    'PROGRAM 3

    '*******************

    'Compute Ex and Ez at points on the y axis, given a charge

    'that goes in a circle in the xz plane.

    '*******************

    Const c As Double = 300000000# 'speed of light

    Const eps0 As Double = 0.00000000000885 'permittivity constant

    Const pi As Double = 3.141592654

    Const q As Double = 1 'charge equals 1 coulomb

    Const A As Double = 1 'amplitude of oscillation equals 1 meter

    Const omega As Double = 0.01 * c / A 'angular frequency (non-relativistic)

    Const steps As Long = 500 'number of iterations

    Const freq As Double = omega / (2 * pi)

    Const tau As Double = 1 / freq

    Const deltat As Double = tau / steps 'time interval between epochs

    Const Py As Double = 5 * A

    Dim i As Long 'loop index

    Dim t(steps) As Double 'present time

    Dim x(steps), z(steps) As Double

    Dim Px, Pz As Double

    Dim Ex(steps) As Double 'electric field components

    Dim Ez(steps) As Double 'electric field component

    Dim dtmin, dtmax, dt As Double

    Dim ux, uy, uz As Double

    Dim drx, dry, drz, dr As Double

    Dim tr As Double 'retarded time

    Dim xr, yr, zr As Double

    Dim vxr, vyr, vzr, vr As Double

    Dim axr, ayr, azr As Double

    For i = 0 To steps - 1

    ' t(i) = i * deltat

    ' x(i) = A * Cos(omega * t(i))

    ' z(i) = A * Sin(omega * t(i))

    'Next i

    'For i = 0 To steps - 1

    Debug.Print i

    t(i) = i * deltat

    dtmin = 0

    dtmax = 5 * (A + Py)

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    xr = A * Cos(omega * tr)

    yr = 0

    'yr = Py

    zr = A * Sin(omega * tr)

    drx = -xr

    dry = Py

    drz = -zr

    dr = Sqr(drx ^ 2 + dry ^ 2 + drz ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    vxr = -omega * A * Sin(omega * tr)

    vyr = 0

    vzr = omega * A * Cos(omega * tr)

    vr = Sqr(vxr ^ 2 + vyr ^ 2 + vzr ^ 2)

    axr = -(omega ^ 2) * A * Cos(omega * tr)

    'axr = -(omega ^ 2) * A * Sin(omega * tr)

    ayr = 0

    azr = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - vxr

    uy = c * dry / dr - vyr

    uz = c * drz / dr - vzr

    Ex(i) = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy + drz * uz) ^ 3 * (ux * (c ^ 2 - vr ^ 2) + dry * (ux * ayr - uy * axr) - drz * (uz * axr - ux * axr))

    Ez(i) = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy + drz * uz) ^ 3 * (uz * (c ^ 2 - vr ^ 2) + drx * (uz * axr - ux * azr) - dry * (uy * azr - uz * ayr))

    Next i

    Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

    For i = 0 To steps - 1

    Write #1, t(i), Ex(i), Ez(i)

    Next i

    Close

    MsgBox ("Ready for plotting")

    Stop

    End Sub

  25. Program 4
  26. Private Sub cmdProgram4_Click()

    'PROGRAM 4

    '*****************

    'Compute the work per cycle that must be done

    'to drive 2 co-oscillating charges along the x-axis.

    'Then compute the energy flux through an enclosing surface.

    '*****************

    'First compute the work per cycle expended to counteract

    'the radiation reaction and the interactive forces.

    Const c As Double = 300000000#

    Const eps0 As Double = 0.00000000000885

    Const pi As Double = 3.141592654

    Const steps As Long = 500

    Const q As Double = 1

    Const A As Double = 1

    Const omega As Double = 0.01 * c / A

    Const d As Double = 1

    Const freq As Double = omega / (2 * pi)

    Const period As Double = 1 / freq

    Const deltat As Double = period / steps

    Dim i As Long

    Dim t(steps) As Double

    Dim x1, x2 As Double

    Dim v2 As Double

    Dim dtmin, dtmax As Double

    Dim dt As Double

    Dim tr As Double

    Dim x1r As Double

    Dim drx, dr As Double

    Dim v1r, a1r As Double

    Dim ux As Double

    Dim Ex As Double

    Dim Fx As Double

    Dim P(steps) As Double

    Dim Work, WorkRR, WorkTotal As Double

    'Compute and display the work per cycle expended to

    'counteract the radiation reaction forces.

    WorkRR = q ^ 2 * omega ^ 3 * A ^ 2 / (6 * eps0 * c ^ 3)

    WorkRR = 2 * WorkRR

    MsgBox ("WorkRR = " & WorkRR)

    'Now compute the work per cycle expended to counteract

    'the interactive forces.

    For i = 0 To steps - 1

    t(i) = i * deltat

    x1 = -d / 2 + A * Sin(omega * t(i))

    x2 = d / 2 + A * Sin(omega * t(i))

    v2 = omega * A * Cos(omega * t(i))

    dtmin = 0

    dtmax = 5 * (A + d) / c

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    x1r = -d / 2 + A * Sin(omega * tr)

    drx = x2 - x1r

    dr = Abs(drx)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    v1r = omega * A * Cos(omega * tr)

    a1r = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - v1r

    Ex = q / (4 * pi * eps0) * dr / (drx * ux) ^ 3 * (ux * (c ^ 2 - v1r ^ 2))

    Fx = -q * Ex

    P(i) = Fx * v2

    Next i

    Work = 0

    For i = 0 To steps - 1

    Work = Work + P(i) * deltat

    Next i

    Work = 2 * Work

    MsgBox ("Interactive Work = " & Work)

    WorkTotal = WorkRR + Work

    MsgBox ("Total Work per Cycle = " & WorkTotal)

    'Now compute the energy flux per cycletime through

    'an enclosing surface.

    Dim TotalFlux As Double

    Dim Eflux(steps) As Double

    Dim j As Long

    Dim theta As Double

    Const dtheta As Double = pi / steps

    Dim EFluxRate(steps) As Double

    Dim nx, ny As Double

    Dim dArea As Double

    Const R As Double = 5 * A

    Dim Px, Py As Double

    Dim dry, uy As Double

    Dim Ey, Bz As Double

    Dim x2r, v2r, a2r As Double

    Dim Sx, Sy, Snormal As Double

    Dim Ex2, Ey2, Bz2 As Double

    For i = 0 To steps - 1

    Debug.Print i

    t(i) = i * deltat

    Eflux(i) = 0

    EFluxRate(i) = 0

    For j = 0 To steps - 1

    theta = (j + 0.5) * dtheta

    nx = Cos(theta)

    ny = Sin(theta)

    Px = R * Cos(theta)

    Py = R * Sin(theta)

    dArea = 2 * pi * Py * R * dtheta

    dtmin = 0

    dtmax = 5 * (R + A)

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    x1r = -d / 2 + A * Sin(omega * tr)

    drx = Px - x1r

    dry = Py

    dr = Sqr(drx ^ 2 + dry ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    v1r = omega * A * Cos(omega * tr)

    a1r = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - v1r

    uy = c * dry / dr

    Ex = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (ux * (c ^ 2 - v1r ^ 2) + dry * (-uy * a1r))

    Ey = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (uy * (c ^ 2 - v1r ^ 2) - drx * (-uy * a1r))

    Bz = 1 / (c * dr) * (drx * Ey - dry * Ex)

    dtmin = 0

    dtmax = 5 * (R + A)

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    x2r = d / 2 + A * Sin(omega * tr)

    drx = Px - x2r

    dry = Py

    dr = Sqr(drx ^ 2 + dry ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    v2r = omega * A * Cos(omega * tr)

    a2r = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - v2r

    uy = c * dry / dr

    Ex2 = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (ux * (c ^ 2 - v2r ^ 2) + dry * (-uy * a2r))

    Ey2 = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (uy * (c ^ 2 - v2r ^ 2) - drx * (-uy * a2r))

    Bz2 = 1 / (c * dr) * (drx * Ey2 - dry * Ex2)

    Ex = Ex + Ex2

    Ey = Ey + Ey2

    Bz = Bz + Bz2

    Sx = eps0 * c ^ 2 * (Ey * Bz)

    Sy = eps0 * c ^ 2 * (-Ex * Bz)

    Snormal = Sx * nx + Sy * ny

    EFluxRate(j) = Snormal * dArea

    Next j

    For j = 0 To steps - 1

    Eflux(i) = Eflux(i) + EFluxRate(j) * deltat

    Next j

    Next i

    TotalFlux = 0

    For i = 0 To steps - 1

    TotalFlux = TotalFlux + Eflux(i)

    Next i

    MsgBox ("Total Flux = " & TotalFlux)

    Stop

    End Sub

  27. Program 5
  28. Private Sub cmdProgram5_Click()

    'PROGRAM 5

    '*****************

    'Compute the work per cycle that must be done

    'to drive 2 diametrically opposed charges around a circle.

    '*****************

    Const c As Double = 300000000#

    Const eps0 As Double = 0.00000000000885

    Const pi As Double = 3.141592654

    Const q As Double = 1

    Const A As Double = 1

    Const omega As Double = 0.01 * c / A

    Dim x2, y2 As Double 'right (subject)charge

    Dim v2 As Double

    Dim dtmin, dtmax As Double

    Dim dt As Double

    Dim t, tr As Double

    Dim x1r, y1r As Double 'left (source) charge

    Dim drx, dry, dr As Double

    Dim v1xr, v1yr, a1xr, a1yr, v1r As Double

    Dim ux, uy As Double

    Dim Ey As Double

    Dim Fy As Double

    Dim Work, WorkRR, WorkInt, WorkTotal As Double

    'Compute the work per cycle to counteract the radiation

    'reaction forces

    WorkRR = q ^ 2 * omega ^ 3 * A ^ 2 / (3 * eps0 * c ^ 3)

    WorkRR = 2 * WorkRR

    MsgBox ("WorkRR = " & WorkRR)

    'Compute Ey at q2

    t = 0

    x2 = A

    y2 = 0

    dtmin = 0

    dtmax = 5 * A / c

    Do

    dt = (dtmin + dtmax) / 2

    tr = t - dt

    x1r = -A * Cos(omega * tr)

    y1r = -A * Sin(omega * tr)

    drx = x2 - x1r

    dry = y2 - y1r

    dr = Sqr(drx ^ 2 + dry ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    v1xr = omega * A * Sin(omega * tr)

    v1yr = -omega * A * Cos(omega * tr)

    v1r = Sqr(v1xr ^ 2 + v1yr ^ 2)

    a1xr = omega ^ 2 * A * Cos(omega * tr)

    a1yr = omega ^ 2 * A * Sin(omega * tr)

    ux = c * drx / dr - v1xr

    uy = c * dry / dr - v1yr

    Ey = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (uy * (c ^ 2 - v1r ^ 2) - drx * (ux * a1yr - uy * a1xr))

    Fy = -q * Ey

    WorkInt = Fy * 2 * pi * A

    WorkInt = 2 * WorkInt

    MsgBox ("WorkInt = " & WorkInt)

    WorkTotal = WorkRR + WorkInt

    MsgBox ("WorkTotal = " & WorkTotal)

    Stop

    End Sub

  29. Program 6
  30. Private Sub cmdProgram6_Click()

    'PROGRAM 6

    '******************

    'Compute the inertial force and work per cycle required to drive

    'a charged particle sinusoidally, at non-relativistic and

    'relativistic values of wA.

    '******************

    Const c As Double = 300000000#

    Const eps0 As Double = 0.00000000000885

    Const pi As Double = 3.141592654

    Const steps As Long = 1500

    Const q As Double = 1

    Const A As Double = 1

    'Const omega As Double = 0.01 * c / A

    Const omega As Double = 0.95 * c / A

    Const freq As Double = omega / (2 * pi)

    Const tau As Double = 1 / freq

    Const deltat As Double = tau / steps

    Const m0 As Double = 1

    Dim WorkPerCycle As Double

    Dim Fx(steps) As Double

    Dim t(steps) As Double

    Dim index As Long

    Dim vx As Double

    Dim ax As Double

    Dim gamma As Double

    WorkPerCycle = 0

    For index = 0 To steps - 1

    t(index) = (index + 0.5) * deltat

    vx = omega * A * Cos(omega * t(index))

    ax = -(omega ^ 2) * A * Sin(omega * t(index))

    gamma = 1 / Sqr(1 - vx ^ 2 / c ^ 2)

    Fx(index) = gamma ^ 3 * m0 * ax

    WorkPerCycle = WorkPerCycle + Fx(index) * vx * deltat

    Next index

    MsgBox ("Work per cycle = " & WorkPerCycle)

    Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

    For index = 0 To steps - 1

    Write #1, t(index), Fx(index)

    Next index

    Close

    MsgBox ("Ready for plotting")

    Stop

    End Sub

  31. Program 7
  32. Private Sub cmdProgram7_Click()

    'PROGRAM 7

    '******************

    'Compute the rad react force and work per cycle required to drive

    'a charged particle sinusoidally, at non-relativistic and

    'relativistic values of wA.

    '******************

    Const c As Double = 300000000#

    Const eps0 As Double = 0.00000000000885

    Const pi As Double = 3.141592654

    Const steps As Long = 500

    Const q As Double = 1

    Const A As Double = 1

    'Const omega As Double = 0.01 * c / A

    Const omega As Double = 0.95 * c / A

    Const freq As Double = omega / (2 * pi)

    Const tau As Double = 1 / freq

    Const deltat As Double = tau / steps

    Const lambda As Double = c / freq

    Const R As Double = 5 * A

    Dim WorkPerCycle As Double

    Dim FRR(steps) As Double

    Dim t(steps) As Double

    Dim i As Long

    Dim vx, ax As Double

    Dim daxdt As Double

    Dim gamma As Double

    WorkPerCycle = 0

    For i = 0 To steps - 1

    t(i) = (i + 0.5) * deltat

    vx = omega * A * Cos(omega * t(i))

    ax = -omega ^ 2 * A * Sin(omega * t(i))

    daxdt = -(omega ^ 3) * A * Cos(omega * t(i))

    gamma = 1 / Sqr(1 - vx ^ 2 / c ^ 2)

    FRR(i) = gamma ^ 4 * q ^ 2 / (6 * pi * eps0 * c ^ 3) * (daxdt + 3 * gamma ^ 2 * vx * ax ^ 2 / c ^ 2)

    WorkPerCycle = WorkPerCycle - FRR(i) * vx * deltat

    Next i

    MsgBox ("Work per cycle = " & WorkPerCycle)

    Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

    For i = 0 To steps - 1

    vx = omega * A * Cos(omega * t(i))

    Write #1, t(i), -FRR(i) * vx

    Next i

    Close

    MsgBox ("Ready for plotting")

    'As a check, compute the energy flux per cycle time through

    'a surrounding surface.

    Const dtheta As Double = pi / steps

    Dim j As Long

    Dim theta As Double

    Dim dArea As Double

    Dim Px, Py As Double

    Dim Ex, Ey As Double 'electric field components

    Dim Bz As Double 'magnetic field component

    Dim Sx, Sy As Double 'Poynting vector

    Dim Snormal As Double

    Dim nx, ny As Double

    Dim dtmin, dtmax, dt As Double

    Dim ux, uy As Double

    Dim drx, dry, dr As Double

    Dim tr As Double 'retarded time

    Dim xr As Double

    Dim vr As Double

    Dim ar As Double

    Dim WorkRR As Double

    Dim EFluxRate(steps) As Double

    Dim Eflux(steps) As Double

    Dim TotalFlux As Double

    TotalFlux = 0

    For i = 0 To steps - 1

    Debug.Print i

    Eflux(i) = 0

    For j = 0 To steps - 1

    EFluxRate(j) = 0

    Next j

    For j = 0 To steps - 1

    theta = (j + 0.5) * dtheta

    nx = Cos(theta)

    ny = Sin(theta)

    Px = R * Cos(theta)

    Py = R * Sin(theta)

    dArea = 2 * pi * Py * R * dtheta

    dtmin = 0

    dtmax = 5 * (R + A)

    Do

    dt = (dtmin + dtmax) / 2

    tr = t(i) - dt

    xr = A * Sin(omega * tr)

    drx = Px - xr

    dry = Py

    dr = Sqr(drx ^ 2 + dry ^ 2)

    If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

    If c * dt - dr > 0 Then

    dtmax = dt

    Else

    dtmin = dt

    End If

    Loop

    vr = omega * A * Cos(omega * tr)

    ar = -(omega ^ 2) * A * Sin(omega * tr)

    ux = c * drx / dr - vr

    uy = c * dry / dr

    Ex = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (ux * (c ^ 2 - vr ^ 2) + dry * (-uy * ar))

    Ey = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (uy * (c ^ 2 - vr ^ 2) - drx * (-uy * ar))

    Bz = 1 / (c * dr) * (drx * Ey - dry * Ex)

    Sx = eps0 * c ^ 2 * (Ey * Bz)

    Sy = eps0 * c ^ 2 * (-Ex * Bz)

    Snormal = Sx * nx + Sy * ny

    EFluxRate(j) = EFluxRate(j) + Snormal * dArea

    Next j

    For j = 0 To steps - 1

    Eflux(i) = Eflux(i) + EFluxRate(j) * deltat

    Next j

    Next i

    Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

    For i = 0 To steps - 1

    Write #1, t(i), Eflux(i)

    Next i

    Close

    MsgBox ("Ready for plotting")

    For i = 0 To steps - 1

    TotalFlux = TotalFlux + Eflux(i)

    Next i

    MsgBox ("Total Flux = " & TotalFlux)

    Stop

    End Sub

  33. Program 8

Private Sub cmdProgram8_Click()

'PROGRAM 8

'*******************

'Compute the work per cycle expended to drive 2 co-oscillating

'charges at each separation in a range of separations.

'*******************

Const c As Double = 300000000#

Const eps0 As Double = 0.00000000000885

Const pi As Double = 3.141592654

Const steps As Long = 500

Const q As Double = 1

Const A As Double = 1

Const omega As Double = 0.01 * c / A

Const freq As Double = omega / (2 * pi)

Const period As Double = 1 / freq

Const deltat As Double = period / steps

Const lambda As Double = c / freq

Const dd As Double = (3 * lambda - (0.5 * lambda)) / 100

Dim d(100) As Double

Dim i, j As Long

Dim t(steps) As Double

Dim x1, x2 As Double

Dim v2 As Double

Dim dtmin, dtmax As Double

Dim dt As Double

Dim tr As Double

Dim x1r As Double

Dim drx, dr As Double

Dim v1r, a1r As Double

Dim ux As Double

Dim Ex As Double

Dim Fx As Double

Dim P(steps) As Double

Dim Work(100) As Double

For j = 0 To 99

d(j) = 0.5 * lambda + j * dd

Next j

'Now compute the works per cycle expended to counteract

'the interactive forces.

For j = 0 To 99

Debug.Print j

For i = 0 To steps - 1

t(i) = i * deltat

x1 = -d(j) / 2 + A * Sin(omega * t(i))

x2 = d(j) / 2 + A * Sin(omega * t(i))

v2 = omega * A * Cos(omega * t(i))

dtmin = 0

dtmax = 5 * (A + d(j)) / c

Do

dt = (dtmin + dtmax) / 2

tr = t(i) - dt

x1r = -d(j) / 2 + A * Sin(omega * tr)

drx = x2 - x1r

dr = Abs(drx)

If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

If c * dt - dr > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

v1r = omega * A * Cos(omega * tr)

a1r = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * drx / dr - v1r

Ex = q / (4 * pi * eps0) * dr / (drx * ux) ^ 3 * (ux * (c ^ 2 - v1r ^ 2))

Fx = -q * Ex

P(i) = Fx * v2

Next i

Work(j) = 0

For i = 0 To steps - 1

Work(j) = Work(j) + P(i) * deltat

Next i

Work(j) = 2 * Work(j)

Next j

Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

For j = 0 To 99

Write #1, d(j) / lambda, Work(j)

Next j

Close

MsgBox ("Ready for plotting")

Dim y1, y2 As Double

Dim v As Double

Dim y1r As Double

Dim dry As Double

Dim uy As Double

Dim WorkRR As Double

For j = 0 To 99

d(j) = 0.5 * lambda + j * dd

Next j

For j = 0 To 99

Debug.Print j

For i = 0 To steps - 1

t(i) = i * deltat

x1 = A * Sin(omega * t(i))

y1 = -d(j) / 2

x2 = A * Sin(omega * t(i))

y2 = d(j) / 2

v = omega * A * Cos(omega * t(i))

dtmin = 0

dtmax = 5 * (A + d(j)) / c

Do

dt = (dtmin + dtmax) / 2

tr = t(i) - dt

x1r = A * Sin(omega * tr)

y1r = y1

drx = x2 - x1r

dry = y2 - y1r

dr = Sqr(drx ^ 2 + dry ^ 2)

If Abs(c * dt - dr) < 2 ^ (-30) Then Exit Do

If c * dt - dr > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

v1r = omega * A * Cos(omega * tr)

a1r = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * drx / dr - v1r

uy = c * dry / dr

Ex = q / (4 * pi * eps0) * dr / (drx * ux + dry * uy) ^ 3 * (ux * (c ^ 2 - v1r ^ 2) + dry * -uy * a1r)

Fx = -q * Ex

P(i) = Fx * v

Next i

Work(j) = 0

For i = 0 To steps - 1

Work(j) = Work(j) + P(i) * deltat

Next i

Work(j) = 2 * Work(j)

Next j

WorkRR = q ^ 2 * omega ^ 3 * A ^ 2 / (6 * eps0 * c ^ 3)

WorkRR = 2 * WorkRR

For j = 0 To 99

If Work(j) < 0 Then

If Abs(Work(j)) > WorkRR Then

MsgBox ("Work = " & Work(j) & " and WorkRR = " & WorkRR)

End If

End If

Next j

Open "c:\WINMCADC\Physics\PosSy.prn" For Output As #1

For j = 0 To 99

Write #1, d(j) / lambda, Work(j)

Next j

Close

MsgBox ("Ready for plotting")

Stop

End Sub