Programming for the "Snapshot" Algorithm

Option Explicit

Private Sub cmdTransform_Click()

Const c As Double = 300000000#

Const u As Double = 0.95 * c 'No distortion in Kprime when set to zero

Const L As Double = 2

'Use this value of omega to compute simple length contraction.

Const omega As Double = 0

'Use this value to compute both length contraction and bending.

'Const omega As Double = 50000

Const N As Long = 100 'Number of parts in right half of rod

Const M As Long = 10000# 'Number of dtprimes

Const deltax As Double = L / 2 / N 'Length of part

Dim xp(N) As Double 'xprime coordinates of parts

Dim yp(N) As Double 'yprime coordinates of parts

Dim t As Double 'K time

Dim x As Double 'K kinematic variables

Dim y As Double

Dim vx As Double

Dim vy As Double

Dim ax As Double

Dim ay As Double

Dim gamma As Double

Dim tp As Double 'Kprime time

Dim deltatp As Double

Dim dtp As Double 'Incremental period for updating Kprime coords.

Dim vxp As Double 'Kprime kinematic variables

Dim vyp As Double

Dim axp As Double

Dim ayp As Double

Dim MIndex As Long

Dim NIndex As Long

Dim hypot As Double 'Triangle hypotenuse

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

'For each part of the rod right half ...

For NIndex = 0 To N - 1

Debug.Print NIndex

t = 0 'Initial time

'...Compute the K space coords;

hypot = (NIndex + 1 / 2) * deltax

x = hypot

y = 0

'compute the amount of Kprime time you must work through

deltatp = -gamma * (t - u * x / c ^ 2)

dtp = deltatp / M

'Then for each part...

For MIndex = 0 To M - 1

'...Compute the K variables;

x = hypot * Cos(omega * t)

y = hypot * Sin(omega * t)

vx = -omega * y

vy = omega * x

ax = -omega ^ 2 * x

ay = omega ^ 2 * y

'Transform the K variables to Kprime variables

tp = gamma * (t - u * x / c ^ 2)

xp(NIndex) = gamma * (x - u * t)

yp(NIndex) = y

vxp = (vx - u) / (1 - u * vx / c ^ 2)

vyp = vy / gamma / (1 - u * vx / c ^ 2)

axp = ax / gamma ^ 3 / (1 - u * vx / c ^ 2) ^ 3

ayp = ay / gamma ^ 2 / (1 - u * vx / c ^ 2) ^ 2 + u * vy * ax / c ^ 2 / gamma ^ 2 / (1 - u * vx / c ^ 2) ^ 3

'Compute follow-on Kprime variables (assuming constant acceleration)

xp(NIndex) = xp(NIndex) + vxp * dtp + 1 / 2 * axp * dtp ^ 2

yp(NIndex) = yp(NIndex) + vyp * dtp + 1 / 2 * axp * dtp ^ 2

tp = tp + dtp

'update K time and iterate.

t = gamma * (tp + u * xp(NIndex) / c ^ 2)

Next MIndex

Next NIndex

'Output the Kprime x and y coordinates for plotting.

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

For NIndex = 0 To N - 1

Write #1, xp(NIndex), yp(NIndex)

Next NIndex

Close

MsgBox ("Ready for plotting")

Stop

End Sub