Private Sub cmdLContractTDilation_Click()

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

'Given a charged satellite, circling an oppositely charged

'central body at constant speed, use Maxwell, Lorentz and

'Newton to compute the satellite's motion relative to another

'inertial frame. Output data reflecting the system shape

'for plotting purposes.

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

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

Const vp As Double = 0.01 * c 'Satellite speed in Kprime

‘Const vp as Double = .95 * c

Const vq As Double = vp 'Central body speed in K

Const Rp As Double = 1 'Orbital radius in Kprime

Const m0 As Double = 1 'Satellite rest mass

'Const taup As Double = 4 * pi * Rp / vp 'Orbital period in Kprime

Const taup As Double = 2 * pi * Rp / vp

'Variables

Dim tau As Double 'Quasi-cycloid period in K

Dim deltat As Double 'Time between motion updates

Dim gammap As Double '1/sqr(1-v2^2/c^2)

Dim q As Double 'Computed charge for circular motion in Kprime

Dim x As Double 'Computed satellite position in K, at time t

Dim xRepo(Steps) As Double 'x-vq*t (used to plot shape in K)

Dim y(Steps) As Double 'satellite y-position in K

Dim deltax, deltay As Double 'Displacements from central body in K

Dim R As Double 'Distance from central body to satellite, in K

Dim theta As Double 'Angle between x-axis and R vector

Dim vx, vy, v, gamma As Double 'Satellite variables

Dim ax, ay As Double 'Satellite variables

Dim t As Double 'Current time in K

Dim Ex As Double 'x-component of electric field vector at satellite in K

Dim Ey As Double 'y-component of electric field vector at satellite in K

Dim Bz As Double 'z-component of magnetic field vector at satellite in K

Dim Fx, Fy As Double 'Lorentz force acting on satellite in K

Dim Denom, Num As Double 'Used in solutions for ax and ay

Dim index As Long

Dim C1 As Double 'Saves compute time

Dim SecondMax As Boolean

Dim ymax As Double 'Second max value of satellite y position

'Try an educated trial value for the quasi-cycloid period in K.

tau = 1 / Sqr(1 - vq ^ 2 / c ^ 2) * taup 'Trial value for period in K

deltat = tau / Steps

gammap = 1 / Sqr(1 - vq ^ 2 / c ^ 2)

'Find what size charge will cause circular motion in Kprime.

q = Sqr(4 * pi * epsilon0 * Rp * gammap * m0 * vp ^ 2)

C1 = q / (4 * pi * epsilon0)

'Initialize satellite variables (for t=0) in K.

x = 0

y(0) = Rp

vx = 0

vy = 0

v = 0

'Then repeatedly compute the repositioned satellite x-position (x-vq*t)

'and y-position for shape plotting.

For index = 0 To Steps - 1

'1. Compute the fields.

t = index * deltat

'Use previously computed satellite positions.

xRepo(index) = x - vq * t

'R is always computed from the central body.

deltax = xRepo(index)

deltay = y(index)

'R might vary in K.

R = Sqr(deltax ^ 2 + deltay ^ 2)

'theta is the angle between the x-axis and the current satellite

'position, with the central body at the junction of the angle legs.

If deltax = 0 Then

If t = 0 Then

theta = pi / 2

Else

theta = 3 * pi / 2

End If

Else

If deltax < 0 And deltay >= 0 Then theta = pi / 2 + Atn(-deltax / deltay)

If deltax < 0 And deltay < 0 Then theta = pi + Atn(deltay / deltax)

If deltax >= 0 And deltay < 0 Then theta = 3 * pi / 2 + Atn(deltax / -deltay)

If deltax >= 0 And deltay >= 0 Then theta = Atn(deltay / deltax)

End If

v = Sqr(vx ^ 2 + vy ^ 2)

gamma = 1 / Sqr(1 - v ^ 2 / c ^ 2)

Ex = C1 * Cos(theta) / R ^ 2 * (1 - vq ^ 2 / c ^ 2) / (Sqr(1 - (vq / c * Sin(theta)) ^ 2)) ^ 3

Ey = C1 * Sin(theta) / R ^ 2 * (1 - vq ^ 2 / c ^ 2) / (Sqr(1 - (vq / c * Sin(theta)) ^ 2)) ^ 3

Bz = 1 / c ^ 2 * vq * Ey

'2. Compute the Lorentz force.

Fx = -q * (Ex + vy * Bz)

Fy = -q * (Ey - vx * Bz)

'3. Compute the acceleration.

'See the article for these formulas.

Num = c ^ 2 * (c ^ 2 + gamma ^ 2 * vy ^ 2) * Fx - gamma ^ 2 * c ^ 2 * vx * vy * Fy

Denom = (c ^ 2 + gamma ^ 2 * vy ^ 2) * (c ^ 2 * gamma * m0 + gamma ^ 3 * m0 * vx ^ 2) - gamma ^ 5 * m0 * vx ^ 2 * vy ^ 2

ax = Num / Denom

Num = c ^ 2 * (c ^ 2 + gamma ^ 2 * vx ^ 2) * Fy - gamma ^ 2 * c ^ 2 * vx * vy * Fx

Denom = (c ^ 2 + gamma ^ 2 * vx ^ 2) * (c ^ 2 * gamma * m0 + gamma ^ 3 * m0 * vy ^ 2) - gamma ^ 5 * m0 * vx ^ 2 * vy ^ 2

ay = Num / Denom

'4. Update the position and velocity.

x = x + vx * deltat + 1 / 2 * ax * deltat ^ 2

If index < Steps - 1 Then y(index + 1) = y(index) + vy * deltat + 1 / 2 * ay * deltat ^ 2

vx = vx + ax * deltat

vy = vy + ay * deltat

'5. Iterate.

Next index

'Show the xRepo value at which the satellite cuts the x-axis.

index = 0

Do

If Abs(y(index + 1)) < Abs(y(index)) And Abs(y(index + 1)) < Abs(y(index + 2)) Then

MsgBox ("xRepo value where satellite cuts x-axis = " & Abs(xRepo(index + 1)))

MsgBox ("xRepo/Rp = " & Abs(xRepo(index)) / Rp)

MsgBox ("sqr(1-vq^2/c^2) = " & Sqr(1 - vq ^ 2 / c ^ 2))

Exit Do

End If

index = index + 1

Loop

'Find out when the satellite again has its maximum positive value.

t = 0

index = 0

SecondMax = False

ymax = R

Do

If y(index + 1) > y(index) Then

ymax = y(index + 1)

tau = t + deltat

End If

index = index + 1

t = index * deltat

Loop Until index = Steps - 1

tau = t

MsgBox ("taup = " & taup & ", tau = " & tau & ", tau/taup = " & tau / taup)

MsgBox ("1/sqr(1-vq^2/c^2) = " & 1 / Sqr(1 - vq ^ 2 / c ^ 2))

'Output the xRepo and y satellite positions for shape plotting.

Open "c:\WINMCAD\Physics\XformElecDyn.prn" For Output As #1

For index = 0 To Steps / 1000 - 1

Debug.Print index

Write #1, xRepo(1000 * index), y(1000 * index)

Next index

Close

MsgBox ("Ready for plotting")

End Sub