Private Sub cmdVarRad_Click()

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

'Compute the x-component of the spherical shell electric field

'at its center when the shell has variable radii.

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

Const c As Double = 299792000

Const eps0 As Double = 0.00000000000885

Const pi As Double = 3.14159265358979

Const Steps As Long = 25

Const Steps1 As Long = 2000

Const q As Double = -1 'charge of spherical shell

Const deltatheta As Double = pi / Steps1 'increments in azimuth angle

Const accel As Double = 1000000#

Const Rmin As Double = 1

Const Rmax As Double = 10

Const dR As Double = (Rmax - Rmin) / Steps

Dim sigma As Double

Dim R(Steps) As Double

Dim Ex(Steps) As Double 'x-component of E at origin

Dim Index, Index1 As Long 'counters

Dim dq As Double 'charge of annulus

Dim x, y, z As Double 'present quantities

Dim xr, yr, zr, vrx, vry, vrz, arx, ary, arz As Double 'retarded quantities

Dim Drx, Dry, Drz, D As Double 'displacement from dq to origin

Dim ux, uy, uz As Double

Dim theta As Double 'azimuth of dq

Dim dtmin, dtmax, dt, tr As Double

'For each radius...

For Index = 0 To Steps - 1

Debug.Print Index

R(Index) = Rmin + Index * dR

sigma = q / (4 * pi * R(Index) ^ 2)

Ex(Index) = 0

'Then for each azimuth angle at that radius...

For Index1 = 0 To Steps1 - 1

theta = Index1 * deltatheta + deltatheta / 2

x = R(Index) * Cos(theta)

y = R(Index) * Sin(theta)

z = 0

dq = 2 * pi * R(Index) * Sin(theta) * R(Index) * deltatheta * sigma

'...compute Ex attributable to the pt chge,

'and add the result to the running total for Ex.

dtmin = 0

dtmax = 5 * R(Index) / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = x + 1 / 2 * accel * tr ^ 2

yr = y

zr = z

Drx = -xr

Dry = -yr

Drz = -zr

D = Sqr(Drx ^ 2 + Dry ^ 2 + Drz ^ 2)

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

If c * dt - D > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx = -accel * tr

vry = 0

vrz = 0

arx = accel

ary = 0

arz = 0

ux = c * Drx / dR - vrx

uy = c * Dry / dR - vry

uz = c * Drz / dR - vrz

Ex(Index) = Ex(Index) + dq / (4 * pi * eps0) * D / (Drx * ux + Dry * uy + Drz * uz) ^ 3 * (ux * (c ^ 2 - vrx ^ 2) + Dry * (ux * ary - uy * arx) - Drz * (uz * arx - ux * arz))

Next Index1

Next Index

'After Ex for each radius has been computed, output

'the values for plotting purposes.

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

For Index = 0 To Steps - 1

Write #1, R(Index), Ex(Index)

Next Index

Close

MsgBox ("Ready for plotting")

End

End Sub

Private Sub cmdModified_Click()

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

'Compute the x-component of the spherical shell electric field

'at its center when the shell has constant accelerations.

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

Const c As Double = 299792000

Const eps0 As Double = 0.00000000000885

Const pi As Double = 3.14159265358979

Const Steps As Long = 25

Const Steps1 As Long = 2000

Const q As Double = -1 'charge of spherical shell

Const Radius As Double = 1

Const deltatheta As Double = pi / Steps1 'increments in azimuth angle

Const sigma As Double = q / (4 * pi * Radius ^ 2) 'shell surface chge density

Const accelmax As Double = 1000000#

Const daccel As Double = accelmax / Steps

Dim accel(Steps) As Double

Dim Ex(Steps) As Double 'x-component of E at origin

Dim Index, Index1 As Long 'counters

Dim dq As Double 'charge of annulus

Dim x, y, z As Double 'present quantities

Dim xr, yr, zr, vrx, vry, vrz, arx, ary, arz As Double 'retarded quantities

Dim Drx, Dry, Drz, D As Double 'displacement from dq to origin

Dim ux, uy, uz As Double

Dim theta As Double 'azimuth of dq

Dim dtmin, dtmax, dt, tr As Double

'For each acceleration...

For Index = 0 To Steps - 1

Debug.Print Index

accel(Index) = Index * daccel

Ex(Index) = 0

'Then for each azimuth angle during that acceleration...

For Index1 = 0 To Steps1 - 1

theta = Index1 * deltatheta + deltatheta / 2

x = Radius * Cos(theta)

y = Radius * Sin(theta)

z = 0

dq = 2 * pi * Radius * Sin(theta) * Radius * deltatheta * sigma

'...compute Ex attributable to the pt chge,

'and add the result to the running total for Ex.

dtmin = 0

dtmax = 5 * Radius / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = x + 1 / 2 * accel(Index) * tr ^ 2

yr = y

zr = z

Drx = -xr

Dry = -yr

Drz = -zr

D = Sqr(Drx ^ 2 + Dry ^ 2 + Drz ^ 2)

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

If c * dt - D > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vrx = -accel(Index) * tr

vry = 0

vrz = 0

arx = accel(Index)

ary = 0

arz = 0

ux = c * Drx / D - vrx

uy = c * Dry / D - vry

uz = c * Drz / D - vrz

Ex(Index) = Ex(Index) + dq / (4 * pi * eps0) * D / (Drx * ux + Dry * uy + Drz * uz) ^ 3 * (ux * (c ^ 2 - vrx ^ 2) + Dry * (ux * ary - uy * arx) - Drz * (uz * arx - ux * arz))

Next Index1

Next Index

'After Ex for each acceleration has been computed, output

'the values for plotting purposes.

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

For Index = 0 To Steps - 1

Write #1, accel(Index), Ex(Index)

Next Index

Close

MsgBox ("Ready for plotting")

End

End Sub