noxid100@cox.net

Option Explicit

Private Sub cmdFig1ab_Click()

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

'Compute (a) the energy flux through a surrounding spherical surface, vs. time,

'and (b) the energy flux per cycle vs. azimuth angle.

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const vmax As Double = 0.9999 * c 'Maximum speed

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda 'Radius of surrounding surface

Const tau As Double = 2 * pi / omega 'Period of oscillation

Const deltat As Double = tau / (Steps) 'Time between epochs

Const deltatheta As Double = pi / Steps 'Angle increment between surface rings

Const q As Double = 1 'Particle charge

'Variables

Dim i, j As Long 'Loop indexes

Dim t(Steps), theta(Steps) As Double 'Epochs and azimuth angles

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 xr As Double 'Retarded position

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

Dim Drx As Double 'x-component of vector Dr

Dim Dry As Double 'y-component of vector Dr

Dim Dr As Double 'Magnitude of vector Dr

Dim ux As Double 'x-component of vector u

Dim uy As Double 'y-component of vector u

Dim Ex As Double 'x-component of electric field vector

Dim Ey As Double 'y-component of electric field vector

Dim Bz As Double 'z-component of magnetic field vector

Dim Sx As Double 'x-component of Poynting vector

Dim Sy As Double 'y-component of Poynting vector

Dim dA As Double 'Area increment

Dim dAx As Double 'x-component of area increment

Dim dAy As Double 'y-component of area increment

Dim Px As Double 'x-component of point on sphere

Dim Py As Double 'y-component of point on sphere

Dim Eflux(Steps) As Double 'Energy fluxes at angles theta, particular epoch

Dim EfluxTime(Steps) As Double 'Energy flux, entire sphere, particular epoch

Dim EfluxAzimuth(Steps) As Double 'Flux per cycle at each theta

'Zero out the variables to be plotted vs. time and azimuth.

For i = 0 To Steps - 1

EfluxTime(i) = 0

EfluxAzimuth(i) = 0

Next i

'For each time epoch ...

For j = 0 To Steps - 1

Debug.Print j

t(j) = j * deltat

'And for each azimuth angle ...

For i = 0 To Steps - 1

'Compute the energy flux through the area between theta and theta+dtheta.

theta(i) = i * deltatheta + deltatheta / 2

Px = R * Cos(theta(i)) 'Point on spherical surface

Py = R * Sin(theta(i)) 'Ditto

Dry = Py 'Dry is independent of retarded position

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = t(j) - dt

xr = A * Sin(omega * tr)

Drx = Px - xr 'Drx depends on retarded position

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then

dtmax = dt

Else

dtmin = dt

End If

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ex = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (ux * (c ^ 2 - vr ^ 2) - Dry * uy * ar))

Ey = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (uy * (c ^ 2 - vr ^ 2) + Drx * uy * ar))

Bz = (1 / (Dr * c)) * (Drx * Ey - Dry * Ex)

Sx = epsilon0 * c ^ 2 * Ey * Bz

Sy = -epsilon0 * c ^ 2 * Ex * Bz

dA = 2 * pi * R * Sin(theta(i)) * R * deltatheta

dAx = dA * Cos(theta(i))

dAy = dA * Sin(theta(i))

Eflux(i) = (Sx * dAx + Sy * dAy) * deltat

'Update the running total energy flux per cycle at this angle.

EfluxAzimuth(i) = EfluxAzimuth(i) + Eflux(i)

Next i

'Sum all the energy fluxes through the area increments, to get the

'total energy flux through the sphere at the current time.

For i = 0 To Steps - 1

EfluxTime(j) = EfluxTime(j) + Eflux(i)

Next i

'Repeat for the next time epoch.

Next j

'Having computed the energy flux through the entire sphere, as a function

'of time, and the energy flux per cycle as a function of azimuth, output

'the values for use by the plotting utility.

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

For i = 0 To Steps - 1

Write #1, t(i), EfluxTime(i), theta(i), EfluxAzimuth(i)

Next i

Close

MsgBox ("Ready for plotting")

Stop

End Sub

Private Sub cmdFig2abc_Click()

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

'Compute the y-component of the electric field above x=R and in the range

'of y where the component is negative.

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const ymax As Double = 0.000264535 'Found by trial

Const deltay As Double = ymax / Steps

Const vmax As Double = 0.9999 * c 'Maximum speed

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda 'Radius of enveloping sphere

Const x As Double = R 'Middle of pulse

Const q As Double = 1 'Particle charge

'Variables

Dim i As Long 'Loop index

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 xr As Double 'Retarded position

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

Dim Drx As Double 'x-component of vector Dr

Dim Dry As Double 'y-component of vector Dr

Dim Dr As Double 'Magnitude of vector Dr

Dim ux As Double 'x-component of vector u

Dim uy As Double 'y-component of vector u

Dim y(Steps) As Double 'point above x-axis

Dim Ey(Steps) As Double 'y-component of electric field vector

'For each value of y above x=R and in the range of y over which Ey is negative ...

For i = 0 To Steps - 1

'Compute Ey.

y(i) = i * deltay

Dry = y(i)

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = A * Sin(omega * tr)

Drx = x - xr

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then dtmax = dt

If c * dt - Dr < 0 Then dtmin = dt

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ey(i) = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (uy * (c ^ 2 - vr ^ 2) + Drx * uy * ar))

Next i

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

For i = 0 To Steps - 1

Write #1, y(i), 0, Ey(i), 0 '0 positions y(i) correctly in vector for plotter

Next i

Close

MsgBox ("Ready for plotting")

Stop

End Sub

Private Sub cmdFig3ab_Click()

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

'Compute Ex and Bz over a pulsewidth of variable x, and at the point above

'x=R where Ey has its maximum negative value.

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

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

'Use this value for Figs. 3a and 3b.

'

'Const vmax As Double = 0.9999 * c 'Maximum speed of charge

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

'Use this value for maximum speed for Figs. 4a and 4b

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

Const vmax As Double = c

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda 'Center of pulse

Const pulsewidth As Double = 3.3770968684621E-05 'Found by trial

Const xmid As Double = R

Const xmin As Double = xmid - pulsewidth / 2

Const xmax As Double = xmid + pulsewidth / 2

Const dx As Double = (xmax - xmin) / Steps

Const Dry As Double = 0.00026454 'Radius of tiny cylinder

Const q As Double = 1

'Variables

Dim i As Long 'Loop index

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 xr As Double 'Retarded position

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

Dim Drx As Double 'x-component of vector Dr

Dim Dr As Double 'magnitude of vector Dr

Dim ux As Double 'x-component of vector u

Dim uy As Double 'y-component of vector u

Dim Ex(Steps), Ey(Steps), Bz(Steps) As Double

Dim Px As Double 'Point on sphere

Dim Py As Double 'Ditto

Dim x(Steps) As Double

'For each value of x ...

For i = 0 To Steps - 1

'Compute Ey and Bz.

x(i) = xmin + i * dx + dx / 2

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = A * Sin(omega * tr)

Drx = x(i) - xr

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then dtmax = dt

If c * dt - Dr < 0 Then dtmin = dt

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

'Ex(i) needed to compute Bz(i)

Ex(i) = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (ux * (c ^ 2 - vr ^ 2) - Dry * uy * ar))

Ey(i) = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (uy * (c ^ 2 - vr ^ 2) + Drx * uy * ar))

Bz(i) = 1 / (Dr * c) * (Drx * Ey(i) - Dry * Ex(i))

Next i

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

For i = 0 To Steps - 1

Write #1, x(i), Ey(i), Bz(i)

Next i

Close

MsgBox ("Ey and Bz ready to plot")

End Sub

Private Sub cmdEqn2abcd_Click()

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

'Compute the electric field fluxes through the small half cylinder end caps,

'through its side wall, and the total flux through the half cylinder.

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

'Physical and mathematical constants

Const Steps As Double = 50000 'Lots of iterations to minimize numerical errors

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const vmax As Double = 0.9999 * c 'Maximum speed of charge

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const R As Double = 10 * lambda 'Center of pulse at t=0

Const pulsewidth As Double = 3.3770968684621E-05

Const xmin As Double = R - pulsewidth / 2

Const xmax As Double = R 'Compute fluxes over left half of cylinder

Const dx As Double = (xmax - xmin) / Steps

Const ymax As Double = 0.00026454

Const dy As Double = ymax / Steps

Const q As Double = 1

'Variables

Dim i As Long 'Loop index

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 xr As Double 'Retarded position

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

Dim Drx As Double 'x-component of vector Dr

Dim Dry As Double 'y-component of vector Dr

Dim Dr As Double 'magnitude of vector Dr

Dim ux As Double 'x-component of vector u

Dim uy As Double 'y-component of vector u

Dim Ex As Double 'x-component of electric field vector

Dim Ey As Double 'y-component of electric field vector

Dim dA As Double 'Area increment

Dim dEflux As Double 'Increment of flux through area increment

Dim EfluxNetLeft, EfluxNetSide, EfluxNetRight, EfluxNet As Double 'Total fluxes

Dim x, y As Double 'Components of point on cylinder

EfluxNetLeft = 0 'Total electric field fluxes to be displayed

EfluxNetSide = 0

EfluxNetRight = 0

EfluxNet = 0

'Compute E Flux through left end cap.

x = xmin

'For each value of y on the end cap ...

For i = 0 To Steps - 1

'Compute Ex and then the electric field flux.

y = i * dy + dy / 2

Dry = y

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = A * Sin(omega * tr)

Drx = x - xr

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then dtmax = dt

If c * dt - Dr < 0 Then dtmin = dt

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ex = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (ux * (c ^ 2 - vr ^ 2) - Dry * uy * ar))

'Negative dA for flux into cylinder.

dA = -pi * ((y + dy / 2) ^ 2 - (y - dy / 2) ^ 2)

dEflux = Ex * dA

EfluxNetLeft = EfluxNetLeft + dEflux

Next i

MsgBox ("EfluxNetLeft = " & EfluxNetLeft)

'Compute E flux through side wall.

y = ymax

Dry = y

'For each x along the side wall ...

For i = 0 To Steps - 1

'Compute Ey and then the electric field flux.

x = xmin + i * dx + dx / 2

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = A * Sin(omega * tr)

Drx = x - xr

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then dtmax = dt

If c * dt - Dr < 0 Then dtmin = dt

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ey = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (uy * (c ^ 2 - vr ^ 2) + Drx * uy * ar))

dA = 2 * pi * y * dx

dEflux = Ey * dA

EfluxNetSide = EfluxNetSide + dEflux

Next i

MsgBox ("EfluxNetSide = " & EfluxNetSide)

'Compute E flux through right end cap.

x = R

'For each value of y on the end cap ...

For i = 0 To Steps - 1

'Compute Ex and then the electric field flux.

y = i * dy + dy / 2

Dry = y

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = A * Sin(omega * tr)

Drx = x - xr

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then dtmax = dt

If c * dt - Dr < 0 Then dtmin = dt

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ex = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (ux * (c ^ 2 - vr ^ 2) - Dry * uy * ar))

dA = pi * ((y + dy / 2) ^ 2 - (y - dy / 2) ^ 2)

dEflux = Ex * dA

EfluxNetRight = EfluxNetRight + dEflux

Next i

MsgBox ("EfluxNetright = " & EfluxNetRight)

'Sum the fluxes through the end caps and the side wall, and display the

'total flux through the cylinder.

EfluxNet = EfluxNetLeft + EfluxNetSide + EfluxNetRight

MsgBox ("Net E Flux = " & EfluxNet)

End Sub

Private Sub cmdFig5abc_Click()

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

'Compute Ex, Ey and Bz above x=R and over a larger range of y than that

'depicted in Fig. 2.

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

'Physical and mathematical constants

Const Steps As Double = 500 'Iterations

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const vmax As Double = 0.9999 * c 'Maximum speed

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = vmax / A 'Angular frequency

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

Const ymax As Double = lambda / 300

Const dy As Double = ymax / Steps

Const R As Double = 10 * lambda 'All points above x=R

Const x As Double = R

Const q As Double = 1 'Particle charge

'Variables

Dim i As Long 'Loop index

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 xr As Double 'Retarded position (on x-axis)

Dim vr As Double 'Retarded velocity

Dim ar As Double 'Retarded acceleration

Dim Drx As Double 'x-component of vector Dr

Dim Dry As Double

Dim Dr As Double 'magnitude of vector Dr

Dim ux As Double 'x-component of vector u

Dim uy As Double 'y-component of vector u

Dim y(Steps), Ex(Steps), Ey(Steps), Bz(Steps) As Double

'For each value of y ...

For i = 0 To Steps - 1

'Compute Ex, Ey, and Bz.

y(i) = i * dy + dy / 2

Dry = y(i)

dtmin = 0

dtmax = 5 * R / c

Do

dt = (dtmax + dtmin) / 2

tr = -dt

xr = A * Sin(omega * tr)

Drx = x - xr

Dr = Sqr(Drx ^ 2 + Dry ^ 2)

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

If c * dt - Dr > 0 Then dtmax = dt

If c * dt - Dr < 0 Then dtmin = dt

Loop

vr = omega * A * Cos(omega * tr)

ar = -(omega ^ 2) * A * Sin(omega * tr)

ux = c * Drx / Dr - vr

uy = c * Dry / Dr

Ex(i) = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (ux * (c ^ 2 - vr ^ 2) - Dry * uy * ar))

Ey(i) = q / (4 * pi * epsilon0) * ((Dr / (Drx * ux + Dry * uy) ^ 3) * (uy * (c ^ 2 - vr ^ 2) + Drx * uy * ar))

Bz(i) = 1 / (Dr * c) * (Drx * Ey(i) - Dry * Ex(i))

Next i

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

For i = 0 To Steps - 1

Write #1, y(i), Ex(i), Ey(i), Bz(i)

Next i

Close

MsgBox ("Ready for plotting")

Stop

End Sub