|
Appendix A Option Explicit Private Sub cmdCompute_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.95 * c 'Satellite speed in Kprime Const Rp As Double = 1 'Orbital radius in Kprime Const m0 As Double = 1 'Satellite rest mass Const vq As Double = 0.95 * c 'Central body speed in K Const taup As Double = 4 * pi * Rp / vp 'Orbital period in Kprime
'Variables Dim tau As Double 'Quasi-cycloid period in K Dim deltat As Double 'Time between motion updates Dim gammap As Double 'sqr(1-.95^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-.95ct (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 - 0.95 ^ 2) * taup 'Trial value for period in K deltat = tau / Steps gammap = 1 / Sqr(1 - 0.95 ^ 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-.95ct) '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 - 0.95 ^ 2) / (Sqr(1 - (0.95 * Sin(theta)) ^ 2)) ^ 3 Ey = C1 * Sin(theta) / R ^ 2 * (1 - 0.95 ^ 2) / (Sqr(1 - (0.95 * 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-.95^2) = " & Sqr(1 - 0.95 ^ 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-.95^2) = " & 1 / Sqr(1 - 0.95 ^ 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 / 10 - 1 Write #1, xRepo(10 * index), y(10 * index) Next index Close MsgBox ("Ready for plotting")
End Sub |