'**********

'Compute the energy flux per cycle through a spherical surface surrounding

'an oscillating pair of charges, separated by a distance L that varies

'from .5 lambda to 3 lambda.

'**********

'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 q As Double = 1

Const lambda As Double = 0.1 'wavelength (meters)

Const omega As Double = 2 * pi * c / lambda 'angular frequency

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

Const tau As Double = 1 / freq 'period

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

Const Amp As Double = 0.1 * c / omega 'amplitude (meters)

Const LMin As Double = 0.5 * lambda 'minimum separation

Const Lmax As Double = 3 * lambda 'maximum separation

Const deltaL As Double = (Lmax - LMin) / Steps 'interval between separations

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

Const Radius As Double = 2 * (Lmax / 2 + Amp) 'spherical radius

'Variables

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, LIndex 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 As Double 'x-components of vectors Dr

Dim Dry1, Dry2 As Double 'y-components of vectors Dr

Dim Dr1, Dr2 As Double 'magnitudes of vectors Dr

Dim ux1, ux2 As Double 'x-components of vectors u

Dim uy1, uy2 As Double 'y-components of vectors u

Dim Ex1, Ex2, Ex As Double 'x-components of electric field vectors

Dim Ey1, Ey2, Ey As Double 'y-components of electric field vectors

Dim Bz1, Bz2, Bz 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 L(Steps) As Double 'separation of q1 and q2

Dim EradL(Steps) As Double 'energy flux per cycle at separation L

Dim RadWork As Double 'Work per cycle to counteract rad react forces

Dim ERadInt(Steps) As Double 'Total work minus RadWork

'Zero out each value of Erad(L)

For LIndex = 0 To Steps - 1

EradL(LIndex) = 0

Next LIndex

'Compute Erad(L) for each separation from the Minimum to the Maximum L.

For LIndex = 0 To Steps - 1

Debug.Print LIndex

'Set the separation for this loop.

L(LIndex) = LMin + LIndex * deltaL

'Zero out the energy fluxes for each epoch, present separation.

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 = Radius * Cos(theta)

Py = Radius * Sin(theta)

dA = 2 * pi * Py * Radius * dtheta

'First do q1.

dtmin = 0

dtmax = 2 * Radius / c

Do

dt = (dtmax + dtmin) / 2

tr = t(TimeIndex) - dt

xr1 = -L(LIndex) / 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 = omega * Amp * Cos(omega * tr)

arx1 = -(omega ^ 2) * Amp * Sin(omega * tr)

'Then do q2.

dtmin = 0

dtmax = 2 * Radius / c

Do

dt = (dtmax + dtmin) / 2

tr = t(TimeIndex) - dt

xr2 = L(LIndex) / 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 = omega * Amp * Cos(omega * tr)

arx2 = -(omega ^ 2) * Amp * 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

Ex1 = q / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (ux1 * (c ^ 2 - vrx1 ^ 2) + Dry1 * (-uy1 * arx1))

Ey1 = q / (4 * pi * epsilon0) * Dr1 / (Drx1 * ux1 + Dry1 * uy1) ^ 3 * (uy1 * (c ^ 2 - vrx1 ^ 2) - Drx1 * (-uy1 * arx1))

Bz1 = 1 / (c * Dr1) * (Drx1 * Ey1 - Dry1 * Ex1)

Ex2 = q / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (ux2 * (c ^ 2 - vrx2 ^ 2) + Dry2 * (-uy2 * arx2))

Ey2 = q / (4 * pi * epsilon0) * Dr2 / (Drx2 * ux2 + Dry2 * uy2) ^ 3 * (uy2 * (c ^ 2 - vrx2 ^ 2) - Drx2 * (-uy2 * arx2))

Bz2 = 1 / (c * Dr2) * (Drx2 * Ey2 - Dry2 * Ex2)

'Sum the field components for the 2 charges.

Ex = Ex1 + Ex2

Ey = Ey1 + Ey2

Bz = Bz1 + Bz2

'Compute the Poynting vector components, using the net E and B

'field components.

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

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

'Find the Poynting vector component that is normal to the

'spherical surface.

Snormal = Sx * unitx + Sy * unity

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 has been computed, add

'the energy fluxes to Erad for the current separation.

For TimeIndex = 0 To Steps - 1

EradL(LIndex) = EradL(LIndex) + EFlux(TimeIndex)

Next TimeIndex

Next LIndex

'Having computed Erad(L), optionally subtract out the (constant) work expended

'to counteract the radiation reaction force.

RadWork = q ^ 2 * omega ^ 3 * Amp ^ 2 / (3 * epsilon0 * c ^ 3)

For LIndex = 0 To Steps - 1

ERadInt(LIndex) = EradL(LIndex) - RadWork

Next LIndex

'After computing all the energy fluxes over the entire range of L,

'output the values for plotting purposes.

Open "c:\winmcad\ERad.prn" For Output As 1

For LIndex = 0 To Steps - 1

Write #1, L(LIndex) / lambda, EradL(LIndex), ERadInt(LIndex)

Next LIndex

Close 1

MsgBox ("Ready for plotting.")

End Sub