Private Sub cmd1Charge_Click()

'Mathematical and physical constants.

Const pi As Double = 3.14159265

Const c As Double = 299792000 'speed of light

Const eps0 As Double = 0.00000000000885 'permittivity constant

Const steps As Long = 100 'Number of time epochs in one cycle

Const deltatheta As Double = pi / steps 'Step size for angle theta

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

'User-definable constants follow. All but the desired value can

'be commented out, or other values can be substituted. The only

'requirements are that each maximum speed be less than .001c and that

'the distance between the charge and points on the surrounding spherical

'surface be >= a half wavelength at all times.

'Const vmax As Double = 0.001 * c 'Maximum charge speed (m/sec)

'Const vmax As Double = 0.0005 * c

Const vmax As Double = 5000

Const freq As Double = 1000000# 'Frequency of oscillation (Hz)

'Const freq As Double = 5000000#

'Const freq As Double = 800000#

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

'Derived constants and other constants follow.

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

Const Amp As Double = vmax / omega 'Amplitude of oscillation (meters)

MsgBox ("Amp = " & Amp)

Const period As Double = 1 / freq 'Period of oscillation (sec)

Const deltat As Double = period / steps 'Time interval between time epochs (sec)

Const lambda As Double = c / freq 'Wavelength of any emitted radiation (meters)

Const q As Double = 0.001 'Electric charge (coul)

Const R As Double = lambda 'Radius of enclosing spherical surface (meters)

'Variables follow.

Dim theta As Double 'Angle subtended by positive x-axis and point on spherical surface

Dim Ex, Ey, Bz As Double 'Fields at point on spherical surface

Dim Sx, Sy, SNormal As Double 'Poynting vector, S, at point on spherical surface

Dim t(steps), tr, dt As Double 'Current and retarded time

Dim ux, uy As Double

Dim Dx, Dy, D As Double 'Displacement between retarded position and field eval point

Dim indext, indextheta As Long 'For loop indexes

Dim xr, vrx, vr, arx As Double 'Retarded quantities

Dim dtmin, dtmax As Double 'Trial values of dt

Dim StripRadius, StripArea As Double 'Strip on spherical surface, concentric to x-axis

Dim Eflux(steps), NetEFlux(steps) As Double 'Energy fluxes

Dim Erad, W As Double 'Energy flux per cycle time through spherical surface, and agent work per cycle

Dim FRadReact As Double 'Radiation reaction force

Dim Px, Py As Double 'Field evaluation point

Dim Term As Double 'Saves compute time

Dim vx(steps), ax(steps), daxdt(steps) As Double 'Pre-computed quantites

'Set up some useful quantities.

For indext = 0 To steps - 1

t(indext) = indext * deltat

vx(indext) = omega * Amp * Cos(omega * t(indext))

daxdt(indext) = -omega ^ 3 * Amp * Cos(omega * t(indext))

Next indext

'For each time epoch...

For indext = 0 To steps - 1

Debug.Print indext

'...and for each angle theta...

For indextheta = 0 To steps - 1

theta = indextheta * deltatheta + deltatheta / 2

'...compute the energy flux through the strip at theta.

Px = R * Cos(theta)

Py = R * Sin(theta)

StripRadius = Py

StripArea = 2 * pi * StripRadius * R * deltatheta

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmin + dtmax) / 2

tr = t(indext) - dt 'Trial retarded time

xr = Amp * Sin(omega * tr) 'Trial retarded x

Dx = Px - xr

Dy = Py

D = Sqr(Dx ^ 2 + Dy ^ 2)

If Abs(c * dt - D) < 2 ^ -30 Then Exit Do 'Test for best retarded time

If (c * dt - D) > 0 Then dtmax = dt

If (c * dt - D) < 0 Then dtmin = dt

Loop

vrx = omega * Amp * Cos(omega * tr)

vr = vrx

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

ux = c * Dx / D - vrx

uy = c * Dy / D

Term = q / (4 * pi * eps0) * D / ((Dx * ux + Dy * uy) ^ 3)

Ex = Term * (ux * (c ^ 2 - vr ^ 2) + Dy * (-uy * arx))

Ey = Term * (uy * (c ^ 2 - vr ^ 2) - Dx * (-uy * arx))

Bz = 1 / (c * D) * (Dx * Ey - Dy * Ex)

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

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

SNormal = Sx * Cos(theta) + Sy * Sin(theta)

Eflux(indextheta) = SNormal * StripArea * deltat

Next indextheta

'Then sum all of the strip energy fluxes to get the energy flux through the entire sphere.

NetEFlux(indext) = 0

For indextheta = 0 To steps - 1

NetEFlux(indext) = NetEFlux(indext) + Eflux(indextheta)

Next indextheta

Next indext

'After computing the energy fluxes between each pair of time epochs, sum them to get the

'energy flux through the sphere over the course of an entire oscillation time.

Erad = 0

For indext = 0 To steps - 1

Erad = Erad + NetEFlux(indext)

Next indext

'Display the result for inclusion in the appropriate article table.

MsgBox ("ERad = " & Erad)

'Now compute the work per cycle expended to counteract the radiation reaction force.

W = 0

For indext = 0 To steps - 1

FRadReact = q ^ 2 / (6 * pi * eps0 * c ^ 3) * daxdt(indext)

W = W - FRadReact * vx(indext) * deltat

Next indext

'Display this result also, for inclusion in the article table.

MsgBox ("W = " & W)

Stop

End Sub