The goal of this lesson is to find the orbit elements from x''' y''' z'''
for two times and the ellapsed time in days. 
given

k=.01720209895    ie=23.441028
 
x31 =-0.266166523942
y31 = 0.4618979603
z31 =-0.138236049035

x33 =-0.265112136054
y33 = 0.399047678984
z33 =-0.346297286144

t3-t1=10 days

Rotate yz through ie.
y21 = y31*Cos[ie] - z31*Cos[ie] =  0.368786535984  
z21 = y31*Cos[ie] + z31*Cos[ie] = -0.310572731996

y23 = y33*Cos[ie] - z33*Cos[ie] =  0.228355494265
z23 = y33*Cos[ie] + z33*Cos[ie] = -0.476460521693

l1 = [x31, y21, z21]
l3 = [x33, y23, z23]

i1 Cross i3 =a1*ii + b1*jj + c1*kk =

| ii   jj   kk  |
| x31, y21, z21 |
| x33, y23, z23 |

a1 = -0.104791235608
b1 = -0.044481240475
c1 =  0.0369891981709

Tan[i] = Sqrt[a1^2 + b1^2]/c1
If c1<0 then i=i+180
i = inclination
i = 72 degrees

Tan[o] = a/-b
If -b<0 then o=o+180
o= Longitude of the Ascending Node
o=293 degrees

Find yb
yb is (area of segment)/(area of triangle) = almost 1

2f = th3-th1, radians
r1 = Sqrt[x31^2 + y31^2 + z31^2]
r3 = Sqrt[x33^2 + y33^2 + z33^2]

Cos[2*f] = (x31*x33 + y31*y33 + z31*z33)/(r1*r3)
f=0.188253874989 radians
m = k^2 * (t3-t1)^2/(2 * Sqrt[r1*r3] * Cos[f])^3

l = (r1 + r3)/(2 * Sqrt[r1*r3] * Cos[f]) - .5

Solve this equation for yb

yb^2  = m/(l + (Sin[g/2])^2)
yb^3-yb^2 = m*(2*g-Sin[2*g])/(Sin[g])^3

HP method

E14
<< RAD CLEAR { << M YB x^2 / L - Sqrt ASIN 2 * 'G' STO G 2 * G 2 *
SIN - M * G SIN 3 y^x / YB x^2 + YB 3 y^x - RE  >> } {YB} {1} ROOT >>

Find P Q Th1 e E T1 W
P = (( r1 * r3 * Sin[2 * f] * YB )/(k^2 * (T2-T1)^2))^2
P = Semilatus Rectum = .5106
q1 = P/r1 - 1 = e * Cos[th1]
q3 = P/r3 - 1
q2 = (q1 * Cos[2* f]-q3) / Sin[2 *f] = e * Sin[th1]
e = Sqrt[ q1^2+q2^2]
e = eccentricity = .2
th1 = True Anomaly time 1
Change to degree mode
th1= ArcTan[q2/q1]
If q1<0 Add 180 degrees to th1
th1= 111.366647133 drgrees

q = p/(1+e)= perihelion distance
a = q/(1-e) = semimajor axis
q = .4255  a= .531875

se=Sin[E] = ( r1 * Sin[th1] ) / (a * Sqrt[1-e^2] )
ce=Cos[E] = (r1 * Cos[th1])/a + e

E = ArcTan[se/ce] in radians
If ce<0 add Pi

E1 = Eccentric Anomaly time 1 = 1.748988271 radians

Mean Anomaly M = k * t1/a^1.5 = E1 - e * Sin[E]
t1 = 35 days

Find W = argument of perigee = angle in plane of the orbit between 
perihelion an ascending node = puncture point

i = 72 
o = 293

x33 = -0.265112136055
y23 =  0.228355494265
z23 = -0.476460521693

x1 = x33 * Cos[o] + y23 * Sin[o]
y1 = y23 * Cos[o] - x33 * Sin[o]

x1 = -.313789904881
y1 = -.154811407975
y = y1/Cos[i] 
y = -.500980239901

y/x1=1.59654670883 = tan[u3]
If x1<0, add 180 degrees

u3 = 237.938952042
w = u3-th3 = 105 degrees