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.
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.
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)
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.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.
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.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.
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)
The momentum density in the electromagnetic field is
g
=Integration produces the result that the total momentum in the electromagnetic field of a spherical shell, moving at constant, non-relativistic velocity v, is

In light of this, the electromagnetic (rest) mass of such a shell is defined to be
![]()
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.
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.
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
.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.)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

![]()
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)
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)
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 .05
l 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=.05
lIn Fig. 4.1_2 Sy is plotted at a distance of .15
l 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=.15
lFinally, 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
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.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.
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.
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.
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)
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.
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.
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
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.
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.
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.
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!
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 q2
w3A2/(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.
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.
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)
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 
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.
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, .5
l 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.
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
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
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
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
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
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
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
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