,

.

If q1 = q2 then a net amount of work per cycle, W = W1+W2 = 2W1 must be done by a driving agent. 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 this 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 equals 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); q2 is absorbing some of the energy emitted by q1.

Private Sub cmdAbsorption_Click()

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

'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

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const Steps As Long = 100

Const q1 As Double = 0.0001

Const q2 As Double = 1

Const Amp As Double = 1 'Amplitude of each oscillation

Const omega As Double = 0.01 * c / Amp 'Angular frequency of each oscillation

Const freq As Double = omega / (2 * pi) 'frequency of each oscillation

Const tau As Double = 1 / freq 'period of each oscillation

Const lambda As Double = c * tau 'wavelength

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

Const L As Double = 1 * lambda 'separation

Const dtheta As Double = pi / Steps 'Interval between angles over sphere

Const Radius As Double = L / 2 'Spherical radius of inner spheres

Const BigRad As Double = L 'spherical radius of outer sphere

'Variables

Dim W1, W2, W As Double

Dim Flux1, Flux2, Flux As Double

Dim x1(Steps), x2(Steps) As Double 'current coordinates of q1 and q2

Dim v(Steps), gamma(Steps), a(Steps), 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, xr2 As Double 'Retarded positions (on x-axis)

Dim vrx1, vrx2 As Double 'Retarded velocities

Dim arx1, arx2 As Double 'Retarded accelerations

Dim Drx1, Drx2, Dry1, Dry2 As Double 'components of vectors Dr

Dim Dr1, Dr2 As Double 'magnitudes of vectors Dr

Dim ux1, ux2, uy1, uy2 As Double 'components of vectors u

Dim Ex1(Steps), Ex2(Steps) As Double 'x-components of electric field vectors

Dim Fx1(Steps), Fx2(Steps), Fx(Steps) As Double 'Interact forces on q1 and q2

Dim FRadReact1(Steps), FRadReact2(Steps) As Double 'Radiation Reaction forces

Dim FAgent1(Steps), FAgent2(Steps) As Double 'Agent counteractions

Dim P1(Steps), P2(Steps) As Double 'Agent expended powers

Dim theta As Double 'Angle between x-axis and point on sphere (in xy-plane)

Dim Py, Px As Double 'Coordinates of point on sphere

Dim TimeIndex, ThetaIndex As Long 'Loop indexes

Dim FluxEx1, FluxEx2, FluxEx As Double 'x-components of electric field vectors

Dim FluxEy1, FluxEy2, FLuxEy As Double 'y-components of electric field vectors

Dim FluxBz1, FLuxBz2, 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, unity, dA As Double 'unit vectors and spherical area increment

Dim EFluxRate(Steps), EFlux(Steps) As Double 'Enery 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

Debug.Print Index

'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)

'Compute the interactive forces.

'Compute the electric field of q1, at the position of q2, and

'the electric field of q2 at the position of q1.

'Compute the needed retarded quantities for q1.

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 = " & W1 & "; W2 = " & W2 & "; W = " & 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

Debug.Print TimeIndex

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 = " & Flux1 & "; Flux1/W1 = " & 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

Debug.Print TimeIndex

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 = " & Flux2 & "; Flux2/W2 = " & 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

Debug.Print TimeIndex

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 = " & Flux & "; Flux/W = " & Flux / W)

Stop

End Sub