Option Base 1
Sub Path()
' program Path
'----------------------------------------------------------------------
'
' Author: Morris G. Anderson
'
'
' See the book "Time, Matter, and Gravity" copyright 2004 by Morris G. Anderson
' for theory behind the equations in this program. This book along with
' supporting information can be found at http://anderson.morris.home.att.net/
'
' This program can be used to calculate the path of an orbit,
' the precession of an orbit, and the path of a trajectory.
'
' Created on August 28, 2002
' Last update August 25, 2006 by Morris G. Anderson
'
' This program was originally written to run on an IBM Type 9112-265 workstation
' compiled as follows:
'
' xlf -qautodbl=dblpad -o path path_aix_march_30_2004.f
'
' The original code is documented in the book Time, Matter, and Gravity.
'
'----------------------------------------------------------------------
' Converted to Excel Visual Basic by Z. Daniel Judd 8/16/04
'
' Further updates have been made by Morris G. Anderson to improve the
' accuracy. However, the fundamental equations are the same as documented
' in the book Time, Matter, and Gravity.
'
'----------------------------------------------------------------------
' Definitions
'
' a = angle of change for the vector rc per iteration
' am = maximum angle for calculation
' b = motion vector. 1 = i, 2 = j, 3 = k, 4 = magnitude of v / c
' c = speed of light at a given point
' cr = the reference speed of light at a distance of rr from m
' cs = standard speed of light
' g = gravitational constant
' m = mass of gravitational body
' ni = number of increments for path calculation
' Pii = 3.14159265359
' r = position vector. 1 = i, 2 = j, 3 = k, 4 = magnitude of r
' rc = curvature vector. 1 = i, 2 = j, 3 = k, 4 = magnitude of rc
' rr = radius from center of m corresponding to cr
' sm = maximum path (length / positon radius) per calculation, if = 0.0 then
' there is no constraint
' v = velocity
' vt = temporary vector used in calculations
' w = rotation vector
' x = x Cartesian coordinate
' y = y Cartesian coordinate
' z = Z Cartesian coordinate
'----------------------------------------------------------------------
Dim tol_1 As Double
tol_1 = 0.00000000000001
' Tol_2 is used to determine which method to use for minimizing
' numerical errors
Dim tol_2 As Double
tol_2 = 0.1
Dim a As Double, ab As Double, am As Double, cr As Double, cs As Double
Dim f0 As Double, f1 As Double, g As Double, m As Double, pii As Double
Dim rr As Double, s As Double, sm As Double, t As Double
Dim v As Double, vx As Double, vy As Double, vz As Double
Dim x0 As Double, x1 As Double, xj As Double
Dim c(260000) As Double, si(260000) As Double, W(4) As Double
Dim x(260000) As Double, y(260000) As Double, z(260000) As Double
Dim b(260000, 4) As Double, r(260000, 4) As Double, rc(260000, 4) As Double
Worksheets("xyz").Activate
With Cells
.Clear
.NumberFormat = "General"
.HorizontalAlignment = xlCenter
.Font.Name = "Arial"
.Font.Size = 10
End With
Worksheets("Input").Activate
Worksheets("Input").Range("icount").Value = ""
Worksheets("Input").Range("am_out").Value = ""
Worksheets("Input").Range("ab_out").Value = ""
Worksheets("Input").Range("curve_warning").Value = ""
Worksheets("Input").Range("t").Value = ""
Worksheets("Input").Range("s").Value = ""
Worksheets("Input").Range("x_end").Value = ""
Worksheets("Input").Range("y_end").Value = ""
Worksheets("Input").Range("z_end").Value = ""
Worksheets("Input").Range("r_max").Value = ""
Worksheets("Input").Range("r_min").Value = ""
Worksheets("Input").Range("beta_max").Value = ""
Worksheets("Input").Range("beta_min").Value = ""
Worksheets("Input").Range("v_max").Value = ""
Worksheets("Input").Range("v_min").Value = ""
Worksheets("Input").Range("c_max").Value = ""
Worksheets("Input").Range("c_min").Value = ""
Worksheets("Input").Range("Iterating").Value = "Iterating"
' Read input file data
' Read in program option. 1 = Orbital Precession, 2 = Trajectory
Opt1 = Worksheets("Input").OptionButton1.Value
If Opt1 = True Then
option1 = 1
Else
option1 = 2
End If
' Read in the standard speed of light, cs, in m/s
cs = Worksheets("Input").Range("cs").Value
' Debug.Print "cs="; cs
pii = Application.pi()
' Read in the reference speed of light in m/s
cr = Worksheets("Input").Range("cr").Value
' Read in the radius from m for cr in meters.
rr = Worksheets("Input").Range("rr").Value
' Read in g, universal gravitational constant. Units = m^3/(kg s)
g = Worksheets("Input").Range("g").Value
' Read in m, This is the mass of the governing body in kilograms
m = Worksheets("Input").Range("m").Value
' Read in x, y, z, Initial location coordinates with respect to
' the center of mass of the governing body, m, in meters.
x(1) = Worksheets("Input").Range("x").Value
y(1) = Worksheets("Input").Range("y").Value
z(1) = Worksheets("Input").Range("z").Value
' Read in vx, vy, vz, Initial velocity in m/s. If Beta_Check = True,
' then these values are treated as vx/c, vy/c, vz/c
Beta_Check = Worksheets("Input").CheckBox1.Value
vx = Worksheets("Input").Range("vx").Value
vy = Worksheets("Input").Range("vy").Value
vz = Worksheets("Input").Range("vz").Value
' Read in number of rotation steps (ni), maximum rotation in
' degrees (am), and the maximum path (length / position radius)
' of rotation (sm). If sm is set to 0.0, then no limit applied
' to the path length of each step.
ni = Worksheets("Input").Range("ni").Value
am = Worksheets("Input").Range("am").Value
sm = Worksheets("Input").Range("sm").Value
am = am * pii / 180#
a = am / ni
'----------------------------------------------------------------------
' Calculate initial values for starting point
'----------------------------------------------------------------------
r(1, 4) = Sqr(x(1) ^ 2 + y(1) ^ 2 + z(1) ^ 2)
r(1, 1) = x(1) / r(1, 4)
r(1, 2) = y(1) / r(1, 4)
r(1, 3) = z(1) / r(1, 4)
'
' For the following calculation:
' See equation 4-10 on page 30 of "Time, Matter, and Gravity"
'
c(1) = cr * Exp((2# * m * g / cs / cs) * ((r(1, 4) - rr) / (rr * r(1, 4))))
If Beta_Check = True Then
b(1, 4) = Sqr(vx ^ 2 + vy ^ 2 + vz ^ 2)
If b(1, 4) = 0 Then
b(1, 1) = -r(1, 1)
b(1, 2) = -r(1, 2)
b(1, 3) = -r(1, 3)
Else
b(1, 1) = vx / b(1, 4)
b(1, 2) = vy / b(1, 4)
b(1, 3) = vz / b(1, 4)
End If
Else
v = Sqr(vx * vx + vy * vy + vz * vz)
b(1, 4) = v / c(1)
If b(1, 4) = 0 Then
b(1, 1) = -r(1, 1)
b(1, 2) = -r(1, 2)
b(1, 3) = -r(1, 3)
Else
b(1, 1) = (vx / c(1)) / b(1, 4)
b(1, 2) = (vy / c(1)) / b(1, 4)
b(1, 3) = (vz / c(1)) / b(1, 4)
End If
End If
If b(1, 4) = 0 And sm = 0 Then
sm = 0.000001
End If
' ---------------------------------------------------------------------
' Calculate direction of rotation vector - use cross product. This
' is a unit vector that is perpendicular to the plane of rotation.
' ---------------------------------------------------------------------
If b(1, 4) = 0 Then
W(1) = 1
W(2) = 0
W(3) = 0
W(4) = 1
Else
W(1) = r(1, 2) * b(1, 3) - r(1, 3) * b(1, 2)
W(2) = r(1, 3) * b(1, 1) - r(1, 1) * b(1, 3)
W(3) = r(1, 1) * b(1, 2) - r(1, 2) * b(1, 1)
If W(1) = 0 And W(2) = 0 And W(3) = 0 Then
W(1) = 1
W(4) = 1
Else
W(4) = Sqr(W(1) ^ 2# + W(2) ^ 2# + W(3) ^ 2#)
W(1) = W(1) / W(4)
W(2) = W(2) / W(4)
W(3) = W(3) / W(4)
W(4) = 1#
End If
End If
' ---------------------------------------------------------------------
' Calculate the direction of the curvature vector at the starting
' point. This is done with a cross product.
' ---------------------------------------------------------------------
rc(1, 1) = b(1, 2) * W(3) - b(1, 3) * W(2)
rc(1, 2) = b(1, 3) * W(1) - b(1, 1) * W(3)
rc(1, 3) = b(1, 1) * W(2) - b(1, 2) * W(1)
' ---------------------------------------------------------------------
' option1 = 1. The program calculates the precession of an orbit.
' If the max angle (am) is input as 360 this is done
' for one revolution. If it is put in as 2 * 360 this
' is done for two revolutions, and so forth.
' ---------------------------------------------------------------------
If option1 <> 1 Then GoTo 100
nmax = 50
icount = 1
x0 = am
f0 = curve(a, b, c, cs, g, m, ni, r, rc, si, sm, tol_1, tol_2, W, x, y, z)
am = am + 0.01
a = am / ni
x1 = am
f1 = curve(a, b, c, cs, g, m, ni, r, rc, si, sm, tol_1, tol_2, W, x, y, z)
For j = 1 To nmax
xj = x1 - f1 * (x1 - x0) / (f1 - f0)
If Abs((xj - x1) / x1) < tol_1 Then GoTo 20
x0 = x1
x1 = xj
f0 = f1
am = x1
a = am / ni
f1 = curve(a, b, c, cs, g, m, ni, r, rc, si, sm, tol_1, tol_2, W, x, y, z)
icount = icount + 1
'Write ( 6,'(''iteration'', i5)') icount
Worksheets("Input").Range("icount").Value = icount
Next j
10 'Continue
20 'Continue
If icount >= nmax Then
Worksheets("Input").Range("warning").Value = "Solution Failed to Converge for Option 1"
End If
ab = b(ni + 1, 1) * b(1, 1) + b(ni + 1, 2) * b(1, 2) + b(ni + 1, 3) * b(1, 3)
ab = Application.Acos(ab) * 180# / pii
Worksheets("Input").Range("am_out").Value = am * 180# / pii
Worksheets("Input").Range("ab_out").Value = ab
' This next command removes the "Iterating" word from the spread sheet
Worksheets("Input").Range("Iterating").Value = ""
'
' Calculate the total rotation time and path length
' Note: Point 1 is the starting location, therefore, si(1) = 0
t = 0#
s = 0#
jrmin = 1
jrmax = 1
For j = 1 To ni
' Calculate the time required from start to finish.
t = t + si(j + 1) / ((b(j, 4) * c(j) + b(j + 1, 4) * c(j + 1)) / 2#)
s = s + si(j + 1)
' Check for min and max radius
If r(j + 1, 4) < r(jrmin, 4) Then jrmin = j + 1
If r(j + 1, 4) > r(jrmax, 4) Then jrmax = j + 1
Next j
GoTo 200
100 'Continue
' ---------------------------------------------------------------------
' option1 = 2. This option calculates the trajectory of motion.
' For example: if the max angle is input as 35 this
' is done for a total change in the path direction of
' 35 degrees.
' ---------------------------------------------------------------------
If option1 <> 2 Then GoTo 200
t = 0#
s = 0#
jrmin = 1
jrmax = 1
f1 = curve(a, b, c, cs, g, m, ni, r, rc, si, sm, tol_1, tol_2, W, x, y, z)
Worksheets("xyz").Activate
Worksheets("xyz").Range("A1").Select
With Selection
.Range("A1").Value = "x"
.Range("B1").Value = "y"
.Range("C1").Value = "z"
.Range("D1").Value = "rc"
.Range("E1").Value = "c"
.Range("F1").Value = "b"
.Range("G1").Value = "t"
End With
With Range("A1:G1")
.Font.Bold = True
.Borders(xlEdgeBottom).LineStyle = xlContinuous
End With
j = 1
Worksheets("xyz").Range("A1").Select
With Selection
.Offset(j, 0).Range("A1").Value = x(j)
.Offset(j, 1).Range("A1").Value = y(j)
.Offset(j, 2).Range("A1").Value = z(j)
.Offset(j, 3).Range("A1").Value = rc(j, 4)
.Offset(j, 4).Range("A1").Value = c(j)
.Offset(j, 5).Range("A1").Value = b(j, 4)
.Offset(j, 6).Range("A1").Value = t
End With
For j = 1 To ni
' Calculate the time required from start to finish.
' Note: si(1) = 0 as point 1 is the starting location
t = t + si(j + 1) / ((b(j, 4) * c(j) + b(j + 1, 4) * c(j + 1)) / 2#)
s = s + si(j + 1)
' Check for min and max radius
If r(j + 1, 4) < r(jrmin, 4) Then jrmin = j + 1
If r(j + 1, 4) > r(jrmax, 4) Then jrmax = j + 1
' Write out information at each step of the trajectory
Worksheets("xyz").Range("A1").Select
With Selection
.Offset(j + 1, 0).Range("A1").Value = x(j + 1)
.Offset(j + 1, 1).Range("A1").Value = y(j + 1)
.Offset(j + 1, 2).Range("A1").Value = z(j + 1)
.Offset(j + 1, 3).Range("A1").Value = rc(j + 1, 4)
.Offset(j + 1, 4).Range("A1").Value = c(j + 1)
.Offset(j + 1, 5).Range("A1").Value = b(j + 1, 4)
.Offset(j + 1, 6).Range("A1").Value = t
End With
Next j
' Calculate the total angle of rotation
ab = b(ni + 1, 1) * b(1, 1) + b(ni + 1, 2) * b(1, 2) + b(ni + 1, 3) * b(1, 3)
ab = Application.Acos(ab) * 180# / pii
Worksheets("Input").Range("ab_out").Value = ab
200 'Continue
Worksheets("Input").Range("t").Value = t
Worksheets("Input").Range("s").Value = s
Worksheets("Input").Range("x_end").Value = x(ni + 1)
Worksheets("Input").Range("y_end").Value = y(ni + 1)
Worksheets("Input").Range("z_end").Value = z(ni + 1)
Worksheets("Input").Range("r_max").Value = r(jrmax, 4)
Worksheets("Input").Range("beta_max").Value = b(jrmax, 4)
Worksheets("Input").Range("v_max").Value = b(jrmax, 4) * c(jrmax)
Worksheets("Input").Range("c_max").Value = c(jrmax)
Worksheets("Input").Range("r_min").Value = r(jrmin, 4)
Worksheets("Input").Range("beta_min").Value = b(jrmin, 4)
Worksheets("Input").Range("v_min").Value = b(jrmin, 4) * c(jrmin)
Worksheets("Input").Range("c_min").Value = c(jrmin)
' This next command removes the "Iterating" word from the spread sheet
Worksheets("Input").Range("Iterating").Value = ""
Worksheets("Input").Activate
End Sub
'----------------------------------------------------------------------
'
' Function curve
'
' This function calculates the average radius of curvature for a
' given amount of orbital rotation and the end point conditions.
'----------------------------------------------------------------------
Function curve(a As Double, b() As Double, c() As Double, cs As Double, _
g As Double, m As Double, ni, r() As Double, rc() As Double, _
si() As Double, sm As Double, tol_1 As Double, tol_2 As Double, _
W() As Double, x() As Double, y() As Double, _
z() As Double) As Double
Dim ai As Double
Dim g1 As Double, g2 As Double, g3 As Double, g4 As Double, g5 As Double
Dim g6 As Double, g7 As Double, g8 As Double
Dim rca As Double, rcheck As Double
Dim vt(4) As Double
si(1) = 0#
For n = 1 To ni
' Calculate path radius of curvature at point n
'
' For the following calculation:
' See equation 7-12 on page 72 of "Time, Matter, and Gravity"
'
'
' Check to see if you have gone more than 1000 bilion light years,
' If so stop. You have gone far enough
If r(n, 4) > 1E+28 Then
ni = n - 1
GoTo 40
End If
'
' Check for maximum radius of curvature. If greater than limiting value assume
' straight line motion.
If (Abs((r(n, 1) * rc(n, 1) + r(n, 2) * rc(n, 2) _
+ r(n, 3) * rc(n, 3))) * 2# * m * g / (r(n, 4) * r(n, 4) * cs ^ 2#) > 1E-99) Then
rc(n, 4) = (r(n, 4) * r(n, 4) * cs ^ 2#) / _
((1# + 1# / b(n, 4) ^ 2#) * m * g) / _
(r(n, 1) * rc(n, 1) + r(n, 2) * rc(n, 2) + r(n, 3) * rc(n, 3))
rca = rc(n, 4)
'
' If rca is less than zero, stop calculation. For a single body
' system, numerical errors can cause rca to become less than zero.
' However, this is incorrect and indicates that motion should actually
' be directed to the center of the governing body with an infinite
' radius of curvature.
'
If rca < 0# Then
ni = n - 1
GoTo 40
End If
' Calculate direction of curvature vector for n+1
ai = a
' Check on limit of path length for iteration
If sm > 0# Then
If (rca * ai) > (sm * r(n, 4)) Then
ai = sm * r(n, 4) / rca
End If
End If
vt(1) = rc(n, 1) + b(n, 1) * Tan(ai)
vt(2) = rc(n, 2) + b(n, 2) * Tan(ai)
vt(3) = rc(n, 3) + b(n, 3) * Tan(ai)
vt(4) = Sqr(vt(1) * vt(1) + vt(2) * vt(2) + vt(3) * vt(3))
rc(n + 1, 1) = vt(1) / vt(4)
rc(n + 1, 2) = vt(2) / vt(4)
rc(n + 1, 3) = vt(3) / vt(4)
b(n + 1, 1) = W(2) * rc(n + 1, 3) - W(3) * rc(n + 1, 2)
b(n + 1, 2) = W(3) * rc(n + 1, 1) - W(1) * rc(n + 1, 3)
b(n + 1, 3) = W(1) * rc(n + 1, 2) - W(2) * rc(n + 1, 1)
' x(n + 1) = x(n) + rca * (rc(n + 1, 1) - rc(n, 1))
' y(n + 1) = y(n) + rca * (rc(n + 1, 2) - rc(n, 2))
' z(n + 1) = z(n) + rca * (rc(n + 1, 3) - rc(n, 3))
' The following method reduces numerical truncation errors over
' the above method that has been commented out.
vt(1) = b(n, 1) + b(n + 1, 1)
vt(2) = b(n, 2) + b(n + 1, 2)
vt(3) = b(n, 3) + b(n + 1, 3)
vt(4) = Sqr(vt(1) * vt(1) + vt(2) * vt(2) + vt(3) * vt(3))
vt(1) = vt(1) / vt(4)
vt(2) = vt(2) / vt(4)
vt(3) = vt(3) / vt(4)
vt(4) = 2 * rca * Sin(ai / 2)
x(n + 1) = x(n) + vt(1) * vt(4)
y(n + 1) = y(n) + vt(2) * vt(4)
z(n + 1) = z(n) + vt(3) * vt(4)
Else
' Set Rc = 1E+99 as a flag to the user but calculate position of next point
' based on straight line motion
rc(n, 4) = 1E+99
rca = rc(n, 4)
'
' For straight line motion the direction of the radius of curvature does not change
'
rc(n + 1, 1) = rc(n, 1)
rc(n + 1, 2) = rc(n, 2)
rc(n + 1, 3) = rc(n, 3)
b(n + 1, 1) = b(n, 1)
b(n + 1, 2) = b(n, 2)
b(n + 1, 3) = b(n, 3)
x(n + 1) = x(n) + b(n, 1) * sm * r(n, 4)
y(n + 1) = y(n) + b(n, 2) * sm * r(n, 4)
z(n + 1) = z(n) + b(n, 3) * sm * r(n, 4)
End If
r(n + 1, 4) = Sqr((x(n + 1)) ^ 2 + (y(n + 1)) ^ 2 + (z(n + 1)) ^ 2)
r(n + 1, 1) = x(n + 1) / r(n + 1, 4)
r(n + 1, 2) = y(n + 1) / r(n + 1, 4)
r(n + 1, 3) = z(n + 1) / r(n + 1, 4)
' -----------------------------------
' Calculate conditions for point n+1.
' -----------------------------------
' Calculate speed of light at point n+1.
'
' For the following calculation:
' See equation 4-10 and 4-14 on page 30 and 31 of "Time, Matter, and Gravity"
'
g1 = ((2# * m * g / cs / cs) * ((r(n + 1, 4) - r(1, 4)) / (r(1, 4) * r(n + 1, 4))))
' Choose the best method to minimize numerical errors
If (Abs(g1) > tol_2) Then
c(n + 1) = c(1) * Exp((2# * m * g / cs / cs) _
* ((r(n + 1, 4) - r(1, 4)) / (r(1, 4) * r(n + 1, 4))))
dc = c(1) - c(n + 1)
Else
g2 = g1 * g1
g3 = g2 * g1
g4 = g3 * g1
g5 = g4 * g1
g6 = g5 * g1
g7 = g6 * g1
g8 = g7 * g1
dc = c(1) * (-g1 - g2 / 2# - g3 / 6# - g4 / 24# - g5 / 120# _
- g6 / 720# - g7 / 5040# - g8 / 40320#)
c(n + 1) = c(1) - dc
End If
' Calculate magnitude of beta at point n+1
b(n + 1, 4) = Sqr((dc + b(1, 4) * b(1, 4) * c(n + 1)) / c(1))
j = 0
10 'Continue
' Jump out of loop if straight line motion
If rc(n, 4) = 1E+99 Then GoTo 20
If rc(n + 1, 4) = 1E+99 Then GoTo 20
j = j + 1
If j > 100 Then
Worksheets("Input").Range("curve_warning").Value = "Solution Failed to Converge for rc at n = " & n & " try reducing sm or increasing ni"
GoTo 20
End If
' Calculate path radius of curvature at point n+1
'
' For the following calculation:
' See equation 7-12 on page 72 of "Time, Matter, and Gravity"
'
' Check for minimum dot product. If less than limiting value assume
' straight line motion.
If (Abs((r(n + 1, 1) * rc(n + 1, 1) + r(n + 1, 2) * rc(n + 1, 2) _
+ r(n + 1, 3) * rc(n + 1, 3))) * 2# * m * g / _
(r(n + 1, 4) * r(n + 1, 4) * cs ^ 2#) > 1E-99) Then
rc(n + 1, 4) = (r(n + 1, 4) * r(n + 1, 4) * cs ^ 2#) / _
((1# + 1# / b(n + 1, 4) ^ 2#) * m * g) / _
(r(n + 1, 1) * rc(n + 1, 1) + r(n + 1, 2) * rc(n + 1, 2) _
+ r(n + 1, 3) * rc(n + 1, 3))
' Check on difference in curvature used to calculate point n+1. If
' greater than tolerance, calculate point n+1 based on the average
' radius of curvature between point n and n+1.
rcheck = (rc(n, 4) + rc(n + 1, 4)) / 2#
If Abs(1# - (rca / rcheck)) < tol_1 Then
GoTo 20
Else
rca = rcheck
End If
' Check on limit of path length for iteration
If sm > 0# Then
If (rca * ai) > (sm * r(n, 4)) Then
ai = sm * r(n, 4) / rca
End If
End If
vt(1) = rc(n, 1) + b(n, 1) * Tan(ai)
vt(2) = rc(n, 2) + b(n, 2) * Tan(ai)
vt(3) = rc(n, 3) + b(n, 3) * Tan(ai)
vt(4) = Sqr(vt(1) * vt(1) + vt(2) * vt(2) + vt(3) * vt(3))
rc(n + 1, 1) = vt(1) / vt(4)
rc(n + 1, 2) = vt(2) / vt(4)
rc(n + 1, 3) = vt(3) / vt(4)
b(n + 1, 1) = W(2) * rc(n + 1, 3) - W(3) * rc(n + 1, 2)
b(n + 1, 2) = W(3) * rc(n + 1, 1) - W(1) * rc(n + 1, 3)
b(n + 1, 3) = W(1) * rc(n + 1, 2) - W(2) * rc(n + 1, 1)
' x(n + 1) = x(n) + rca * (rc(n + 1, 1) - rc(n, 1))
' y(n + 1) = y(n) + rca * (rc(n + 1, 2) - rc(n, 2))
' z(n + 1) = z(n) + rca * (rc(n + 1, 3) - rc(n, 3))
' The following method reduces numerical truncation errors over
' the above method that has been commented out.
vt(1) = b(n, 1) + b(n + 1, 1)
vt(2) = b(n, 2) + b(n + 1, 2)
vt(3) = b(n, 3) + b(n + 1, 3)
vt(4) = Sqr(vt(1) * vt(1) + vt(2) * vt(2) + vt(3) * vt(3))
vt(1) = vt(1) / vt(4)
vt(2) = vt(2) / vt(4)
vt(3) = vt(3) / vt(4)
vt(4) = 2 * rca * Sin(ai / 2)
x(n + 1) = x(n) + vt(1) * vt(4)
y(n + 1) = y(n) + vt(2) * vt(4)
z(n + 1) = z(n) + vt(3) * vt(4)
Else
' Set Rc = 1E+99 as a flag to the user but calculate position of next point
' based on straight line motion
rc(n + 1, 4) = 1E+99
rc(n + 1, 1) = rc(n, 1)
rc(n + 1, 2) = rc(n, 2)
rc(n + 1, 3) = rc(n, 3)
b(n + 1, 1) = b(n, 1)
b(n + 1, 2) = b(n, 2)
b(n + 1, 3) = b(n, 3)
x(n + 1) = x(n) + b(n, 1) * sm * r(n, 4)
y(n + 1) = y(n) + b(n, 2) * sm * r(n, 4)
z(n + 1) = z(n) + b(n, 3) * sm * r(n, 4)
End If
r(n + 1, 4) = Sqr((x(n + 1)) ^ 2 + (y(n + 1)) ^ 2 + (z(n + 1)) ^ 2)
r(n + 1, 1) = x(n + 1) / r(n + 1, 4)
r(n + 1, 2) = y(n + 1) / r(n + 1, 4)
r(n + 1, 3) = z(n + 1) / r(n + 1, 4)
' -----------------------------------
' Calculate conditions for point n+1.
' -----------------------------------
' Calculate speed of light at point n+1.
'
' For the following calculation:
' See equation 4-10 and 4-14 on page 30 and 31 of "Time, Matter, and Gravity"
g1 = ((2# * m * g / cs / cs) * ((r(n + 1, 4) - r(1, 4)) / (r(1, 4) * r(n + 1, 4))))
' Choose the best method to minimize numerical errors
If (Abs(g1) > tol_2) Then
c(n + 1) = c(1) * Exp((2# * m * g / cs / cs) _
* ((r(n + 1, 4) - r(1, 4)) / (r(1, 4) * r(n + 1, 4))))
dc = c(1) - c(n + 1)
Else
g2 = g1 * g1
g3 = g2 * g1
g4 = g3 * g1
g5 = g4 * g1
g6 = g5 * g1
g7 = g6 * g1
g8 = g7 * g1
dc = c(1) * (-g1 - g2 / 2# - g3 / 6# - g4 / 24# - g5 / 120# _
- g6 / 720# - g7 / 5040# - g8 / 40320#)
c(n + 1) = c(1) - dc
End If
' Calculate magnitude of beta at point n+1
'
' For the following calculation:
' See equation 4-18 on page 32 of "Time, Matter, and Gravity"
b(n + 1, 4) = Sqr((dc + b(1, 4) * b(1, 4) * c(n + 1)) / c(1))
GoTo 10
20 'Continue
' Calculate path length for iteration
If rca = 1E+99 Or rc(n, 4) = 1E+99 Or rc(n + 1, 4) = 1E+99 Then
' Assume straight line motion
si(n + 1) = Sqr((x(n + 1) - x(n)) ^ 2# _
+ (y(n + 1) - y(n)) ^ 2# _
+ (z(n + 1) - z(n)) ^ 2#)
Else
si(n + 1) = rcheck * ai
End If
' Calculate dot product between position unit vector and direction
' unit vector at point n+1.
curve = b(n + 1, 1) * r(n + 1, 1) + b(n + 1, 2) * r(n + 1, 2) + b(n + 1, 3) * r(n + 1, 3)
30 'Continue
Next n
40 'Continue
End Function
Copyright:
These documents cannot be distributed without the written permission of the author. This permission will be normally granted, upon request.
Time, Matter, and Gravity home page