Home

An Example of Field Energy Absorption

Abstract. The energy flux through surfaces containing unequal oscillating charges is shown to equal the work done to drive each charge.

Let us imagine that positive point charges q1 and q2 are separated by a distance L equal to one wavelength. The charges oscillate in phase on the x-axis. Their motions are:

,

.

Assuming some agent holds the separation L constant, and assuming this agent counteracts the radiation reaction forces and the interactive forces, an amount of work W must be done each cycle. If q1 = q2 then W = W1+W2 = 2W1. And if we compute the energy flux per cycle time through a spherical surface with both charges oscillating inside, then that energy flux equals W:

.

It is a simple matter to make q1 unequal to q2, say q1 = .0001 q2. And in this case it turns out that W1 is actually negative, whereas W2 is still positive. Furthermore, if we construct a surface enclosing only q1 and compute the energy flux per cycle time through that surface, then that energy flux is also negative and equal to W1:

.

Note in this case that S is integrated over the smaller surface. But the net electric and magnetic fields of both charges are computed prior to calculating S.

The energy flux over a surface containing only q2 is again positive and equal to W2. And the energy flux through a surface enclosing both charges is again equal to W1 + W2, but now (W1 + W2) < W2. Evidently some of the energy flowing away from q2 is "diverted" to q1 (and its driving agent); q1 is absorbing some of the energy emitted by q2.

Following is a Power Basic program that computes the W's and the field energy fluxes through spherical surfaces enclosing each charge, and a surface enclosing both. The computation of the electromagnetic fields at the charges and over the surfaces is relativistically rigorous. The values computed by the program are:

W1=-.142931399, Flux1=-.142831940

W2=18833.6599, Flux2=18833.6597

W=18833.5169, Flux=18833.5153

 

Program

#COMPILE EXE
#DIM ALL
FUNCTION PBMAIN () AS LONG
'*****************
'Two charges, q1=.0001 coulombs and q2=1 coulomb, oscillate along the x-axis.
'They are separated by one wavelength. The objective is to compute the
'work per cycle expended each cycle by an agent who counteracts the
'radiation reaction and interactive forces. Having done this, the field
'energy flux through enclosing surfaces is computed. Energy conservation is
'then confirmed.
'****************
'W1 is the work expended to drive q1.
'W2 is the work expended to drive q2.
'W is the work expended to drive both charges.
'Flux1 is the energy flux through a surface surrounding only q1.
'Flux2 is the energy flux through a surface surrounding only q2.
'Flux is the energy flux through a surface surrounding both charges.
'Compute W1, W2, and W.
'Physical and mathematical constants
DIM c AS DOUBLE
c= 299792000# 'Speed of light
DIM epsilon0 AS DOUBLE
epsilon0 = 0.00000000000885 'Permittivity constant
DIM pi AS DOUBLE
pi = 3.14159265358979
DIM Steps AS LONG
Steps = 100
DIM q1 AS DOUBLE
q1 = 0.0001
DIM q2 AS DOUBLE
q2 = 1
DIM Amp AS DOUBLE
Amp = 1 'Amplitude of each oscillation
DIM omega AS DOUBLE
omega = 0.01 * c / Amp 'Angular frequency of each oscillation
DIM freq AS DOUBLE
freq = omega / (2 * pi) 'frequency of each oscillation
DIM tau AS DOUBLE
tau = 1 / freq 'period of each oscillation
DIM lambda AS DOUBLE
lambda = c * tau 'wavelength
DIM deltat AS DOUBLE
deltat = tau / Steps 'time interval between epochs
DIM L AS DOUBLE
L = 1 * lambda 'separation
DIM dtheta AS DOUBLE
dtheta = pi / Steps 'Interval between angles over sphere
DIM Radius AS DOUBLE
Radius = L / 2 'Spherical radius of inner spheres
DIM BigRad AS DOUBLE
BigRad = L 'spherical radius of outer sphere
DIM W1 AS DOUBLE
DIM W2 AS DOUBLE
DIM W AS DOUBLE
DIM Flux1 AS DOUBLE
DIM Flux2 AS DOUBLE
DIM Flux AS DOUBLE
DIM x1(Steps) AS DOUBLE
DIM x2(Steps) AS DOUBLE 'current coordinates of q1 and q2
DIM v(Steps) AS DOUBLE
DIM gamma(Steps) AS DOUBLE
DIM a(Steps) AS DOUBLE
DIM dadt(Steps) AS DOUBLE
DIM Index AS LONG 'Loop index
DIM t(Steps) AS DOUBLE 'Current time
DIM tr AS DOUBLE 'Retarded time
DIM dt AS DOUBLE 't - tr
DIM dtmin AS DOUBLE 'Minimum value for dt
DIM dtmax AS DOUBLE 'Maximum value for dt
DIM xr1 AS DOUBLE
DIM xr2 AS DOUBLE 'Retarded positions (on x-axis)
DIM vrx1 AS DOUBLE
DIM vrx2 AS DOUBLE 'Retarded velocities
DIM arx1 AS DOUBLE
DIM arx2 AS DOUBLE 'Retarded accelerations
DIM Drx1 AS DOUBLE
DIM Drx2 AS DOUBLE
DIM Dry1 AS DOUBLE
DIM Dry2 AS DOUBLE 'components of vectors Dr
DIM Dr1 AS DOUBLE
DIM Dr2 AS DOUBLE 'magnitudes of vectors Dr
DIM ux1 AS DOUBLE
DIM ux2 AS DOUBLE
DIM uy1 AS DOUBLE
DIM uy2 AS DOUBLE 'components of vectors u
DIM Ex1(Steps) AS DOUBLE
DIM Ex2(Steps) AS DOUBLE 'x-components of electric field vectors
DIM Fx1(Steps) AS DOUBLE
DIM Fx2(Steps) AS DOUBLE
DIM Fx(Steps) AS DOUBLE 'Interact forces on q1 and q2
DIM FRadReact1(Steps) AS DOUBLE
DIM FRadReact2(Steps) AS DOUBLE 'Radiation Reaction forces
DIM FAgent1(Steps) AS DOUBLE
DIM FAgent2(Steps) AS DOUBLE 'Agent counteractions
DIM P1(Steps) AS DOUBLE
DIM P2(Steps) AS DOUBLE 'Agent expended powers
DIM theta AS DOUBLE 'Angle between x-axis and point on sphere (in xy-plane)
DIM Py AS DOUBLE
DIM Px AS DOUBLE 'Coordinates of point on sphere
DIM TimeIndex AS LONG
DIM ThetaIndex AS LONG 'Loop indexes
DIM FluxEx1 AS DOUBLE
DIM FluxEx2 AS DOUBLE
DIM FluxEx AS DOUBLE 'x-components of electric field vectors
DIM FluxEy1 AS DOUBLE
DIM FluxEy2 AS DOUBLE
DIM FLuxEy AS DOUBLE 'y-components of electric field vectors
DIM FluxBz1 AS DOUBLE
DIM FLuxBz2 AS DOUBLE
DIM FLuxBz AS DOUBLE 'z-components of magnetic field vectors
DIM Sx AS DOUBLE 'x-component of Poynting vector at point on sphere
DIM Sy AS DOUBLE 'y-component of Poynting vector
DIM unitx AS DOUBLE
DIM unity AS DOUBLE
DIM dA AS DOUBLE 'unit vectors and spherical area increment
DIM EFluxRate(Steps) AS DOUBLE
DIM EFlux(Steps) AS DOUBLE 'Energy flux at different angles
DIM Snormal AS DOUBLE 'normal component of Poynting vector
DIM ERad AS DOUBLE 'energy flux per cycle at separation L

FOR Index = 0 TO Steps - 1
'Compute the current positions and derivatives for q1 and q2.
t(Index) = Index * deltat
x1(Index) = -L / 2 + Amp * SIN(omega * t(Index))
x2(Index) = L / 2 + Amp * SIN(omega * t(Index))
v(Index) = omega * Amp * COS(omega * t(Index))
gamma(Index) = 1 / SQR(1 - v(Index) ^ 2 / c ^ 2)
a(Index) = -(omega ^ 2) * Amp * SIN(omega * t(Index))
dadt(Index) = -omega ^ 3 * Amp * COS(omega * t(Index))
'Compute the radiation reaction force on each charge.
FRadReact1(Index) = q1 ^ 2 / (6 * pi * epsilon0 * c ^ 3) * gamma(Index) ^ 4 * (dadt(Index) + 3 * gamma(Index) ^ 2 * v(Index) * a(Index) ^ 2 / c ^ 2)
FRadReact2(Index) = q2 ^ 2 / q1 ^ 2 * FRadReact1(Index)
dtmin = 0
dtmax = 5 * L / c
DO
dt = (dtmax + dtmin) / 2
tr = t(Index) - dt
xr1 = -L / 2 + Amp * SIN(omega * tr)
Drx1 = x2(Index) - xr1
Dr1 = ABS(Drx1)
IF ABS(c * dt - Dr1) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr1 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx1 = omega * Amp * COS(omega * tr)
arx1 = -(omega ^ 2) * Amp * SIN(omega * tr)
'Compute the needed retarded quantities for q2.
dtmin = 0
dtmax = 5 * L / c
DO
dt = (dtmax + dtmin) / 2
tr = t(Index) - dt
xr2 = L / 2 + Amp * SIN(omega * tr)
Drx2 = x1(Index) - xr2
Dr2 = ABS(Drx2)
IF ABS(c * dt - Dr2) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr2 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx2 = omega * Amp * COS(omega * tr)
arx2 = -(omega ^ 2) * Amp * SIN(omega * tr)
'Compute the components of vectors u.
ux1 = c * Drx1 / Dr1 - vrx1
ux2 = c * Drx2 / Dr2 - vrx2
'Compute the electric and magnetic field components
Ex1(Index) = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1) ^ 3 * (ux1 * (c ^ 2 - vrx1 ^ 2))
Ex2(Index) = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2) ^ 3 * (ux2 * (c ^ 2 - vrx2 ^ 2))
NEXT Index
'Now compute the electric (interactive) forces on q1 and q2.
FOR Index = 0 TO Steps - 1
Fx1(Index) = q1 * Ex2(Index)
Fx2(Index) = q2 * Ex1(Index)
NEXT Index
'Then compute W1, W2 and W.
FOR Index = 0 TO Steps - 1
FAgent1(Index) = -FRadReact1(Index) - Fx1(Index)
P1(Index) = FAgent1(Index) * v(Index)
FAgent2(Index) = -FRadReact2(Index) - Fx2(Index)
P2(Index) = FAgent2(Index) * v(Index)
NEXT Index
W1 = 0
W2 = 0
W = 0
FOR Index = 0 TO Steps - 1
W1 = W1 + P1(Index) * deltat
W2 = W2 + P2(Index) * deltat
NEXT Index
W = W1 + W2
MSGBOX ("W1 = " & STR$(W1) & "; W2 = " & STR$(W2) & "; W = " & STR$(W))
'Now compute the energy fluxes through the 3 surfaces.
'First compute Flux1.
ERad = 0
'Zero out the energy fluxes for each epoch.
FOR TimeIndex = 0 TO Steps - 1
EFlux(TimeIndex) = 0
NEXT TimeIndex
'Compute the energy flux through the spherical surface, for each time
'interval.
FOR TimeIndex = 0 TO Steps - 1 'time epochs
t(TimeIndex) = TimeIndex * deltat
'Zero out the energy flux rates for each theta interval.
FOR ThetaIndex = 0 TO Steps - 1
EFluxRate(ThetaIndex) = 0
NEXT ThetaIndex
'Compute the energy flux rate through each theta interval.
FOR ThetaIndex = 0 TO Steps - 1
theta = ThetaIndex * dtheta + dtheta / 2
unitx = COS(theta)
unity = SIN(theta)
Px = -L / 2 + Radius * COS(theta)
Py = Radius * SIN(theta)
dA = 2 * pi * Py * Radius * dtheta
'First do q1.
dtmin = 0
dtmax = 5 * Radius / c
DO
dt = (dtmax + dtmin) / 2
tr = t(TimeIndex) - dt
xr1 = -L / 2 + Amp * (SIN(omega * tr))
Drx1 = Px - xr1
Dry1 = Py
Dr1 = SQR(Drx1 ^ 2 + Dry1 ^ 2)
IF ABS(c * dt - Dr1) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr1 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx1 = Amp * omega * COS(omega * tr)
arx1 = -Amp * omega ^ 2 * SIN(omega * tr)
'Then do q2.
dtmin = 0
dtmax = 5 * Radius / c
DO
dt = (dtmax + dtmin) / 2
tr = t(TimeIndex) - dt
xr2 = L / 2 + Amp * SIN(omega * tr)
Drx2 = Px - xr2
Dry2 = Py
Dr2 = SQR(Drx2 ^ 2 + Dry2 ^ 2)
IF ABS(c * dt - Dr2) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr2 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx2 = Amp * omega * COS(omega * tr)
arx2 = -Amp * omega ^ 2 * SIN(omega * tr)
'Compute the components of vectors u.
ux1 = c * Drx1 / Dr1 - vrx1
uy1 = c * Dry1 / Dr1
ux2 = c * Drx2 / Dr2 - vrx2
uy2 = c * Dry2 / Dr2
'Compute the electric and magnetic field components
FluxEx1 = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (ux1 * (c ^ 2 - vrx1 ^ 2) + Dry1 * (-uy1 * arx1))
FluxEy1 = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (uy1 * (c ^ 2 - vrx1 ^ 2) - Drx1 * (-uy1 * arx1))
FluxBz1 = 1 / (c * Dr1) * (Drx1 * FluxEy1 - Dry1 * FluxEx1)
FluxEx2 = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (ux2 * (c ^ 2 - vrx2 ^ 2) + Dry2 * (-uy2 * arx2))
FluxEy2 = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (uy2 * (c ^ 2 - vrx2 ^ 2) - Drx2 * (-uy2 * arx2))
FLuxBz2 = 1 / (c * Dr2) * (Drx2 * FluxEy2 - Dry2 * FluxEx2)
'Sum the field components for the 2 charges.
FluxEx = FluxEx1 + FluxEx2
FLuxEy = FluxEy1 + FluxEy2
FLuxBz = FluxBz1 + FLuxBz2
'Compute the Poynting vector components, using the net E and B
'field components.
Sx = epsilon0 * c ^ 2 * (FLuxEy * FLuxBz)
Sy = epsilon0 * c ^ 2 * (-FluxEx * FLuxBz)
'Find the Poynting vector component that is normal to the
'spherical surface.
Snormal = Sx * unitx + Sy * unity
'Use it to update the energy flux from this value of theta.
EFluxRate(ThetaIndex) = EFluxRate(ThetaIndex) + Snormal * dA
'Repeat for the next angle theta.
NEXT ThetaIndex
'Now update the running sum of energy fluxes.
FOR ThetaIndex = 0 TO Steps - 1
EFlux(TimeIndex) = EFlux(TimeIndex) + EFluxRate(ThetaIndex) * deltat
NEXT ThetaIndex
NEXT TimeIndex
'Once the energy fluxes for each time interval have been computed, add
'the energy fluxes to Erad for the current separation.
FOR TimeIndex = 0 TO Steps - 1
ERad = ERad + EFlux(TimeIndex)
NEXT TimeIndex
Flux1 = ERad
'Display the results for comparison with W1.
MSGBOX ("Flux1 = " & STR$(Flux1) & "; Flux1/W1 = " & STR$(Flux1 / W1))

'Next compute Flux2.
ERad = 0
'Zero out the energy fluxes for each epoch.
FOR TimeIndex = 0 TO Steps - 1
EFlux(TimeIndex) = 0
NEXT TimeIndex
'Compute the energy flux through the spherical surface, for each time
'interval.
FOR TimeIndex = 0 TO Steps - 1 'time epochs
t(TimeIndex) = TimeIndex * deltat
'Zero out the energy flux rates for each theta interval.
FOR ThetaIndex = 0 TO Steps - 1
EFluxRate(ThetaIndex) = 0
NEXT ThetaIndex
'Compute the energy flux rate through each theta interval.
FOR ThetaIndex = 0 TO Steps - 1
theta = ThetaIndex * dtheta + dtheta / 2
unitx = COS(theta)
unity = SIN(theta)
Px = L / 2 + Radius * COS(theta)
Py = Radius * SIN(theta)
dA = 2 * pi * Py * Radius * dtheta
'First do q1.
dtmin = 0
dtmax = 5 * Radius / c
DO
dt = (dtmax + dtmin) / 2
tr = t(TimeIndex) - dt
xr1 = -L / 2 + Amp * (SIN(omega * tr))
Drx1 = Px - xr1
Dry1 = Py
Dr1 = SQR(Drx1 ^ 2 + Dry1 ^ 2)
IF ABS(c * dt - Dr1) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr1 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx1 = Amp * omega * COS(omega * tr)
arx1 = -Amp * omega ^ 2 * SIN(omega * tr)
'Then do q2.
dtmin = 0
dtmax = 5 * Radius / c
DO
dt = (dtmax + dtmin) / 2
tr = t(TimeIndex) - dt
xr2 = L / 2 + Amp * SIN(omega * tr)
Drx2 = Px - xr2
Dry2 = Py
Dr2 = SQR(Drx2 ^ 2 + Dry2 ^ 2)
IF ABS(c * dt - Dr2) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr2 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx2 = Amp * omega * COS(omega * tr)
arx2 = -Amp * omega ^ 2 * SIN(omega * tr)
'Compute the components of vectors u.
ux1 = c * Drx1 / Dr1 - vrx1
uy1 = c * Dry1 / Dr1
ux2 = c * Drx2 / Dr2 - vrx2
uy2 = c * Dry2 / Dr2
'Compute the electric and magnetic field components
FluxEx1 = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (ux1 * (c ^ 2 - vrx1 ^ 2) + Dry1 * (-uy1 * arx1))
FluxEy1 = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (uy1 * (c ^ 2 - vrx1 ^ 2) - Drx1 * (-uy1 * arx1))
FluxBz1 = 1 / (c * Dr1) * (Drx1 * FluxEy1 - Dry1 * FluxEx1)
FluxEx2 = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (ux2 * (c ^ 2 - vrx2 ^ 2) + Dry2 * (-uy2 * arx2))
FluxEy2 = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (uy2 * (c ^ 2 - vrx2 ^ 2) - Drx2 * (-uy2 * arx2))
FLuxBz2 = 1 / (c * Dr2) * (Drx2 * FluxEy2 - Dry2 * FluxEx2)
'Sum the field components for the 2 charges.
FluxEx = FluxEx1 + FluxEx2
FLuxEy = FluxEy1 + FluxEy2
FLuxBz = FluxBz1 + FLuxBz2
'Compute the Poynting vector components, using the net E and B
'field components.
Sx = epsilon0 * c ^ 2 * (FLuxEy * FLuxBz)
Sy = epsilon0 * c ^ 2 * (-FluxEx * FLuxBz)
'Find the Poynting vector component that is normal to the
'spherical surface.
Snormal = Sx * unitx + Sy * unity
'Use it to update the energy flux from this value of theta.
EFluxRate(ThetaIndex) = EFluxRate(ThetaIndex) + Snormal * dA
'Repeat for the next angle theta.
NEXT ThetaIndex
'Now update the running sum of energy fluxes.
FOR ThetaIndex = 0 TO Steps - 1
EFlux(TimeIndex) = EFlux(TimeIndex) + EFluxRate(ThetaIndex) * deltat
NEXT ThetaIndex
NEXT TimeIndex
'Once the energy fluxes for each time interval have been computed, add
'the energy fluxes to Erad for the current separation.
ERad = 0
FOR TimeIndex = 0 TO Steps - 1
ERad = ERad + EFlux(TimeIndex)
NEXT TimeIndex
Flux2 = ERad
'Display the results for comparison with W1.
MSGBOX ("Flux2 = " & STR$(Flux2) & "; Flux2/W2 = " & STR$(Flux2 / W2))

'Finally compute Flux.
ERad = 0
'Zero out the energy fluxes for each epoch.
FOR TimeIndex = 0 TO Steps - 1
EFlux(TimeIndex) = 0
NEXT TimeIndex
'Compute the energy flux through the spherical surface, for each time
'interval.
FOR TimeIndex = 0 TO Steps - 1 'time epochs
t(TimeIndex) = TimeIndex * deltat
'Zero out the energy flux rates for each theta interval.
FOR ThetaIndex = 0 TO Steps - 1
EFluxRate(ThetaIndex) = 0
NEXT ThetaIndex
'Compute the energy flux rate through each theta interval.
FOR ThetaIndex = 0 TO Steps - 1
theta = ThetaIndex * dtheta + dtheta / 2
unitx = COS(theta)
unity = SIN(theta)
Px = BigRad * COS(theta)
Py = BigRad * SIN(theta)
dA = 2 * pi * Py * BigRad * dtheta
'First do q1.
dtmin = 0
dtmax = 5 * BigRad / c
DO
dt = (dtmax + dtmin) / 2
tr = t(TimeIndex) - dt
xr1 = -L / 2 + Amp * (SIN(omega * tr))
Drx1 = Px - xr1
Dry1 = Py
Dr1 = SQR(Drx1 ^ 2 + Dry1 ^ 2)
IF ABS(c * dt - Dr1) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr1 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx1 = Amp * omega * COS(omega * tr)
arx1 = -Amp * omega ^ 2 * SIN(omega * tr)
'Then do q2.
dtmin = 0
dtmax = 5 * Radius / c
DO
dt = (dtmax + dtmin) / 2
tr = t(TimeIndex) - dt
xr2 = L / 2 + Amp * SIN(omega * tr)
Drx2 = Px - xr2
Dry2 = Py
Dr2 = SQR(Drx2 ^ 2 + Dry2 ^ 2)
IF ABS(c * dt - Dr2) < 2 ^ (-30) THEN EXIT DO
IF c * dt - Dr2 > 0 THEN
dtmax = dt
ELSE
dtmin = dt
END IF
LOOP
vrx2 = Amp * omega * COS(omega * tr)
arx2 = -Amp * omega ^ 2 * SIN(omega * tr)
'Compute the components of vectors u.
ux1 = c * Drx1 / Dr1 - vrx1
uy1 = c * Dry1 / Dr1
ux2 = c * Drx2 / Dr2 - vrx2
uy2 = c * Dry2 / Dr2
'Compute the electric and magnetic field components
FluxEx1 = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (ux1 * (c ^ 2 - vrx1 ^ 2) + Dry1 * (-uy1 * arx1))
FluxEy1 = q1 / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (uy1 * (c ^ 2 - vrx1 ^ 2) - Drx1 * (-uy1 * arx1))
FluxBz1 = 1 / (c * Dr1) * (Drx1 * FluxEy1 - Dry1 * FluxEx1)
FluxEx2 = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (ux2 * (c ^ 2 - vrx2 ^ 2) + Dry2 * (-uy2 * arx2))
FluxEy2 = q2 / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (uy2 * (c ^ 2 - vrx2 ^ 2) - Drx2 * (-uy2 * arx2))
FLuxBz2 = 1 / (c * Dr2) * (Drx2 * FluxEy2 - Dry2 * FluxEx2)
'Sum the field components for the 2 charges.
FluxEx = FluxEx1 + FluxEx2
FLuxEy = FluxEy1 + FluxEy2
FLuxBz = FluxBz1 + FLuxBz2
'Compute the Poynting vector components, using the net E and B
'field components.
Sx = epsilon0 * c ^ 2 * (FLuxEy * FLuxBz)
Sy = epsilon0 * c ^ 2 * (-FluxEx * FLuxBz)
'Find the Poynting vector component that is normal to the
'spherical surface.
Snormal = Sx * unitx + Sy * unity
'Use it to update the energy flux from this value of theta.
EFluxRate(ThetaIndex) = EFluxRate(ThetaIndex) + Snormal * dA
'Repeat for the next angle theta.
NEXT ThetaIndex
'Now update the running sum of energy fluxes.
FOR ThetaIndex = 0 TO Steps - 1
EFlux(TimeIndex) = EFlux(TimeIndex) + EFluxRate(ThetaIndex) * deltat
NEXT ThetaIndex
NEXT TimeIndex
'Once the energy fluxes for each time interval have been computed, add
'the energy fluxes to Erad for the current separation.
ERad = 0
FOR TimeIndex = 0 TO Steps - 1
ERad = ERad + EFlux(TimeIndex)
NEXT TimeIndex
Flux = ERad
'Display the results for comparison with W1.
MSGBOX ("Flux = " & STR$(Flux) & "; Flux/W = " & STR$(Flux / W))
MSGBOX('End of program')
END FUNCTION
'End Sub