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

'Given an oscillating spherical shell of charge,

'Compute the total work per cycle done by the driving force, and the

'total energy flux per cycle through a spherical surface enclosing the

'oscillation site. Display the two results for comparison.

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

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 = 1000

Const ThetaSteps As Long = 100

Const Amp As Double = 1 'amplitude

Const omega As Double = 0.01 * c / Amp 'angular frequency

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

Const tau As Double = 1 / freq 'period

Const deltat As Double = tau / Steps 'time increment

Const q As Double = 1 'shell charge

Const BigRadius As Double = 10 * Amp 'radius of sphere for energy flux

Const Radius As Double = 0.005 'radius of spherical shell

Const dtheta As Double = pi / ThetaSteps 'step in azimuth

Const m As Double = q ^ 2 / (6 * pi * epsilon0 * Radius * c ^ 2) 'elecmag mass

Dim theta As Double

Dim x, v, gamma, a, dadt, Py, Px As Double

Dim i, j As Long 'indexes

Dim t As Double

Dim tr As Double 'retarded time

Dim dt As Double

Dim dtmin As Double

Dim dtmax As Double

Dim xr As Double

Dim vr As Double

Dim ar As Double

Dim Drx As Double

Dim Dry As Double

Dim Dr As Double

Dim ux As Double

Dim uy As Double

Dim Ex As Double

Dim Ey As Double

Dim Bz As Double

Dim Sx As Double

Dim Sy As Double

Dim Snormal As Double

Dim unitx, unity, dA As Double

Dim EFlux As Double

Dim F, P As Double

Dim WorkPerCycle As Double

WorkPerCycle = 0

EFlux = 0

'For each time epoch...

For i = 0 To Steps - 1

'...compute the work done by the agent force during that time increment...

t = i * deltat

x = Amp * Sin(omega * t)

v = omega * Amp * Cos(omega * t)

a = -omega ^ 2 * Amp * Sin(omega * t)

dadt = -omega ^ 3 * Amp * Cos(omega * t)

gamma = 1 / Sqr(1 - v ^ 2 / c ^ 2)

F = gamma ^ 3 * m * a - gamma ^ 4 * q ^ 2 / (6 * pi * epsilon0 * c ^ 3) * (dadt + 3 * gamma ^ 2 * v * a ^ 2 / c ^ 2)

P = F * v

'...and update the running total.

WorkPerCycle = WorkPerCycle + P * deltat

Next i

'Then for each time epoch...

For i = 0 To Steps - 1

Debug.Print i

t = i * deltat

'...and for each azimuth angle

For j = 0 To ThetaSteps - 1

'...set up for flux dot product with area increment...

theta = j * dtheta + dtheta / 2

unitx = Cos(theta)

unity = Sin(theta)

Px = BigRadius * Cos(theta) 'point of Poynting evauation

Py = BigRadius * Sin(theta)

dA = 2 * pi * Py * BigRadius * dtheta 'area of strip

'...and find the fields.

dtmin = 0

dtmax = 4 * (Amp + BigRadius) / c

Do

dt = (dtmax + dtmin) / 2

tr = t - dt

xr = Amp * 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 * Amp * Cos(omega * tr)

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

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ex = q / (4 * pi * epsilon0) * Dr / (Drx * ux + Dry * uy) ^ 3 * (ux * (c ^ 2 - vr ^ 2) - Dry * uy * ar)

Ey = q / (4 * pi * epsilon0) * Dr / (Drx * ux + Dry * uy) ^ 3 * (uy * (c ^ 2 - vr ^ 2) + Drx * uy * ar)

Bz = 1 / (c * Dr) * (Drx * Ey - Dry * Ex)

'Then compute the components of the Poynting vector...

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

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

'...and form the dot product of it with the area strip

Snormal = Sx * unitx + Sy * unity

'Update the total energy flux through the spherical surface.

EFlux = EFlux + Snormal * dA * deltat

Next j

Next i

MsgBox ("W = " & WorkPerCycle & ", Erad =" & EFlux)

Stop

End Sub