In order to control the excution of the program following control parameters are used:
Test  Ecc 1  l  m  C  h  vr  vu  k  eps  Sl R 
1.1  .7  100  2  .007  14.1421  .09899  .00699  .994534  .060  99.99 
1.2  .5  100  2  .005  14.1421  .07071  .00500  .994534  .060  99.99 
1.3  .3  100  2  .003  14.1421  .04242  .00300  .994534  .060  99.99 
1.4  .7  150  2  .007  17.3205  .09899  .00698  .994534  .040  149.99 
1.5  .7  200  2  .007  20.0000  .08082  .00466  .994534  .030  199.99 
1.6  .7  100  1.5  .007  12.2474  .08573  .00699  .994534  .044  100.00 
1.7  .7  100  3  .007  17.3205  .12124  .00699  .994534  .090  99.99 
d^2u mu  + U =  d phi^2 h^2 Equation 15.10 


For i = 0 To imax For j = 1 To di Select Case test Case Is = 2, 3 ' 15.25 page 196 ac = m / h2 + C3 * m * u * u  u vu = vu + ac * ddPhi u = u + vu * ddPhi r = 1 / u phi = Phi00 * pi / 180 + i * dPhi + j * ddPhi TESTMAXR eccr If t2 > 0 Then t2 = 0 End Select next j ' Display xxx = r * Cos(phi) yyy = r * Sin(phi) px = px0 + fac * xxx py = py0  fac * yyy Form2.PSet (px, py), QBColor(kleur) next i 

For i = 0 To imax For j = 1 To di Select Case test Case Is = 4 ' 15.24 page 196 C2 = 2 k1 = (k2  1) / h2 ' Eq 15.24 u2 = u * u fac3 = k1 + u * (2 * m / h2  u + C2 * m * u2) If fac3 < 0 And Count = 0 Then sign1 = sign1 * 1 Count = 1 End If If sign1 = 1 Then kleur = LGreen Else kleur = Yellow If Count > 0 Then Count = Count + 1: If Count >= maxcount Then Count = 0: If sign1 = 1 Then kleur = LRed Else kleur = Black End If vu = sign1 * Sqr(Abs(fac3)) u = u + vu * ddPhi r = 1 / u phi = Phi00 * pi / 180 + i * dPhi + j * ddPhi TESTMAXR eccr If t2 > 0 Then t2 = 0 



Test 5 is a combination of the equations 15.21, 15.22 and 15.23. Equation 15.23 is the most difficult. This equation contains the two derivatives dr/dt^2 and dhi/dt^2. dhi/dt is calculated in equation 15.21. The derivative dr/dt^2 is called fac2. That means that dr/dt = sqr(fac2) The problem is that fac2 can become negative, which has to be taken into account. 
Test  Ecc  l  m  C  h  vr  vu  k  Epsilon  Precession 
5.1  .7  100  2  .007  14.1421  .09899  .00699  .994534  .060  25.65 
5.2  .6  100  2  .006  14.1421  .08485  .00600  .99318  .060  25.65 
5.3  .5  100  2  .005  14.1421  .07071  .00500  .992045  .059  25.65 
5.4  .4  100  2  .004  14.1421  .05656  .00400  .991112  .059  25.65 
5.5  .3  100  2  .003  14.1421  .04242  .00300  .990412  .060  25.65 
5.6  .7  100  3  .007  17.3205  .12124  .00699  .99145  .090  43.80 
5.7  .7  100  4  .007  20.0000  .14000  .00699  .98802  .120  67.41 
5.8  .7  100  1  .007  10.0000  .07000  .00699  .997301  .030  11.76 
m^2 Eps = 3 *  h^2 2 d u m h^2* U^2  + U =  + Eps *  d phi^2 h^2 m Equation 15.27 2 d u m h^2* U^2  =  + Eps *   u d phi^2 h^2 m Equation 15.27 



mu U =  + C * cos(Phi  Phi0) h^2 Equation 15.11 

l  = 1 + Ecc * (cos (Phi  Phi0)) R Equation 15.12 


For i = 1 To imax phi = Phi0 + i * dPhi + pi / 2 Select Case subtest Case Is = 3 ' equation 15.31 u0 = m / h2 * (1 + ecc * Cos(phi)) u1 = aa + bb * phi * Sin(phi) + CC * Cos(2 * phi) u = u0 + Eps * u1 r = 1 / u End Select TESTMAXR eccr If t2 > 0 Then t2 = 0 ' Display px = px0 + fac * r * Cos(phi): py = py0  fac * r * Sin(phi) Form2.PSet (px, py), QBColor(kleur) If maxcount = irev Then GoTo end1 If stop1 = 1 Then GoTo end2 Next i 
m d^2u h^2 u0 =  * (1 + Ecc * cos(phi))  + u =  * u0^2 h^2 dt^2 m Equation 15.33 

m u0 =  * [1 + Ecc * cos(phi) + Eps*Ecc*phi*sin(phi)] h^2 Equation 15.34A 

m u0 =  * {1 + Ecc * cos[phi(1Eps)]} h^2 Equation 15.34B 


What makes numerical calculation difficult, specific for example equation 15.23 that the starting point should be to calculate the parameters and initial conditions using observations.
Specific you have to calculate the masses involved and the parameters h and k. The two parameters h and k are called constants, that may be true for each object (mass) involved, but in fact they are very important because they define the trajectory of each individual object. I expect you need as many parameter combinations (h,k) as objects involved
In the above I have calculated h directly from eccentricity but in general you cannot follow such an approach.
Observations in GR are specific difficult because proper time is used. GR does not use the concept of simultaneity, which makes the whole issue of observation "handling" difficult.
Back to my home page: Contents of This Document