Private Sub cmdNonRelOnly_Click()

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

'For a non-relativistically oscillating, small spherical shell of charge,

'Compute (1) F(t) and -q<Ex>, and (2) P(t) and Rate of Flux as a function of t.

'Output the values for plotting.

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

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 ThetaSteps As Long = 1000

Const Amp As Double = 1 'Amplitude of oscillation

Const omega As Double = 0.001 * c / Amp 'Angular frequency

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

Const tau As Double = 1 / freq 'Oscillation period

Const deltat As Double = tau / Steps 'Time increment

Const deltatfine As Double = deltat / 1000000

Const q As Double = 1 'Charge of one coulomb

Const Radius As Double = 0.005 'Spherical shell radius

Const sigma As Double = q / (4 * pi * Radius ^ 2) 'Surface charge density

Const dtheta As Double = pi / ThetaSteps 'Angle increment

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

Dim theta As Double 'Angle from x-axis to point on the spherical shell

Dim x As Double 'Present position of point chge at shell center

Dim v As Double 'Present velocity of pt chge and spherical shell

Dim gamma As Double

Dim dgamma4dt As Double

Dim a As Double 'Present acceleration of pt chge and spherical shell

Dim dadt As Double 'Time derivative of a

Dim Py As Double 'Present y-position of pt on spherical shell

Dim Px As Double 'Present x-position of pt on spherical shell

Dim i, j As Long 'Loop indexes

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 xr As Double 'Retarded position

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

Dim dardt As Double

Dim Drx As Double 'x-component of vector Dr

Dim Dry As Double 'y-component of vector Dr

Dim Dr As Double 'magnitude of vector Dr

Dim ux As Double 'x-component of vector u

Dim uy As Double 'y-component of vector u

Dim Ex As Double 'x-component of electric field vector

Dim Ey As Double 'y-component of electric field vector

Dim Bz As Double 'z-component of magnetic field vector

Dim Sx As Double 'x-component of Poynting vector

Dim Sy As Double 'y-component of Poynting vector

Dim Snormal As Double 'Component normal to spherical surface

Dim unitx, unity As Double 'Unit vectors

Dim dA As Double 'Area of annular stripe on spherical surface

Dim EFlux(Steps) As Double 'Rate of energy flux through spherical surface

Dim FTheory(Steps) As Double 'Force needed to drive spherical shell of chge

Dim FActual(Steps) As Double '-q<Ex>

Dim PTheory(Steps) As Double 'Rate that FTheory does work

Dim PActual(Steps) As Double ' Rate of energy flux through spherical surface

'Clear the Eflux array.

For i = 0 To Steps - 1

EFlux(i) = 0

Next i

'For each epoch in time ...

For i = 0 To Steps - 1

Debug.Print i

'Use this value of the time to compute things close to t=0.

't(i) = -Steps / 2 * deltatfine + i * deltatfine

'Compute the present variables.

t(i) = i * deltat

x = Amp * Sin(omega * t(i))

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

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

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

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

dgamma4dt = 4 * gamma ^ 3 * gamma ^ 3 * v * a / c ^ 2

'Then compute the driving force and the power it expends.

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

PTheory(i) = FTheory(i) * v

'Now compute -q<Ex> and the rate of energy flux, for comparison purposes.

FActual(i) = 0

PActual(i) = 0

'For each point on the spherical surface (in the xy-plane)...

For j = 0 To ThetaSteps - 1

'Compute the variables.

theta = j * dtheta + dtheta / 2

unitx = Cos(theta)

unity = Sin(theta)

Px = x + Radius * Cos(theta) 'Point on spherical surface

Py = Radius * Sin(theta)

dA = 2 * pi * Py * Radius * dtheta

'Then determine the retarded variables.

dtmin = 0

dtmax = 4 * (Amp + Radius) / c

Do

dt = (dtmax + dtmin) / 2

tr = t(i) - 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

'Then compute the fields at the point on the spherical surface.

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

FActual(i) = FActual(i) - Ex * sigma * dA

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)

'Compute the rate of energy flux through the surface area increment.

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

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

Snormal = Sx * unitx + Sy * unity

EFlux(i) = EFlux(i) + Snormal * dA

Next j

'Update the rate of energy flux through the entire spherical surface.

PActual(i) = EFlux(i)

Next i

'Output the data for use by the graphing software.

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

Open "c:\winmcad\Physics\FActual.prn" For Output As 2

For i = 0 To Steps - 1

'Use these two instructions to plot F(t) and -q<Ex>.

Write #1, t(i), FTheory(i)

Write #2, t(i), FActual(i)

'Use these two instructions to plot the rate of energy flux through

'the spherical surface, and the rate at which the driving force

'does work.

'Write #1, t(i), PTheory(i)

'Write #2, t(i), PActual(i)

Next i

Close

MsgBox ("Ready for plotting.")

End Sub

Private Sub cmdEnergyConChk_Click()

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

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

'total energy flux through a large spherical surface enclosing the

'oscillation site.

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

Const c As Double = 299792000#

Const epsilon0 As Double = 0.00000000000885

Const pi As Double = 3.14159265358979

Const Steps As Long = 100

Const Stepsp As Long = 1000000#

Const ThetaSteps As Long = 1000

Const Amp As Double = 1

Const omega As Double = 0.001 * c / Amp

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

Const tau As Double = 1 / freq

Const deltat As Double = tau / Steps

Const deltatp As Double = tau / Stepsp

Const q As Double = 1

Const BigRadius As Double = 5

Const Radius As Double = 0.005

Const dtheta As Double = pi / ThetaSteps

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

Dim theta As Double

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

Dim i, j As Long

Dim t(Steps), epoch As Double

Dim tr As Double

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 dardt 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(Steps) As Double

Dim FTheor, PTheor As Double

Dim WorkPerCycle, EfluxTotal As Double

Dim PActual(Steps) As Double

WorkPerCycle = 0

EfluxTotal = 0

For i = 0 To Steps - 1

EFlux(i) = 0

Next i

Debug.Print "A"

For i = 0 To Stepsp - 1

epoch = i * deltatp

x = Amp * Sin(omega * epoch)

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

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

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

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

dgamma4dt = 4 * gamma ^ 3 * gamma ^ 3 * v * a / c ^ 2

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

PTheor = FTheor * v

WorkPerCycle = WorkPerCycle + PTheor * deltatp

Next i

For i = 0 To Steps - 1

Debug.Print i

t(i) = i * deltat

x = Amp * Sin(omega * t(i))

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

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

PActual(i) = 0

For j = 0 To ThetaSteps - 1

theta = j * dtheta + dtheta / 2

unitx = Cos(theta)

unity = Sin(theta)

Px = BigRadius * Cos(theta)

Py = BigRadius * Sin(theta)

dA = 2 * pi * Py * BigRadius * dtheta

dtmin = 0

dtmax = 4 * (Amp + BigRadius) / c

Do

dt = (dtmax + dtmin) / 2

tr = t(i) - 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)

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

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

Snormal = Sx * unitx + Sy * unity

EFlux(i) = EFlux(i) + Snormal * dA

Next j

PActual(i) = EFlux(i)

EfluxTotal = EfluxTotal + PActual(i) * deltat

Next i

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

End Sub