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

'Compute the orbital points for a negatively charged satellite,

'orbiting a fixed, positively charged central body.

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

'Constants follow.

Const Orbits As Long = 6 'Compute the equivalent of 6 circular orbits

Const c As Double = 299792000# 'Speed of light

Const epsilon0 As Double = 0.00000000000885 'Permittivity constant

Const pi As Double = 3.14159265358979

Const m As Double = 1

Const q As Double = 1

'**********

'Select whichever circular orbit speed is appropriate.

Const circv As Double = 0.001 * c

'Const circv As Double = 0.05 * c

Const steps As Long = 100000 'Granularity of computation

Const alpha As Double = q ^ 2 / (4 * pi * epsilon0) 'Saves compute time

'Variables follow.

Dim CircR As Double 'Radius of Circular Orbit

Dim omega As Double 'Angular orbital frequency of Circular Orbit

Dim freq As Double 'Frequency of Circular Orbit

Dim tau As Double 'Period of Circular Orbit

Dim deltat As Double ‘Time between epochs

Dim x(Orbits * steps), y(Orbits * steps) As Double 'Computed Orbital Points

Dim vx, vxhalf, vy, vyhalf, v, gamma, ax, ay As Double 'Computed kinematics

Dim F, Fx, Fy As Double 'Electrostatic Force on orbiting particle

Dim D As Double 'Distance from central body to satellite

Dim Denom, Num As Double 'Used in computing acceleration components

Dim i As Long 'index

Dim UCirc, TCirc, Ecirc, UElip, TElip As Double 'Energies

Dim vmax As Double 'Maximum speed attained

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

CircR = q ^ 2 / (4 * pi * epsilon0 * gamma * m * circv ^ 2) 'Circular Orbit Radius

UCirc = -alpha / CircR 'Circular orbit potential energy

TCirc = m * (gamma - 1) * c ^ 2 'CIrcular orbit kinetic energy

Ecirc = UCirc + TCirc 'Circular orbit total energy

omega = circv / CircR 'Circular orbit anglar rate

freq = omega / (2 * pi) 'Circular orbit frequency

tau = 1 / freq 'Circular orbit period

deltat = tau / steps 'Time between epochs

x(0) = 1.9 * CircR 'Fairly eccentric elliptic orbit

y(0) = 0 'Start with satellite crossing x-axis

UElip = -alpha / x(0) 'Elliptic orbit initial potential energy

TElip = Ecirc - UElip 'Elliptic orbit initial kinetic energy

v = c * Sqr(1 - (m ^ 2 * c ^ 4 / (TElip + m * c ^ 2) ^ 2)) 'Initial speed

D = x(0) 'Initial distance from central body to satellite

vx = 0

vy = v

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

Fx = -alpha / D ^ 2

Fy = 0

ax = Fx / (gamma * m)

ay = 0

vxhalf = vx + ax * deltat / 2 'vx at t=deltat/2

vyhalf = vy + ay * deltat / 2

vmax = 0

For i = 1 To Orbits * steps - 1

x(i) = x(i - 1) + vxhalf * deltat

y(i) = y(i - 1) + vyhalf * deltat

D = Sqr(x(i) ^ 2 + y(i) ^ 2)

vx = vxhalf + ax * deltat / 2

vy = vyhalf + ay * deltat / 2

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

If v > vmax Then vmax = v

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

F = alpha / D ^ 2

Fx = -F * x(i) / D

Fy = -F * y(i) / D

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

Num = (c ^ 2 + gamma ^ 2 * vy ^ 2) * (c ^ 2 * gamma * m + gamma ^ 3 * m * vx ^ 2) - gamma ^ 5 * m * vx ^ 2 * vy ^ 2

ax = Denom / Num

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

Num = (c ^ 2 + gamma ^ 2 * vx ^ 2) * (c ^ 2 * gamma * m + gamma ^ 3 * m * vy ^ 2) - gamma ^ 5 * m * vx ^ 2 * vy ^ 2

ay = Denom / Num

vxhalf = vxhalf + ax * deltat

vyhalf = vyhalf + ay * deltat

Next i

MsgBox ("vmax/c = " & vmax / c)

'Output x and y coordinates for plotting.

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

For i = 0 To steps - 1

Write #1, x(Orbits * i), y(Orbits * i)

Next i

Close

MsgBox ("Ready for plotting")

End Sub