'**********

'Compute the y component of the Poynting Vector, at a point on the y-axis,

'at times throughout an oscillating charge cycle time.

'**********

'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 = 500 'Steps in an oscillation

Const A As Double = 1 'Amplitude of oscillation

Const omega As Double = 0.01 * c / A 'Angular frequency

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

Const lambda As Double = c * tau 'Wavelength

'Use following value of Py

Const Py As Double = 1 * lambda

'Const Py As Double = 0.05 * lambda 'y-component of observation point P

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

Const q As Double = 1 'Charge of one coulomb

'Variables

Dim i 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 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 '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(Steps) As Double 'x-component of electric field vector

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

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

Dim Sy(Steps) As Double 'y-component of the Poynting vector

'Turn on the hourglass.

Screen.MousePointer = vbHourglass

'Set up the epochs.

For i = 0 To Steps - 1

    t(i) = i * deltat

Next i

'Compute retarded times, positions, velocities and accelerations.

For i = 0 To Steps - 1

    dtmin = 0

    dtmax = Sqr(A ^ 2 + Py ^ 2) / c

    Do

        dt = (dtmax + dtmin) / 2

        tr = t(i) - dt

        xr = A * Sin(omega * tr)

        If Abs(c * dt - Sqr(xr ^ 2 + Py ^ 2)) < 2 ^ (-30) Then Exit Do

        If c * dt - Sqr(xr ^ 2 + Py ^ 2) > 0 Then

            dtmax = dt

        Else

            dtmin = dt

        End If

    Loop

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

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

    'Compute the components and magnitude of vector Dr.

    Drx = -xr

    Dry = Py

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

    'Compute the components of vector u.

    ux = c * Drx / Dr - vr

    uy = c * Dry / Dr

    'Compute the electric and magnetic field components

    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))

    'Compute the y-component of the Poynting Vector.

    Sy(i) = epsilon0 * c ^ 2 * (-Ex(i) * Bz(i))

Next i

'Write values to a file for plotting.

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

For i = 0 To Steps - 1

    Write #1, t(i), Sy(i)

Next i

Close 1

'Turn off the hourglass

Screen.MousePointer = vbArrow

MsgBox ("Ready for plotting.")

End Sub