Private Sub cmdComputeWorkNonRel

***Appendix, Visual Basic Programs for the article "Relativistic Reaction Forces" ***

Private Sub cmdComputeWorkNonRel_Click()

'**********

'Compute the work per cycle that an agent must do

'to nonrelativistically drive a charge sinusoidally.

'**********

'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 TableRows As Long = 10 'Number of data points or frequencies

Const Steps As Long = 1000 'Steps in interations for each frequency

Const A As Double = 1 'Amplitude of oscillation

Const q As Double = 1 'Charge

Const omegamin As Double = 0.01 * c 'Minimum frequency computed

Const omegamax As Double = 0.1 * c 'Maximum frequency computed

Const domega As Double = (omegamax - omegamin) / TableRows

'Variables.

Dim t, freq, tau, deltat As Double

Dim omega(TableRows), W(TableRows), P As Double

Dim F, v As Double

Dim iomega, itime As Long

'For each frequency ...

For iomega = 0 To TableRows - 1

omega(iomega) = omegamin + iomega * domega

freq = omega(iomega) / (2 * pi)

tau = 1 / freq

deltat = tau / Steps

W(iomega) = 0

'Compute the expended agent work per cycle.

For itime = 0 To Steps - 1

t = itime * deltat

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

F = -q ^ 2 / (6 * pi * epsilon0 * c ^ 3) * (-omega(iomega) ^ 3) * A * Cos(omega(iomega) * t)

P = F * v

W(iomega) = W(iomega) + P * deltat

Next itime

Next iomega

'Output the frequencies and works to a file.

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

For iomega = 0 To TableRows - 1

Write #1, omega(iomega), W(iomega)

Next iomega

MsgBox ("Ready for plotting")

End Sub

Private Sub cmdComputeFluxNonRel_Click()

'**********

'Compute the energy flux per cycle through an enclosing spherical surface

'of radius R=10, for a charge q=1 oscillating on the x-axis. Output the

'energy flux vs. the angular frequency for plotting purposes.

'**********

'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 TableRows As Long = 10 'Number of data points or frequencies

Const Steps As Long = 1000 'Iterations per frequency

Const A As Double = 1 'Amplitude of oscillation

Const q As Double = 1 'Charge

Const R As Double = 10 'Spherical surface radius

Const dtmaxinit As Double = 5 * (A + R) / c 'Initial value for dtmax

Const dphi As Double = pi / Steps 'Phi is angle with x-axis

Const omegamin As Double = 0.01 * c 'Minimum frequency

Const omegamax As Double = 0.1 * c 'Maximum frequency

Const domega = (omegamax - omegamin) / TableRows

Const epsfactor As Double = q / (4 * pi * epsilon0) 'Saves compute time

Const epscsq As Double = epsilon0 * c ^ 2 'Saves compute time

'Variables.

Dim iomega As Long

Dim omega(TableRows), tau, freq, dtau As Double

Dim t, tr As Double

Dim xr, vrx, arx As Double

Dim Px, Py, nx, ny As Double

Dim dAx, dAy, dA As Double

Dim dt As Double

Dim Eflux(TableRows) As Double

Dim phi As Double

Dim itime, iphi As Long

Dim dtmin As Double

Dim dtmax As Double

Dim Drx, Dry, Dr As Double

Dim DdotUcube As Double

Dim ux, uy As Double

Dim Ex, Ey As Double

Dim Bz As Double

Dim Sx, Sy As Double

Dim factor As Double

'For each frequency ...

For iomega = 0 To TableRows - 1

Eflux(iomega) = 0

omega(iomega) = omegamin + iomega * domega

freq = omega(iomega) / (2 * pi)

tau = 1 / freq

dtau = tau / Steps

'And for each epoch at that frequency ...

For itime = 0 To Steps - 1

t = itime * dtau

'And for each angle phi on the spherical surface ...

For iphi = 0 To Steps - 1

phi = iphi * dphi

'Compute the area of the band at that angle phi;

Px = R * Cos(phi + dphi / 2)

Py = R * Sin(phi + dphi / 2)

dA = R * dphi * 2 * pi * Py

nx = Px / R

ny = Py / R

dAx = dA * nx

dAy = dA * ny

'Compute the fields at that the center of that band;

dtmin = 0

dtmax = dtmaxinit

dt = (dtmax + dtmin) / 2

tr = t - dt

xr = A * Sin(omega(iomega) * tr)

Drx = Px - xr

Dry = Py

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

factor = c * dt - Dr

If Abs(factor) < 2 ^ (-30) Then Exit Do

If factor > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx = omega(iomega) * A * Cos(omega(iomega) * tr)

arx = -omega(iomega) ^ 2 * A * Sin(omega(iomega) * tr)

ux = c * Drx / Dr - vrx

uy = c * Dry / Dr

DdotUcube = (Drx * ux + Dry * uy) ^ 3

Ex = epsfactor * (Dr / DdotUcube) * (ux * (c ^ 2 - vrx ^ 2) + Dry * (-uy * arx))

Ey = epsfactor * (Dr / DdotUcube) * (uy * (c ^ 2 - vrx ^ 2) - Drx * (-uy * arx))

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

'Compute the energy flux for an increment of time through that band

Sx = epscsq * (Ey * Bz)

Sy = epscsq * (-Ex * Bz)

Eflux(iomega) = Eflux(iomega) + (Sx * dAx + Sy * dAy) * dtau

Next iphi

Next itime

Next iomega

'Output the frequencies and energy fluxes per cycle.

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

For iomega = 0 To TableRows - 1

Write #1, omega(iomega), Eflux(iomega)

Next iomega

MsgBox ("Ready for plotting")

End Sub

Private Sub cmdGenTable_Click()

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

'Compute the percent differences, using the work and energy flux per

'cycle computed above. Output frequency, work per cycle, energy flux

'per cycle and percent difference for copying into article table.

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

Const TableRows As Long = 10

Dim iomega As Long

Dim omega, work, flux, PercentDiff As Double

Open "c:\WINMCAD\Physics\QWork.prn" For Input As #1

Open "c:\WINMCAD\Physics\QFlux.prn" For Input As #2

Open "c:\WINMCAD\Physics\QTAble.prn" For Output As #3

For iomega = 0 To TableRows - 1

Input #1, omega, work

Input #2, omega, flux

PercentDiff = (flux - work) / work

Write #3, omega, work, flux, PercentDiff

Next iomega

MsgBox ("Table data ready.")

End Sub

Private Sub cmdComputeWorkRel_Click()

'**********

'Compute the work per cycle that an agent must do

'to relativistically drive a charge sinusoidally.

'**********

'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 TableRows As Long = 10 'Number of data points or frequencies

Const Steps As Long = 1000 'Steps in interations for each frequency

Const A As Double = 1 'Amplitude of oscillation

Const q As Double = 1 'Charge

Const omegamin As Double = 0.85 * c 'Minimum frequency computed

Const omegamax As Double = 0.95 * c 'Maximum frequency computed

Const domega As Double = (omegamax - omegamin) / TableRows

'Variables.

Dim t, freq, tau, deltat As Double

Dim omega(TableRows), W(TableRows), P As Double

Dim F, v, gamma, ax, daxdt As Double

Dim iomega, itime As Long

'For each frequency ...

For iomega = 0 To TableRows - 1

omega(iomega) = omegamin + iomega * domega

freq = omega(iomega) / (2 * pi)

tau = 1 / freq

deltat = tau / Steps

W(iomega) = 0

'Compute the expended agent work per cycle.

For itime = 0 To Steps - 1

t = itime * deltat

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

ax = -omega(iomega) ^ 2 * A * Sin(omega(iomega) * t)

daxdt = -omega(iomega) ^ 3 * A * Cos(omega(iomega) * t)

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

F = -gamma ^ 4 * q ^ 2 / (6 * pi * epsilon0 * c ^ 3) * (daxdt + 3 * gamma ^ 2 * v * ax ^ 2 / c ^ 2)

P = F * v

W(iomega) = W(iomega) + P * deltat

Next itime

Next iomega

'Output the frequencies and works to a file.

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

For iomega = 0 To TableRows - 1

Write #1, omega(iomega), W(iomega)

Next iomega

MsgBox ("Ready for plotting")

End Sub