Newtonian motion in 3D space, Question

[PB666]

This is one for all the computer and physics whizz kids out there.

I have been writing some code and came across what should be a standard problem, but I cannot find the answer and I have look all over the place. Wiki may have the answer but the symbology they are using is confused with other math symbols an it is difficult to follow. When you apply accelerations to a particle in 3D space you kick it from keplerian universe (rails, soft rails, whatever) to a newtonian universe, while it is fairly east to work the position, velocity and acceleration vectors is the newtonian state and also a bit harder to go from the newtonian state in the XY back to the Keplerian state. The three dimensional state that includes motion in the Z dimension however creates a problem.

The problem comes when applying motion along the Z axis. So I came up with the following solution (Edit: V is velocity, T is along the tangent to P, radial is along the parallel to r, r = <(P - O)_x,y,z>, O is the focal points position) Assuming the CB is large and the particle is small then r = ||P|| and points in <P> direction.

In order to calculate w which is V_T/r the relation ship of V_∙ to V = SQRT(V_T² + V_radial²) this is the same in 2 and 3 dimensions since velocities are in the polar coordinate system.

V_radial = P_u∙V where Pu denotes the normalized position vector whose magnitude is a moment of radius (both vectors have i, j and k components that are perpendicular to each other)

It seems like this is the correct answer but I have yet to find verification on the internet and it proves rather difficult to verify the conclusions through the programming since I need to create an artificial angle, Θ, that describes the objects position relative to the plane of motion. There are two positions when the vector jumps off the rails. At is apoapsis and at the SIGN(motion) * -Burnspan/2. Which leads to the second problem . . . where are the rails (notably in z).

Does anyone have a way to convert from a x, y, z=0 system back to a x, y, z system?

 

 

Edited December 29, 2017 by PB666

[mikegarrison]

It seems like you are trying to describe an orbit that is not in the ecliptic plane in a rectangular coordinate system that is based in the ecliptic plane, is that correct? Any particular reason why you have to use that particular coord system? Any coord system can be converted to any other, but that doesn't mean it should be.

Edited December 29, 2017 by mikegarrison

[PB666]

The opposite question is why not use R³ space completely, the reason for not using it completely is that during part of an orbit that is specifically non-inertial , and in general this is going to occur near the apoapsis, the newtonian coordinate translations introduce unnecessary inaccuracies and increase computational time, and since the calculations occur per interval of time, as e increases the error also increases. Therefore long periods where no thrust is applied its best to jump on rails. The reason is that to go off the rails only takes one more application of a subroutine no matter how much time has expired. Although the elliptical is less computation intensive it introduces more inaccuracies when if follows non-inertial time intervals (this is due to the sin/cos/tan functions). When thrusts are applied, and frankly its easier to factor space-time, thrust/m and any other acceleration that might need to be added (see other threads on rotation) in the R³ space than in Keplarian ellipticals.

The reason to jump  in and out of R3 space is particular to the problem of using ION drives, notable low energy or mass flow techniques for retaining Pe. In this circumstance thrust might not be applied all the time, or along the elliptical with radial vector of velocity (such radially down before Pe or radially up after Pe). Hence if you are on the elliptical and there are no forces applied on the axis normal to travel then there is no reason to rotate the elliptical. So that the elliptical can be used to calculate the Θ relative to the velocity vector for which create acceleration vectors (remember that we don't really care about force, but specific force, we assume that mass has been factored into the acceleration) and observe the result in order to tune the program. The problem is that at some point you may need to jump off the rails and see where you have been.

 

 

Public Sub Ellipticize(X As Double, Y As Double, Z As Double, vX As Double, vY As Double, vZ As Double, v As Double, SME As Double, a As Double, e As Double, l As Double, Pe As Double, Apo As Double, Optional theta as Double,  Optional phi As Double)
Dim vD As Double, vR As Double, vT As Double, Period As Double, Sweep As Double, Area As Double, XY As Double, posRatio As Double, velRatio As Double, DZ As Double
Dim Xu As Double, Yu As Double, Zu As Double, r As Double, vXu As Double, vYu As Double, vZu As Double, pvu As Double
Dim vRad As Double, vTan As Double, omega As Double, QueryAngle As Double
  a = -stdGP / (2 * SME)
  Period = (2 * Pi * a ^ 1.5) / uR
  Normalize X, Y, Z, r, Xu, Yu, Zu <===== unit vector definition
  vRad = DotProduct(Xu, Yu, Zu, vX, vY, vZ) <===== This is the Pu*V
  vTan = SQRT(v^2 - vRad^2) <====tangential velocity
  omega = vTan / r <==== angular velocity 
  Sweep = (omega * r ^ 2) / 2
  Area = Sweep * Period
  b = Area / (Pi * a)
  l = b ^ 2 / a
  e = (1 - b ^ 2 / a ^ 2) ^ 0.5
  Pe = a * (1 - e): Apo = a * (1 + e)
  If vRad = 0 Then <===If then else end if defines the keplerian theta. 
    If (r - 0.001) < Pe And Pe < (r + 0.001) Then Theta = 0 Else Theta = Pi <=== eliminate solutions with only one answer
   Else
    Theta = arccos((1 - l / r) / e) <===== get a preliminary solution 
    If vRad < 0 Then Theta = 2 * Pi - Theta <=====convert (reverse signs) of all solutions in which radial velocity is negative 
  End If
  ' need to add stuff here fore the cartesian relationals.
End Sub

 

Edited December 29, 2017 by PB666

[K^2]

No sane simulation is going to jump between Cartesian and Keplerian coordinate systems, because precision losses on each burst of your engines are going to make it suck. This is particularly bad when gravity is still a dominant force, and you are integrating it in Cartesian. But more importantly, the typical solution is actually a lot easier.

The reason Keplerian Elements show up is because they are constants of motion for a 1/R central potential. However, they aren't the canonical set. If you start with Hamiltonian, in place of semi-major axis and eccentricity, you'll end up with energy and angular momentum. This is super useful, because these two have very simple equations of motion under perturbation.

dE/dt = v·F

dL/dt = r⨯F

Energy gives you semi-major axis directly. Magnitude of angular momentum with respect to it gives you eccentricity. Orientation of angular momentum vector gives you inclination and ascending node. The only things missing are argument of periapsis and anomaly. These have equations of motion as well, but they are a pain to work with, and tend to be numerically unstable anyways.

Instead, we keep track of Cartesian location of the rocket, r. Given current true anomaly and periapsis, we compute expected position of the rocket and its velocity v. This already takes into account any influence of gravitational forces. We can now apply perturbation force F to compute v', r', E', and L' using your favorite integration method. We include velocity here for numerical stability. For this step, we treat E, L, and r as independent variables. Finally, using these new quantities, we recompute Keplerian Elements.

The chief advantage here is that we keep changes due to external forces separate from influence of gravity. That allows keeping track of a nice, clean orbit without nearly as many errors. The other advantage is that you simply don't have to worry about complex computations. The only things you have to do is convert between anomalies and compute position of the periapsis from these. That's just as easy to do in 3D as it is in 2D. You get inclination and argument of ascending node gratis.

[PB666]

I don't know about jumping but moment of force application in the Keplerian system was problematic. The force application problem in the R³ is easy like you say, very easy when using unit vectors, I just don't know how to rail it like the Keplerian. Although I must say I still have to make three passes over a time interval to settle at the right position and velocity when two forces are applied at once. But still its just a repetition replacing the old p1 with the new p1. Just to make sure P1 is equilibrated.

I can do this in R². The problem is RSS is that they have rigged the entire system to the tilt axis of the Earth, and although I want to test in the Earth elliptical,  eventually I want to model outside of Earth and model along exit inclinations that favor planetary intercepts and so I was lusting after an R³ model.

The KE-PE equivalency is nice but the problem is sorting Z and pXY vectos  out over long coasts, if I could do it I would not switch out, coasting in Keplarian system is super easy, its just so inaccurate on the jump you don't want to do it more than twice, the problem in this system is that jumping out requires using the Arctan (Arcsin or Arccos) function, which introduces about 3 or 4 digits of error. Not to bad per orbit but terrible per burn moment.

Lets get to the specific, the R³ ellipitical has a rotation from systemic E_x which with r and  a rotation in Z from plane XY from that defines Pe. Thus all Θ in the set of pXY have a unique r based on l / (1 - e cos (Θ - Θ_^pe)). pY/pX can be defined but not Z and it cannot be derived from the magnitude r because many Z and _K * (pY/pX) can define equal magnitudes so the answer is vague in R3 on the jump. Some aspect of directions of P and V have to be carried though the jump. Since the Keplarian angular moment vectors always point in the Z direction Z is locked and tan of pY/pX is always (Θ - Θ_^pe). It seems that what is missing here in R³ is a way to define Z=0 in terms of its pXY relative to Θ_^pe. Shouldn't we be able to define the parameters for the ellipse as iX²/? + jY²/? + kZ²/? = 1

The problem is that I haven't yet found a way to equate pZ to (PXY or δ_{from Pe}Θ). If only I can get an equivilancy going then problem is easily solved. If you put the ship on rails for say 45 degrees ||P|| ||V|| pops out ||Vrad| and ||Vtan|| if you know Θ relative to pE and motion then you know the sign of Vrad and if you know that then Radial velocity component is simply sign(Vrad) * unit vectors of P * ||Vrad|| But you need to position vectors for V and you need the position vectors for Vtan.

I am assuming that these are derived from the Angular moment, I assumed this was a scalar, but it seems like you are telling me its a vector quantity. How to derive unit vectors of P and V in the future (or at a future angle) based on last known physics?

I think I may have the answer provided I can convert a C centric elliptoid calculation to a f centric format.

 

    x(t) = Cx + a cos(t) Ux + b sin(t) Vx
    y(t) = Cy + a cos(t) Uy + b sin(t) Vy
    z(t) = Cz + a cos(t) Uz + b sin(t) Vz

The problem with this formula is that if C is the center we only know t at Pe and Apo; x, y, and z at (?); ||a||  but not the components<U>, ||b|| same problem, know the magnitude but not the components<V>. So this wont work.

 

Edited December 30, 2017 by PB666

[PB666]

Well here it is if anyone wants it, its a flow chart from hell but. . . .

All you are given in this example are std gravitational parameter for an object, <P> and <V>.,  ">" means they are 3-D vectors. These are in R3 space. I will say with all the math it might be better to change coordinate systems and switch back when finished.  The first thing that we know we have is a normal to <P> and <V> is <T> the tangent to both (we assume that no rocket can fly absolutely strait up or down.
  <P> x <V> = <T>    [Note that when crossing vectors that crosses in the same direction are 0 so the Px * Vx = 0 x . .  in the R3 universe we use right hand rule. P goes in as 1's, V goes in as 2' and T comes out as S

Public Sub DefineTangent(X1 As Double, Y1 As Double, Z1 As Double, X2 As Double, Y2 As Double, Z2 As Double, S1 As Double, S2 As Double, S3 As Double)
  S1 = Y1 * Z2 - Z1 * Y2: S2 = Z1 * X2 - X1 * Z2: S3 = X1 * Y2 - Y1 * X2
End Sub

(1) Once we know the normal we have the Equations for all points in the plane of the elliptical Tx * X + Ty * Y + Tz * Z = 0. In this way we have cut the size of the universe to ∞^2/3

We can define any point ellipse by its  Θ as r = l / (1 + e cos Θ). We have to realize that extraction of pXY does not suffice. Even though we know an elliptical in R³ space it is an ovoid manifold where each r is a circle about the semi-major axis only one of these planes with its velocity vectors fits the above equity 1 and only one plane has coordinates that are in the XY plane, and unless the orbit is equitorial, they are unlikely the same. As inferred above calculate l and e which immediately fall out of equations already listed above. One could use a squeeze between a recursive computational regime to get a set of coordinates whereby ||query|| = r.

The equation for any point of an ellipse from its center is given by r = (a  cos Φ)² + (b sin Φ)²  https://www.mathopenref.com/coordparamellipse.html [quoted because X and Y are not in our coordinate system but in the elliptical] however we are starting from F1 and are only given a measure of distance from the center.

x = h + a  cosΘ
y = k + b  sinΘ 

So lets review what I have F1 is on (0, 0, 0) and C is some point on the R³ that is -ea distance from 0, 0, 0 such that X² + Y² + Z², <T> is tangential and is pi raidans from the periapsis. From C the position can be mapped in the elliptical.  Im not going to worry about the position of the Pe for now, it will pop out later. I am just going to assume that i have. In an ellipse the X axis is defined as the semi-major axis. the F1 and C points  anglular begins on the ellipse are along the X-axis with the point a that crosses the elliptical closest to F1 therefore F1 and Pe lie Θ₀ direction from C and Apo.  C and Apo lie Θ_π from F1 and Pe.

Next the angle we are give to jump to is given as a Θ, lets call it Θ_j   (See above). We want to calculate the coordinates of this angle, the problem is that Θ_j will only give coordinates in the X,Y coordinate system. and to get them in x,y,z we need to use the center based program and the Angle defined by J C and F1. To do that we need the unit vectors for a [_Au below] (we already have these to get C, and if we didn't all you would have to do is normalize them), for b the is the normalize vectors of <a> x <t> crossproduct (defined above just plug the variables in _Bu). Φ is the C centered angle of J (defined above).

    Jx = Cx + a cosΦj xAu + b sinΦj xBu
    Jy = Cy + a cosΦj yAu + b sinΦj yBu
    Jz = Cz + a cosΦj zAu + b sinΦj zBu

Finally after much work there is one half the solution. We need a velocity vector, we already have the magnitude, SME = V²/2 - PE and we know SME and PE (µ/r) so that we know V. We also know V has a negative radial velocity because it is traveling from Apo to Pe (this we have just left P or Apo and jumping to a burn point. So the magnitude of the tangential velocity to <j> is 2 * Sweep/ (r * v_tan). Again this is pependicular to normalized crossproduct <J> X <T>. The remainder is the -SQRT(V²- V_tan) as V_rad and . <V_rad> = V_rad * <Ju>  <--- unit vectors of J. Sum <V_tan> + <V_rad> and now <L> is known.

So the last missing part of this are the Keplerian parameters that will be needed above.

  a = -µ / (2 * SME)
  Period = (2 * Pi * a ^ 1.5) / µ^2
  Normalize X, Y, Z, r, Xu, Yu, Zu <===== unit vector definition
  vRad = DotProduct(Xu, Yu, Zu, vX, vY, vZ) <===== This is the Pu*V
  vTan = SQRT(v^2 - vRad^2) <====tangential velocity
  omega = vTan / r <==== angular velocity 
  Sweep = (omega * r ^ 2) / 2
  Area = Sweep * Period
  b = Area / (Pi * a)
  l = b ^ 2 / a
  e = (1 - b ^ 2 / a ^ 2) ^ 0.5
  Pe = a * (1 - e): Apo = a * (1 + e) 

So with these defined the only thing that really remains is determining is the position of the center. We actually know sort of where the Pe is based we can figure this out also. This is rather a hard part, knowing how far Pe is from the Origin and From P. and from this we know where P is from the line F₁P. This creates a circle that crosses the elliptical plane once. So that we only accept the two points that cross the elliptical. The line PF₁ and that once passed its -Pu (negative of Ps unit vectors) * r where r  = l / (1 + e Cos (Θ - pi)) There is also a perpendicular to <P> at 0,0,0 that is r= 1/ (1 + e Cos (Θ - Pi/2)) and so given the two closest angles to 0 it is possible to estimate since Θ is between two angles which are 90 degrees. At least we know that Z between these is not hilariously different, such as on the North and South Pole of the geodesic, and in general Z should be in a linear range. But since we can use scalars to project vectors we can project both the reversed and perpendicular and because one knows the angle one can extend them so that the line between the two is tangential to Pe and that line is being strait the Position of Pe can be predicted.

So that the unit vectors of Pe, we can call Au since the are also the unit vectors of A (Needed above). Then the unit -e * a * Au = <C> thats everything that is required. Its done we jump. Note: about one thing, for P that are low Θ, ||P|| is less than L, P and the perpendicular to P at 0 can be used to perform the trig above instead of (K) The only catch here is when ||P|| = l, the <M> is Pe. For P very close to l the Z component of M can be used to estimate its contribution, then its fairly strait forward to guess PXY bases upon P lambda to the global coordinate system. The calculations will cause significant deviation if performed frequently, however during one cycle the max this is performed once to jump near Apo and consituitively jump to the burn initiation point. Looking at the figure above the spacecraft begins its burn at J and ends at P, then jumps to close to Apo, at very high ISP burns it _can_ correct Pe and then Jump to J. To precent the computer from slowing down the computer has a preset period of time that it test to see if it reaches P, it can easily do this by determining the normalized X,Y, and Z coordinates of the Jump initiation point (where P presumably is when it Jumps to J) When the normal vectors of the ship pass those points (generally going to be arccos Xu unless polar orbit) the ship will go into the above routined and warp to the next burn point.