Private Sub cmdComputeW_Click()

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

'Compute W(D) and output the results for plotting purposes.

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

Const pi As Double = 3.14159265

Const c As Double = 299792000

Const eps0 As Double = 0.00000000000885

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

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

Const Amp As Double = 1

Const omega As Double = vmax / Amp

Const freq As Double = omega / 2 / pi

Const period = 1 / freq

Const deltat As Double = period / steps

Const q As Double = 1

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

Const deltaSep As Double = (3 * lambda - lambda) / steps

Const x01 = -lambda / 2

Const x02 = lambda / 2

Dim Sep(steps) As Double

Dim Ex As Double

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

Dim ux As Double

Dim Dx, D As Double

Dim indext, indexSep As Long 'For loop indexes

Dim xr, vrx, vr, arx As Double

Dim dtmin, dtmax As Double 'Trial values of dt

Dim W(steps) As Double

Dim FRadReact1(steps), FRadReact2(steps) As Double

Dim FInteract1(steps), FInteract2(steps) As Double

Dim NetF1(steps), NetF2(steps) As Double

Dim Px, Py As Double

Dim Term As Double

Dim vx(steps), daxdt(steps) As Double

Dim WIsolated As Double

Dim SepFirstMin, SepFirstMax As Double

'Precompute some useful variables.

For indext = 0 To steps - 1

t(indext) = indext * deltat

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

daxdt(indext) = -omega ^ 2 * vx(indext)

Sep(indext) = x02 - x01 + indext * deltaSep

Next indext

'Compute the work per cycle needed to drive both charges when they

'are infinitely separated.

WIsolated = 2 * q ^ 2 * Amp ^ 2 * omega ^ 3 / (6 * eps0 * c ^ 3)

'Compute the radiation reaction forces, which are the same at all Separations.

For indext = 0 To steps - 1

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

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

Next indext

'Then for each charge separation ...

For indexSep = 0 To steps - 1

Debug.Print indexSep

'And at each time epoch ...

For indext = 0 To steps - 1

'... compute the Interactive forces, first on the right hand charge,

'and then on the left.

Px = Sep(indexSep) / 2 + Amp * Sin(omega * t(indext))

dtmin = 0

dtmax = 8 * lambda / c

Do

dt = (dtmin + dtmax) / 2

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

xr = -Sep(indexSep) / 2 + Amp * Sin(omega * tr) 'Trial retarded x

Dx = Px - xr

D = Sqr(Dx ^ 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

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

Ex = Term * (ux * (c ^ 2 - vr ^ 2))

FInteract2(indext) = Ex * q

Px = -Sep(indexSep) / 2 + Amp * Sin(omega * t(indext))

dtmin = 0

dtmax = 8 * lambda / c

Do

dt = (dtmin + dtmax) / 2

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

xr = Sep(indexSep) / 2 + Amp * Sin(omega * tr) 'Trial retarded x

Dx = Px - xr

D = Sqr(Dx ^ 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

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

Ex = Term * (ux * (c ^ 2 - vr ^ 2))

FInteract1(indext) = Ex * q

Next indext

'Sum the RadReact and Interactive forces. The negative is the force

'the driving agent must exert at each time epoch.

For indext = 0 To steps - 1

NetF1(indext) = -FRadReact1(indext) - FInteract1(indext)

NetF2(indext) = -FRadReact2(indext) - FInteract2(indext)

Next indext

'Compute the net work per cycle that the driving agent must do at the

'current charge separation.

W(indexSep) = 0

For indext = 0 To steps - 1

W(indexSep) = W(indexSep) + NetF1(indext) * vx(indext) * deltat

W(indexSep) = W(indexSep) + NetF2(indext) * vx(indext) * deltat

Next indext

Next indexSep

'Find the separations at which W is maximum, and after that the separation

'at which W returns to a minimum.

indexSep = 0

Do

indexSep = indexSep + 1

If W(indexSep) > W(indexSep - 1) Then

SepFirstMax = lambda + (indexSep) * deltaSep

End If

Loop Until W(indexSep) < W(indexSep - 1)

MsgBox ("First maximum occurs at " & SepFirstMax / lambda & " wavelengths.")

Do

indexSep = indexSep + 1

Loop Until W(indexSep) > W(indexSep - 1)

SepFirstMin = lambda + (indexSep - 1) * deltaSep

MsgBox ("Second minimum occurs at " & SepFirstMin / lambda & " wavelengths.")

'Output the values for plotting.

Open "c:\WINMCAD\Physics\Sect8.PRN" For Output As #1

For indexSep = 0 To steps - 1

Write #1, Sep(indexSep) / lambda, W(indexSep), WIsolated

Next indexSep

Close

MsgBox ("Ready for plotting")

Stop

End Sub