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