Private Sub cmdComputeForces_Click()

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

'Compute the interactive forces acting on q1 and q2.

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

'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 L As Double = 0.5 * lambda 'separation

'Variables

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

Dim Index 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 'Interact forces on q1 and q2

For Index = 0 To Steps - 1

'Compute the current positions of q1 and q2.

t(Index) = Index * deltat

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

Px2 = L / 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 * L / c

Do

dt = (dtmax + dtmin) / 2

tr = t(Index) - dt

xr1 = -L / 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 * L / c

Do

dt = (dtmax + dtmin) / 2

tr = t(Index) - dt

xr2 = L / 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 and magnetic 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

'Output the times and interactive forces for plotting purposes.

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

For Index = 0 To Steps - 1

Write #1, t(Index), Fx1(Index), Fx2(Index), Fx(Index)

Next Index

Close 1

MsgBox ("Ready for plotting.")

End Sub