Using linear Transformations for modeling spaceflight, how-to

[PB666]

First off, let me define what the need is.

When you launch from a planet you are basically using a spherical cartesian system. Longitude and latitude are replaced by a non-newtonian coordinate system. Once you get into space your orbit has no fixation relative to the surface unless your are in GTO. These orbits are typically oriented to the position of the sun at the northern vernal equinox. This system can have two orientations, the one orientation is relative to earths poles, the second can be relative to the normal of a surface pointing in the direction of the sun (essentially consider a sphere of radius r, its the path along the surface that intercepts a line between the Earths' center and the sun). Within the EM system there can also be other coordinate systems. There could be a lunar-centric/earth-referenced system, lunar-centric/solar-referenced, L1centric, L2-centric. Then there are the myriad of planetary and satellite coordinate systems. Going between points in the solar system can be cumbersome if you use earth rotational coordinate system marked on the vernal equinox because everything (including Earth's orbital inclination, by ~ 23'). There is a preferential plane of motion (marked largely by the motion of Jupiter) that is used to define the plane of the system). And some day a spacecraft might voyager might use a galactic coordinate system based on an orientation to some far off galactic center close to the plane of the galaxy.

Each coordinate system has a basis. The earth's X,Y, Z has a line that goes through Greenwhich and some point in the Atlantic near West Africa marks the X=r point Z = r close to the geographic north pole.

If you read the wiki page on Change of Basis you are given the following: c = cosine, sin = sine. α , ß, gamma are the three angles by which any vector can be defined by a change of x -> x' (rotation around z), z->z' (rotation around x'), and x'->x'' (rotation around z') all counterclockwise.

Source wikipedia - change of basis.

While this is correct, this basis, quoting wikipedia.  "R be a new basis given by its Euler angles. The matrix of the basis will have as columns the components of each vector. Therefore, this matrix will be [above matrix]". But in fact this is not very useful, because by the fact that you have the angles that define the rotations, you pretty much already know where R is with respect E³. This matrix tells you what you already know. By saying its not very useful, its useful on the 'return trip back to the previous (home) coordinate system, but not so useful on the leg of a journey that creates the need for the return home conversion. 

But what you really, really badly want to know is how to map a vector already mapped in E³ in your new basis  R and here is what the wiki article says. " Again, any vector of the space can be changed to this new basis by left-multiplying its components by the inverse of this matrix."

This single quotations underlies a fairly complex problem that can be simplified with certain information that wikipedia article does not readily provide. 
The first bit of information one needs to know the determinant of the Matrix above. To know that you have to how that matrix was derived.

The matrix was derived from three matrices which are principally defined here;https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations  and converted into symbolism, defined above, below.

https://i.imgur.com/bwHo4O9.png

One law of determinants is that if all the determinants are 1 then the multiplication of the Matrices is also 1. https://en.wikipedia.org/wiki/Determinant

With this rather innocuous fact one can derived the inverse matrix.

https://i.imgur.com/8ftGcsI.png

This is of course not what the math actually looks like.

The only way to know for certain if the proposed R-1 is to multiply R^-1R =  I (the identity matrix, I₃ below)

 

This is the first stage after multiplying everything out and getting rid of extra parenthesis.

I_{row , column} or I_i,j i=row# and j= column#

I_1,1 = c²1c²3 + s²1c²2s²3 - 2cs1c2cs3 + c²1s²3+s'1c²2c²3 + 2cs1c2cs3 + s²1s²2
I_1,2 = c²1c2cs3-cs1c²2s²3 - s²1c2cs3+cs1c²3 + cs1s²3 + s²1c2cs3 - c²1c2cs3 - cs1c²'2c²3 - cs1s²2    
I_1,3 = c1s2cs3 - s1cs2s²3 - c1s2cs3 - s1cs2c²3 + s1cs2
I_2,1 = cs1c²3 + c²1c2cs3 - cs1c²2s²3 - s²1c2cs3 + s²1c2cs3+cs1s²3 - c²1c2cs3-cs1c²2c²3 - s²2cs1    
I_2,2 = s²1c²3 + c²1c²2s²3 + 2cs1c2cs3 + c²1c²2c²3 + s²1s²3 - 2cs1c2cs3 + c²1s²2    
I_2,3 = c1cs2s²3+s1s2cs3-s1s2cs3+c1cs2c²3-c1cs2
I_3,1 = c1s2cs3 - s1cs2s²3 - c1s2cs3 - s1cs2c²3 + s1cs2
I_3,2 = c1cs2s²3 + s1s2cs3 - s1s2cs3 + c1cs2c²3 - c1cs2    
I_3,3 = s²2s²3 + s²2c'3 + c²2
https://en.wikipedia.org/wiki/File:Matrix.svg

Source: Wikipedia - Identity Matrix
An Example (the easiest).

I_3,3 = s²2s²3 + s²2c²3 + c²2 (converting the identity subscripts to angles, since excel does not handle text markup well)
I_3,3 = s²_ßs²_Γ + s²_ßc²_Γ + c²_ß (rearranging)
I_3,3 = s²_ßc²_Γ + s²_ßs²_Γ + c²_ß  (factoring)
I_3,3 = s²_ß(c²_Γ + s²_Γ)+ c²_ß  (Replacing, c²Θ+s²Θ = 1)
I_3,3 = s²_ß(1)+ c²_ß  (Reducing)
I_3,3 = s²_ß+ c²_ß  (Replacing, c²Θ+s²Θ = 1)
I_3,3 = 1

Pretty much terms disappear in the diagonal finally as replacing a final c²Θ+s²Θ as 1, and off diagonal by a final subtraction of x + -x = 0 resulting in the identity matrix.

An Example (the tie for hardest, witchcraft no less)

I_1,2 = c²1c2cs3-cs1c²2s²3 - s²1c2cs3+cs1c²3 + cs1s²3 + s²1c2cs3 - c²1c2cs3 - cs1c²'2c²3 - cs1s²2 (getting rid of ambiguous symbolism and markups, c = cos, s = sin)
I_1,2 = c²_αc_ßc_Γs_Γ - c_αs_αc²_ßs²_Γ - s²_αc_ßc_Γs_Γ + c_αs_αc²_Γ + c_αs_αs²_Γ + s²_αc_ßc_Γs_Γ - c²_αc_ßc_Γs_Γ - c_αs_αc²_ßc²_Γ - c_αs_αs²_ß  (rearranging)
I_1,2 = c²_αc_ßc_Γs_Γ + s²_αc_ßc_Γs_Γ + c_αs_αc²_Γ + c_αs_αs²_Γ - c²_αc_ßc_Γs_Γ - s²_αc_ßc_Γs_Γ  - c_αs_αc²_ßc²_Γ- c_αs_αc²_ßs²_Γ - c_αs_αs²_ß  (factoring)
I_1,2 = (c²_α + s²_α)c_ßc_Γs_Γ + c_αs_α(c²_Γ + s²_Γ) - (c²_α+ s²_α)c_ßc_Γs_Γ - c_αs_αc²_ß(c²_Γ + s²_Γ) - c_αs_αs²_ß (Replacing, c²Θ+s²Θ = 1)
I_1,2 = (1)c_ßc_Γs_Γ + c_αs_α(1) - (1)c_ßc_Γs_Γ - c_αs_αc²_ß(1) - c_αs_αs²_ß (Reducing)
I_1,2 = c_ßc_Γs_Γ + c_αs_α - c_ßc_Γs_Γ - c_αs_αc²_ß - c_αs_αs²_ß (rearranging)
I_1,2 = (c_ßc_Γs_Γ - c_ßc_Γs_Γ)+ c_αs_α  - c_αs_α(c²_ß + s²_ß ) (rearranging & factoring)
I_1,2 = (c_ßc_Γs_Γ - c_ßc_Γs_Γ)+ c_αs_α - c_αs_α(c²_ß + s²_ß ) (replacing, c²Θ+s²Θ = 1; x + -x = 0 )
I_1,2 = (0)+ c_αs_α - c_αs_α(1)  (reducing)
I_1,2 =  c_αs_α - c_αs_α  (replacing; x + -x = 0 )
I_1,2 =  0

[If you were having trouble sleeping the above 12 lines of linear algebra is much more effective than any sleep medication, I discovered this while watching lectures on quantum entanglement; if that still doesn't work try solving the unsolved elements for 1 or 0]
 

So how does one use the inverted Basis Matrix (R^-1), wait, why would someone, anyone want to use this. Chances are you already are. When you reach the edge of an SOI in KSP, you are going to move to a Kerbols SOI.

Lets give an example, you are heading radially outward along a prograde AoA from Kerbin, north is above your head, and Kerbol is to your right. You reach the SOI, now instead of heading radially at V_kerbin  you are now heading you are traveling a different velocity, circum Kerbol at a Velocity in the opposite direction (Kerbins orbital velocity - the last radial velocity inside kerbin upon exit). You also might notice you are accelerating toward kerbol, this is because your new velocity creates a centripedal acceleration that is below the current force of gravity from Kerbol. (this would be a useful problem except KSP likes to map the Kerbols coordinate system versus the Kerbins, so . . .we need a real system)

There's two sets of problems here though, acceleration is a moment problem that we do not have to worry about at the change of basis, the velocity 'blip' can be converted using the inverse Matrix and the position vector can be determined by determining the current position of planet  from its star, the distance of the ship from the planet, the path the planet travels. and the current angle of the vessel from the path at the paths exit. Basic trigonometry.  Since velocity in a 3-D coordinate system can be a vector independent of the coordinate system, we can superimpose the stars coordinate system on planets center, ascertain the change of basis angles, the apply the inverse matrix to the velocity vector to get its new velocity vector in the stars coordinate system (the magnitude of the velocity is still referenced to the planet) so the planets velocity vector relative to the star needs to be added to the new velocity vector to deduce the ships velocity relative to the star while still inside the planets SOI.  Since it is easier on an escape vector to calculate both the planets position (and therefore its gravity vectors) the stars position (and therefore its gravity vectors) one can simply change coordinates once an escape trajectory is realized. One way to do that is to change basis as the ship crosses the line from the planet to the star of from the star through the planet and out the other side (outer system transfers). Anyway

 ^ev = Xe₁, Ye₂, Ze₃ and ^rv = X'r₁, Y'r₂, Z'r₃

You start with your planets X, Y , Z coordinates in the local system they would looke like this https://www.khanacademy.org/math/precalculus/precalc-matrices/matrices-as-transformations/e/multiplying_a_matrix_by_a_vectorThis is a process otherwise known as a transformation, but in the new coordinate system we have a new set of axis (so in the video you would draw another set of axis and plot your transform in that new set of axis)
X' = Xc_αc_Γ - Xs_αc_ßs_Γ - Ys_αc_Γ + Yc_αc_ßs_Γ + Zs_ßs_Γ
Y' = -Xc_αs_Γ - Xs_αc_ßc_Γ - Ys_αs_Γ + Yc_αc_ßc_Γ + Zs_ßc_Γ
Z' =  Xs_αs_ß - Yc_αs_ß + Zc_ß

And  written finally out in a fashion that even a college graduate might understand.

Everyone wants a summary.
X' = X * cos α *cos _Γ   -   X * sin α * cos ß * sin _Γ   -   Y * sinα * cos_Γ  +   Y * cos α * cos ß * sin_Γ   +   Z * sinß * sin_Γ
Y' = -X * cos α * sin_Γ   -    X * sinα * cos ß * cos_Γ    -   Y * sinα * sin_Γ   +   Y * cos α * cos ß * cos_Γ   +   Z * sinß * cos_Γ
Z' =  X * sin α * sin ß  -    Y * cos α * sin ß               +  Z * cos ß

Suppose that you only rotated around one axis, say the Z axis. The above matrix is still good: not that element R^-1₃₃ = cos ß if no rotation along X' axis. If then ß = 0 and the cosine of zero is 1 while the sin 0 is zero. So whatever Z value you have, Z' = X * 0 + Y * 0 + Z * 1, in actuality, if you never make a rotation ß around the X' then basically you apply the A = matrix above in which the α = (your α + _Γ).
Does this happen? if your thrust changes the position of your periapsis but not the inclination relative to your external reference frame, then before the thrust your elliptical reference frame is at angle α relative to the basis (E³) and the movement of the periapsis move _Γ  relative to the old elliptical, to place the crafts position and velocity vectors in the new elliptical reference frame you could use the equation above applying α and _Γ.

Summary, this is as hairy as the math gets, but using it you can follow position and velocity in two coordinate systems (or even three as the above shows) at one time, moving flawlessly between different systems of plotting. I think its pretty obvious how to got to a R coordinate system back to the E basis system (the equations are very close).

 

 

 

 

 

 

 

 

 

 

 

[Nuke]

in before quaternions.

[PB666]

lets not go there, matrix multiplication and minor determinants are bad enough.

[PakledHostage]

2 hours ago, PB666 said:

lets not go there, matrix multiplication and minor determinants are bad enough.

Quaternions are easier in a lot of ways, IMHO. The math isn't that complicated if you want to "roll your own", but there are plenty of libraries out there too.

[YNM]

Why not spherical coordinates ? Surface rotation is just a sphere with some angular velocity. Newtonian derrivation of Kepler laws are done in spherical coordinates as well.

Quaternions is useful if you want to have rotations too. In fact, quaternions is the only way you can have rotations in 3D or so. I think this manifests somewhat in orbital parameters - you need radius, inclination, arg. of periapse and then anomalies.

Edited March 6, 2018 by YNM

[PB666]

5 hours ago, YNM said:

Why not spherical coordinates ? Surface rotation is just a sphere with some angular velocity. Newtonian derrivation of Kepler laws are done in spherical coordinates as well.

Quaternions is useful if you want to have rotations too. In fact, quaternions is the only way you can have rotations in 3D or so. I think this manifests somewhat in orbital parameters - you need radius, inclination, arg. of periapse and then anomalies.

When ß is non-zero, α is essentially (or can be transformed easily to) the longditude of the ascending node, ß is the inclination, Γ is the argument of the perisapsis all amp to the reference frame. If you know Γ and you know your crafts position in R³ then you also know your position in E³ . The difference is that in E³ you map to a spherical coordinate system, and in R³ you map to a polar coordinate system as long as you know the E³ coordiantes of your perisapsis. Once you have ß fixed, then α globally references all possible ellipses with inclination ß and so set of all Γ can map there semi-major axes comprehensively. Imagine it like this ß is a rotation about an X_αe₁, Y_αe₂. Γ₀ is Xe₁, Ye₂ in the sense that X'r₁ = 1 and Y'r₂  = 0, for all coordinate in the plane defined by Z'r₃ = 0. If you know where Γ₀ is you can know where any point in the plane is, including the argument of periapsis (say Γ₁) and the true anomoly (say Γ₂) and all past periapsis in the same plane (say Γ₃ - Γ_n)  as long as you can map back to Γ₀. Since these are all unit vector then radius is simply a position vector in the 2-D plane orthogonal to r₃.  Orbital parameters are all preserved.

So why use this, the reason is accuracy, on the computer if you are constantly using the keplarian laws of motion, then every process is a function of sin, cosine followed by inverse sin and inverse cosine, this is highly inaccurate if use for each moment of acceleration. It is better to use X, Y, Z (newtonian) coordinate system to prevent drift of the orbit, i've  tested this and the drift is something like 10-100 times faster. The problem is when you are working with piddly ION drives where the moment of thrust is tiny and the burns can last years the relative variance is problematic. IOW, once you are in a distinct keplarian orbit your are fine, and if you are perform short  bursts of thrust you probably don't have a problem. But with ION driven spacecraft, you do not want to be altering your periapsis continuously using a keplarian equations, it will get filthy. So that once we have placed the craft in a Z = 0 plane, the math greatly simplifies and as long as we stay in that plane we can continue mapping the coordinates of our periapsis Γ relative to X_αe₁, Y_αe₂.

Quaternions are an artifact that turn out to be useful in computer graphics, the problem is the goal here is not to create computer graphics or simulate physics on the computer in 3D. The specific goal is to use a set of coordinates that are most easily tested using the La grangian (the least action required) in this case to achieve desired orbit with regard to energy spent and time. Because of this we want drift to be minimal per small unit of time. 

Let me give an example. A goal of burning would be to project a vector in space in a given direction. In the case we can set an exit position in E³ (say a Mars orbit intercepting Hohmann transfer ellipse that projects off of Earth's SOI (or hill sphere). There are many small referencing problems here. But the first problem is that this exit position in E³ maps to a coordinate position defined at the solar system. We can define the points along the Hohmann transfer in the heliocentric coordinate system let's call it H³. Xh₁, Y'h₂, Zh₃ And like the motion of the Earth. We can also project onto the sphere of Earths SOI where the Hohmann ellipse will intersect that sphere. So these are relatively easy maniupulations and can be done anyway you like, but then you need to convert to Earth's coordinate system, again this can be done by any number of techniques, it's not important. The reason these operation are 'trivial' is that you only do them once, if you were to repeat them 100 times a second for 10 years, the method you choose would not be trivial. But then there is a single ellipse within that system (E³ is infinite but SOI is not) that projects out of the SOI.

But the ellipse that leaves, given that impulse is not ever instantaneous in the real world, a set of evolving ellipses with time, to be accurate many different ellipses in Z'=0 all with very small intervals of time. So given any small range of acceleration (i.e. a =10^-4 to 10^-2) a large amount of dV. The ISP, dV and amount of time burning can change. The question is when is best to burn, along what angles, and at what ISP (N = KW * eff./exhaust velocity) to achieve the fastest exit times with the least amount of fuel burnt. So for any given point along an ellipse we need to momentarily optimize parameters to achieve the highest benefit for the lowest cost. These transformations put me where I want to be, this is what I am going to use, I just thought I would present them here. Take it or leave it.