Option Explicit

Private Sub cmdComputeErads_Click()

'Computes the energy fluxes per cycle, and the values of Erad according to

'Larmor and Abraham-Lorentz. This data is used to populate Table 1-1 in the

'article.

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const dtheta As Double = pi / Steps 'Angle increment

'********

'Fill in vmax with the desired value prior to running.

'*******

Const vmax As Double = 0.01 * c 'Maximum speed of spherical shell

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda

Const tau As Double = 2 * pi / omega 'Period of oscillation

Const deltat As Double = tau / (Steps) 'Time between epochs

Const q As Double = 1 'Particle charge

Const Term As Double = q ^ 2 / (6 * pi * epsilon0 * c ^ 3)

'Variables

Dim i, j As Long 'Loop index

Dim t(Steps) As Double 'Time at which dEflux is computed

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 (on x-axis)

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

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 dA As Double 'Area increment

Dim dAx As Double 'x-component of area increment

Dim dAy As Double 'y-component of area increment

Dim theta(Steps) As Double 'Angle with x-axis

Dim Px As Double 'Point on sphere

Dim Py As Double 'Ditto

Dim dEflux As Double 'Increment to be added to Eflux(i)

Dim Eflux(Steps) As Double 'Energy fluxes at angles theta

Dim EfluxTime(Steps) 'Energy fluxes through entire sphere

Dim EnergyFlux As Double

Dim PRadLarmor(Steps), PRadAL(Steps) As Double

Dim v, accel, dadt As Double

Dim EradLarmor, EradAL As Double

'First compute the works per cycle, according to Larmor and Abraham-Lorentz,

'using the non-relativistic formulas for Prad.

For i = 0 To Steps - 1

t(i) = i * deltat

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

accel = -omega ^ 2 * A * Sin(omega * t(i))

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

PRadLarmor(i) = Term * accel ^ 2

PRadAL(i) = -Term * dadt * v

Next i

EradLarmor = 0

EradAL = 0

For i = 0 To Steps - 1

EradLarmor = EradLarmor + PRadLarmor(i) * deltat

EradAL = EradAL + PRadAL(i) * deltat

Next i

'Then compute the energy flux per cycle through the enclosing spherical

'surface.

For j = 0 To Steps - 1

Debug.Print j

t(j) = j * deltat

'Compute the energy flux through each slice between theta and theta+dtheta.

For i = 0 To Steps - 1

Eflux(i) = 0

theta(i) = i * dtheta + dtheta / 2

Px = R * Cos(theta(i))

Py = R * Sin(theta(i))

Dry = Py

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = t(j) - dt

xr = A * Sin(omega * tr)

Drx = Px - xr

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 * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * 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 / (Dr * c)) * (Drx * Ey - Dry * Ex)

Sx = epsilon0 * c ^ 2 * Ey * Bz

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

dA = 2 * pi * R * Sin(theta(i)) * R * dtheta

dAx = dA * Cos(theta(i))

dAy = dA * Sin(theta(i))

Eflux(i) = (Sx * dAx + Sy * dAy) * deltat

Next i

'Sum all the energy fluxes through the individual slices, to get the

'total energy flux through the sphere at time t.

EfluxTime(j) = 0

For i = 0 To Steps - 1

EfluxTime(j) = EfluxTime(j) + Eflux(i)

Next i

'Repeat for the next time.

Next j

EnergyFlux = 0

'Sum all of the energy fluxes for the various times, to find the energy

'flux per cycle through the spherical surface.

For j = 0 To Steps - 1

EnergyFlux = EnergyFlux + EfluxTime(j)

Next j

'Display the results for inclusion in Table 1-1.

MsgBox ("vmax = " & vmax / c & ", Erad(Larmor) = " & EradLarmor & ", Erad(AL) = " & EradAL & ", Energy Flux per Cycle = " & EnergyFlux)

End Sub

Private Sub cmdFigs1_1_Click()

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

'Compute the data for Figs. 1-1a and 1-1b

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const dtheta As Double = pi / Steps 'Angle increment

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

'Comment out all vmax constants except the one desired.

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

'Const vmax As Double = 0.9999 * c

'Const vmax As Double = c

Const vmax As Double = 0.01 * c 'Maximum speed of spherical shell

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda

Const tau As Double = 2 * pi / omega 'Period of oscillation

MsgBox ("tau = " & tau)

Const deltat As Double = tau / (Steps) 'Time between epochs

Const q As Double = 1 'Particle charge

Const Term As Double = q ^ 2 / (6 * pi * epsilon0 * c ^ 3)

'Variables

Dim i, j As Long 'Loop index

Dim t(Steps) As Double 'Time at which dEflux is computed

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 (on x-axis)

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

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 dA As Double 'Area increment

Dim dAx As Double 'x-component of area increment

Dim dAy As Double 'y-component of area increment

Dim theta(Steps) As Double 'Angle with x-axis

Dim Px As Double 'Point on sphere

Dim Py As Double 'Ditto

Dim dEflux As Double 'Increment to be added to Eflux(i)

Dim Eflux(Steps) As Double 'Energy fluxes at angles theta, one epoch

Dim EfluxNet(Steps) As Double 'Total energy flux at theta, complete cycle

Dim EfluxTime(Steps) 'Energy fluxes through entire sphere, one epoch

For j = 0 To Steps - 1

EfluxTime(i) = 0

EfluxNet(i) = 0

Next j

For j = 0 To Steps - 1

Debug.Print j

t(j) = j * deltat

For i = 0 To Steps - 1

Eflux(i) = 0

theta(i) = i * dtheta + dtheta / 2

Px = R * Cos(theta(i))

Py = R * Sin(theta(i))

Dry = Py

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = t(j) - dt

xr = A * Sin(omega * tr)

Drx = Px - xr

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 * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * 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 / (Dr * c)) * (Drx * Ey - Dry * Ex)

Sx = epsilon0 * c ^ 2 * Ey * Bz

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

dA = 2 * pi * R * Sin(theta(i)) * R * dtheta

dAx = dA * Cos(theta(i))

dAy = dA * Sin(theta(i))

dEflux = (Sx * dAx + Sy * dAy) * deltat

Eflux(i) = dEflux

EfluxNet(i) = EfluxNet(i) + Eflux(i)

Next i

For i = 0 To Steps - 1

EfluxTime(j) = EfluxTime(j) + Eflux(i)

Next i

Next j

Open "c:\WINMCAD\Physics\POSEfluxTime.prn" For Output As 1

For i = 0 To Steps - 1

Write #1, t(i), EfluxTime(i), theta(i), EfluxNet(i)

Next i

Close

MsgBox ("Energy fluxes computed.")

End Sub

Private Sub cmdTable2_1_Click()

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

'In Order to cut down on numerical errors, the number of steps in this

'program is raised from 500 to 5000. It takes a while to complete the

'program.

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

'Computes the energy fluxes per cycle, and the values of Erad according to

'Larmor and Abraham-Lorentz. This data is used to populate Table 1-1 in the

'article.

'Physical and mathematical constants

Const Steps As Double = 5000 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const dtheta As Double = pi / Steps 'Angle increment

'********

'Fill in vmax with the desired value prior to running.

'*******

Const vmax As Double = 0.99 * c 'Maximum speed of spherical shell

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda

Const tau As Double = 2 * pi / omega 'Period of oscillation

Const deltat As Double = tau / (Steps) 'Time between epochs

Const q As Double = 1 'Particle charge

Const Term As Double = q ^ 2 / (6 * pi * epsilon0 * c ^ 3)

'Variables

Dim i, j As Long 'Loop index

Dim t(Steps) As Double 'Time at which dEflux is computed

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 (on x-axis)

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

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 dA As Double 'Area increment

Dim dAx As Double 'x-component of area increment

Dim dAy As Double 'y-component of area increment

Dim theta(Steps) As Double 'Angle with x-axis

Dim Px As Double 'Point on sphere

Dim Py As Double 'Ditto

Dim dEflux As Double 'Increment to be added to Eflux(i)

Dim Eflux(Steps) As Double 'Energy fluxes at angles theta

Dim EfluxTime(Steps) 'Energy fluxes through entire sphere

Dim EnergyFlux As Double

Dim PRadLarmor(Steps), PRadAL(Steps) As Double

Dim v, accel, dadt, gamma, dgammadt As Double

Dim EradLarmor, EradAL As Double

'First compute the works per cycle, according to Larmor and Abraham-Lorentz,

'using the relativistic formulas for Prad.

For i = 0 To Steps - 1

t(i) = i * deltat

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

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

accel = -omega ^ 2 * A * Sin(omega * t(i))

dgammadt = gamma ^ 3 * v * accel / c ^ 2

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

PRadLarmor(i) = gamma ^ 6 * Term * accel ^ 2

PRadAL(i) = -Term * (gamma ^ 4 * dadt + accel * 3 / 4 * 4 * gamma ^ 3 * dgammadt) * v

Next i

EradLarmor = 0

EradAL = 0

For i = 0 To Steps - 1

EradLarmor = EradLarmor + PRadLarmor(i) * deltat

EradAL = EradAL + PRadAL(i) * deltat

Next i

'Then compute the energy flux per cycle through the enclosing spherical

'surface.

For j = 0 To Steps - 1

Debug.Print j

t(j) = j * deltat

'Compute the energy flux through each slice between theta and theta+dtheta.

For i = 0 To Steps - 1

Eflux(i) = 0

theta(i) = i * dtheta + dtheta / 2

Px = R * Cos(theta(i))

Py = R * Sin(theta(i))

Dry = Py

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = t(j) - dt

xr = A * Sin(omega * tr)

Drx = Px - xr

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 * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * 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 / (Dr * c)) * (Drx * Ey - Dry * Ex)

Sx = epsilon0 * c ^ 2 * Ey * Bz

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

dA = 2 * pi * R * Sin(theta(i)) * R * dtheta

dAx = dA * Cos(theta(i))

dAy = dA * Sin(theta(i))

Eflux(i) = (Sx * dAx + Sy * dAy) * deltat

Next i

'Sum all the energy fluxes through the individual slices, to get the

'total energy flux through the sphere at time t.

EfluxTime(j) = 0

For i = 0 To Steps - 1

EfluxTime(j) = EfluxTime(j) + Eflux(i)

Next i

'Repeat for the next time.

Next j

EnergyFlux = 0

'Sum all of the energy fluxes for the various times, to find the energy

'flux per cycle through the spherical surface.

For j = 0 To Steps - 1

EnergyFlux = EnergyFlux + EfluxTime(j)

Next j

'Display the results for inclusion in Table 1-1.

MsgBox ("vmax = " & vmax / c & ", Erad(Larmor) = " & EradLarmor & ", Erad(AL) = " & EradAL & ", Energy Flux per Cycle = " & EnergyFlux)

End Sub

Private Sub cmdFigs2_2_Click()

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

'Compute the Data for Figs. 2-2a and 2-2b.

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const vmax As Double = 0.9999 * c 'Maximum speed of spherical shell

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

Const tau As Double = 2 * pi / omega 'Period of oscillation

Const deltat As Double = tau / (Steps) 'Time between epochs

Const q As Double = 1 'Particle charge

Const Term As Double = q ^ 2 / (6 * pi * epsilon0 * c ^ 3)

'Variables

Dim i As Long 'Loop index

Dim t(Steps) As Double 'Time at which dEflux is computed

Dim v, accel, dadt, gamma, dgammadt As Double

Dim PRadLarmor(Steps) As Double

Dim PRadAL(Steps) As Double

For i = 0 To Steps - 1

t(i) = i * deltat

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

accel = -omega ^ 2 * A * Sin(omega * t(i))

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

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

dgammadt = gamma ^ 3 * v * accel / c ^ 2

PRadLarmor(i) = Term * gamma ^ 6 * accel ^ 2

PRadAL(i) = -Term * v * (gamma ^ 4 * dadt + 3 / 4 * accel * 4 * gamma ^ 3 * dgammadt)

Next i

Open "c:\WINMCAD\Physics\POSPrad.prn" For Output As #1

For i = 0 To Steps - 1

Write #1, t(i), PRadLarmor(i), PRadAL(i)

Next i

Close

MsgBox ("Ready for plotting PRad's")

End Sub

Private Sub cmdFigs2_3_Click()

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

'Compute the data for Figs. 2-3a and 2-3b and, with a change of vmax,

'for Figs. 2-4a and b.

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const dtheta As Double = pi / Steps 'Angle increment

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

'Comment out all but the desired vmax.

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

Const vmax As Double = 0.9999 * c 'Maximum speed of spherical shell

'Const vmax As Double = c

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda

Const tau As Double = 2 * pi / omega 'Period of oscillation

Const deltat As Double = tau / (Steps) 'Time between epochs

Const q As Double = 1 'Particle charge

Const Term As Double = q ^ 2 / (6 * pi * epsilon0 * c ^ 3)

'Variables

Dim i, j As Long 'Loop index

Dim t(Steps) As Double 'Time at which dEflux is computed

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 (on x-axis)

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

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 dA As Double 'Area increment

Dim dAx As Double 'x-component of area increment

Dim dAy As Double 'y-component of area increment

Dim theta(Steps) As Double 'Angle with x-axis

Dim Px As Double 'Point on sphere

Dim Py As Double 'Ditto

Dim dEflux As Double 'Increment to be added to Eflux(i)

Dim Eflux(Steps) As Double 'Energy fluxes at angles theta

Dim EfluxNet(Steps) As Double

Dim EfluxTime(Steps) 'Energy fluxes through entire sphere

For j = 0 To Steps - 1

EfluxTime(i) = 0

EfluxNet(i) = 0

Next j

For j = 0 To Steps - 1

Debug.Print j

t(j) = j * deltat

For i = 0 To Steps - 1

Eflux(i) = 0

theta(i) = i * dtheta + dtheta / 2

Px = R * Cos(theta(i))

Py = R * Sin(theta(i))

Dry = Py

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = t(j) - dt

xr = A * Sin(omega * tr)

Drx = Px - xr

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 * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * 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 / (Dr * c)) * (Drx * Ey - Dry * Ex)

Sx = epsilon0 * c ^ 2 * Ey * Bz

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

dA = 2 * pi * R * Sin(theta(i)) * R * dtheta

dAx = dA * Cos(theta(i))

dAy = dA * Sin(theta(i))

dEflux = (Sx * dAx + Sy * dAy) * deltat

Eflux(i) = dEflux

EfluxNet(i) = EfluxNet(i) + Eflux(i)

Next i

For i = 0 To Steps - 1

EfluxTime(j) = EfluxTime(j) + Eflux(i)

Next i

Next j

Open "c:\WINMCAD\Physics\POSEfluxTime.prn" For Output As 1

For i = 0 To Steps - 1

Write #1, t(i), EfluxTime(i), theta(i), EfluxNet(i)

Next i

Close

MsgBox ("Energy fluxes computed.")

End Sub