Option Explicit

Private Sub cmd8PtChges_Click()

Const N As Double = 90 'Number of time epochs

Const pi As Double = 3.141592654

Const eps0 As Double = 0.00000000000885

Const c As Double = 299792000

Const R As Double = 1 'Spherical shell radius

Const lambda As Double = 100 * R 'Suggested wavelength of any emitted radiation

Const nu As Double = c / lambda 'Frequency

Const omega As Double = 2 * pi * nu 'Angular frequency

Const vmax As Double = 0.001 * c 'Maximum, non-relativistic speed

Const Amp As Double = vmax / omega 'Amplitude of oscillation

Const q As Double = 0.001 'Charge of individual, quasi-point charge

Const Totalq As Double = 6 * q 'Charge of entire spherical shell

Const period As Double = 1 / nu 'Oscillation period

Const deltat As Double = period / N 'Time between epochs

Const mElecMag As Double = Totalq ^ 2 / (6 * pi * eps0 * R * c ^ 2) 'Electromgnetic mass of shell

Dim xs(6), ys(6), zs(6) As Double 'Positions at t=0

Dim timeindex, srcindex, subindex As Long 'Loop counters

Dim t(N), dt, tr, dtmax, dtmin As Double 'Present and retarded time variables

Dim sourcex, sourcey, sourcez As Double 'Present position of field source charge

Dim subjectx, subjecty, subjectz As Double 'Present position of acted upon charge

Dim xr, yr, zr As Double 'Retarded position of source charge

Dim Dx, Dy, Dz, D As Double 'Vector from source retarded position to subject charge present position

Dim vrx, vry, vrz, vr As Double 'Retarded velocity

Dim arx, ary, arz As Double 'Retarded acceleration

Dim ux, uy, uz As Double 'Working vector in point charge field solution

Dim Ex As Double 'Electric field at subject charge

Dim F(N) As Double 'Total force experienced by all 6 quasi-point charges at time N*deltat

Dim FTheory(N) As Double 'Theoretical total force by spherical shell at time N*deltat

Dim TotalReact, InertReact, RadReact As Double 'Theoretical reaction forces

Dim SphereCenter(N) As Double 'Position of spherical shell center at time N*deltat

'Set up the quasi-point charge coordinates at time t=0.

xs(0) = R

ys(0) = 0

zs(0) = 0

xs(1) = 0

ys(1) = 0

zs(1) = -R

xs(2) = -R

ys(2) = 0

zs(2) = 0

xs(3) = 0

ys(3) = 0

zs(3) = R

xs(4) = 0

ys(4) = R

zs(4) = 0

xs(5) = 0

ys(5) = -R

zs(5) = 0

'Initialize variable arrays that will be output for plotting purposes.

For timeindex = 0 To N - 1

t(timeindex) = timeindex * deltat

SphereCenter(timeindex) = Amp * Sin(omega * t(timeindex))

F(timeindex) = 0

Next timeindex

'For each time epoch ...

For timeindex = 0 To N - 1

'... compute the interactive forces;

For srcindex = 0 To 5

sourcex = SphereCenter(timeindex) + xs(srcindex)

sourcey = ys(srcindex)

sourcez = zs(srcindex)

yr = sourcey 'Source retarded y and z are same as current y and z

zr = sourcez

For subindex = 0 To 5

If subindex = srcindex Then GoTo SameCharge 'Skip when source and subject are the same

subjectx = SphereCenter(timeindex) + xs(subindex) 'All subject ring chges have same x-coord

subjecty = ys(subindex)

subjectz = zs(subindex)

dtmin = 0

dtmax = (5 * R) / c

Do 'Retarded quantitites depend on subject pt chge position

dt = (dtmin + dtmax) / 2

tr = t(timeindex) - dt

xr = sourcex + Amp * Sin(omega * tr)

Dx = subjectx - xr

Dy = subjecty - yr

Dz = subjectz - zr

D = Sqr(Dx ^ 2 + Dy ^ 2 + Dz ^ 2)

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

If (c * dt - D) > 0 Then dtmax = dt

If (c * dt - D) < 0 Then dtmin = dt

Loop

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

vry = 0

vrz = 0

vr = Sqr(vrx ^ 2 + vry ^ 2 + vrz ^ 2)

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

ary = 0

arz = 0

ux = c * Dx / D - vrx

uy = c * Dy / D - vry

uz = c * Dz / D - vrz

Ex = q / (4 * pi * eps0) * D / ((Dx * ux + Dy * uy + Dz * uz) ^ 3) * (ux * (c ^ 2 - vr ^ 2) + (Dy * (ux * ary - uy * arx) - Dz * (uz * arx - ux * arz)))

F(timeindex) = F(timeindex) + q * Ex

SameCharge: Next subindex

Next srcindex

'Now compute the theoretical force experienced by the spherical shell.

InertReact = -mElecMag * (-omega ^ 2 * Amp * Sin(omega * t(timeindex))) 'Newton third

RadReact = Totalq ^ 2 / (6 * pi * eps0 * c ^ 3) * (-omega ^ 3 * Amp * Cos(omega * t(timeindex))) 'Ab/Lor

TotalReact = InertReact + RadReact 'Newton and Abraham/Lorentz

FTheory(timeindex) = TotalReact

Next timeindex

'Output the time, theoretical and computed net forces for plotting purposes.

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

For timeindex = 0 To N - 1

Write #1, t(timeindex) / period, FTheory(timeindex), F(timeindex)

Next timeindex

Close

MsgBox ("Ready for plotting")

Stop

End Sub