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

'Compute the work per cycle expended to counteract the interactive

'forces on q1 and q2, at a range of separations.

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

'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 Steps As Long = 100

Const q As Double = 1

Const lambda As Double = 0.1 'wavelength

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

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

Const tau As Double = 1 / freq 'period

Const deltat As Double = tau / Steps 'time interval between epochs

Const Amp As Double = 0.1 * c / omega 'amplitude, in meters

Const LMin As Double = 0.5 * lambda 'minimum separation

Const Lmax As Double = 3 * lambda 'maximum separation

Const deltaL As Double = (Lmax - LMin) / Steps 'interval between separations

'Variables

Dim Px1, Px2 As Double 'coordinates of q1 and q2

Dim Index, LIndex As Long 'Loop index

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 xr1, xr2 As Double 'Retarded positions (on x-axis)

Dim vrx1, vrx2 As Double 'Retarded velocities

Dim arx1, arx2 As Double 'Retarded accelerations

Dim Drx1, Drx2 As Double 'x-components of vectors Dr

Dim Dr1, Dr2 As Double 'magnitudes of vectors Dr

Dim ux1, ux2 As Double 'x-components of vectors u

Dim Ex1(Steps), Ex2(Steps) As Double 'x-components of electric field vectors

Dim Fx1(Steps), Fx2(Steps), Fx(Steps) As Double 'interactive forces

Dim WInt(Steps) As Double 'Work per cycle

Dim L(Steps) As Double 'charge separation

Dim vx As Double 'present charge velocity

For LIndex = 0 To Steps - 1

WInt(LIndex) = 0

Next LIndex

'For each separation...

For LIndex = 0 To Steps - 1

Debug.Print LIndex

L(LIndex) = LMin + LIndex * deltaL

'Compute the work per cycle to counteract the interactive forces.

For Index = 0 To Steps - 1

'Compute the current positions of q1 and q2.

t(Index) = Index * deltat

Px1 = -L(LIndex) / 2 + Amp * Sin(omega * t(Index))

Px2 = L(LIndex) / 2 + Amp * Sin(omega * t(Index))

'Compute the electric field of q1, at the position of q2, and

'the electric field of q2 at the position of q1.

'Compute the needed retarded quantities for q1.

dtmin = 0

dtmax = 2 * Lmax / c

Do

dt = (dtmax + dtmin) / 2

tr = t(Index) - dt

xr1 = -L(LIndex) / 2 + Amp * Sin(omega * tr)

Drx1 = Px2 - xr1

Dr1 = Abs(Drx1)

If Abs(c * dt - Dr1) < 2 ^ (-30) Then Exit Do

If c * dt - Dr1 > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx1 = omega * Amp * Cos(omega * tr)

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

'Compute the needed retarded quantities for q2.

dtmin = 0

dtmax = 2 * Lmax / c

Do

dt = (dtmax + dtmin) / 2

tr = t(Index) - dt

xr2 = L(LIndex) / 2 + Amp * Sin(omega * tr)

Drx2 = Px1 - xr2

Dr2 = Abs(Drx2)

If Abs(c * dt - Dr2) < 2 ^ (-30) Then Exit Do

If c * dt - Dr2 > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx2 = omega * Amp * Cos(omega * tr)

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

'Compute the components of vectors u.

ux1 = c * Drx1 / Dr1 - vrx1

ux2 = c * Drx2 / Dr2 - vrx2

'Compute the electric field components

Ex1(Index) = q / (4 * pi * epsilon0) * Dr1 / ((Drx1 * ux1) ^ 3) * (ux1 * (c ^ 2 - vrx1 ^ 2))

Ex2(Index) = q / (4 * pi * epsilon0) * Dr2 / ((Drx2 * ux2) ^ 3) * (ux2 * (c ^ 2 - vrx2 ^ 2))

Next Index

'Now compute the electric (interactive) forces on q1 and q2.

For Index = 0 To Steps - 1

Fx1(Index) = q * Ex2(Index)

Fx2(Index) = q * Ex1(Index)

Fx(Index) = Fx1(Index) + Fx2(Index)

Next Index

'Compute the work per cycle expended to counteract the interactive

'forces, at the current separation.

For Index = 0 To Steps - 1

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

WInt(LIndex) = WInt(LIndex) - Fx(Index) * vx * deltat

Next Index

Next LIndex

'Output the works per cycle for plotting purposes.

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

For Index = 0 To Steps - 1

Write #1, L(Index) / lambda, WInt(Index)

Next Index

Close 1

MsgBox ("Ready for plotting.")

End Sub