Question about craft attitude / coordinate transformations

[Bill Phil]

No, that doesn't work because you don't have those digits to be able to round. That's the thing, the numbers are only exact up to so many digits / Bits, after which they just stop. If you had the additional digits to round, then you could just use those in the first place.

Maybe add one more place for the sake of rounding? Could you do that?

[Kertherina]

Maybe add one more place for the sake of rounding? Could you do that?

If i can add more digits (what i can actually do thanks to the MPFR library, although it is quite the effort compared to using a double), why should i stop at one extra? If i do this i will use something like 25 digits, if only because that is already the precision i use for position vectors and it's easiest to use one precision for all arbitrary-precision variables.

First i will try fixing the orientation with some cross products now and then though, that would be the best solution in terms of computation power needed.

[Bill Phil]

If i can add more digits (what i can actually do thanks to the MPFR library, although it is quite the effort compared to using a double), why should i stop at one extra? If i do this i will use something like 25 digits, if only because that is already the precision i use for position vectors and it's easiest to use one precision for all arbitrary-precision variables.

First i will try fixing the orientation with some cross products now and then though, that would be the best solution in terms of computation power needed.

Does the problem come from the computer using combatting values or something?

[K^2]

I do use GNU MPFR for arbitrary precision arithmetic so i can just use a higher precision variable

Be careful with that. Regardless of numerical precision, integrating equations of motion produces errors. It has to do with finite time steps. Typically, double precision computations are more than sufficient for any simulation, so you can stick to that for performance reasons. But you can only go so small on time steps before your machine can no longer handle computations. There are numerous methods for dealing with finite time steps. Velocity Verlet is one of these. It's a special case of Runge-Kutta methods. If this is the first time you hear of these, it's worth taking a look at.

Gravity problems are notorious for causing problems with numerical integration. There is no symplectic integrator for gravity. In simplest terms, it means that there is no method that perfectly prevents energy drifts in simulation. The errors can be reduced to within reason. There are a number of methods developed specifically to deal with gravity. Unfortunately, I'm not an expert in these. For my purposes, I've never had to go beyond implicit RK methods. And if you aren't looking for NASA-grade precision, you will probably be fine with a higher order explicit Runge-Kutta. RK4 will probably be more than enough.

The main reason I'm mentioning this is because everybody's first instinct when writing a simulation is to use the simple forward Euler, and that tends to produce horrible results for gravity.

[PB666]

Be careful with that. Regardless of numerical precision, integrating equations of motion produces errors. It has to do with finite time steps. Typically, double precision computations are more than sufficient for any simulation, so you can stick to that for performance reasons. But you can only go so small on time steps before your machine can no longer handle computations. There are numerous methods for dealing with finite time steps. Velocity Verlet is one of these. It's a special case of Runge-Kutta methods. If this is the first time you hear of these, it's worth taking a look at.

Gravity problems are notorious for causing problems with numerical integration. There is no symplectic integrator for gravity. In simplest terms, it means that there is no method that perfectly prevents energy drifts in simulation. The errors can be reduced to within reason. There are a number of methods developed specifically to deal with gravity. Unfortunately, I'm not an expert in these. For my purposes, I've never had to go beyond implicit RK methods. And if you aren't looking for NASA-grade precision, you will probably be fine with a higher order explicit Runge-Kutta. RK4 will probably be more than enough.

The main reason I'm mentioning this is because everybody's first instinct when writing a simulation is to use the simple forward Euler, and that tends to produce horrible results for gravity.

From first hand experience Euler equation will fail. In an N-body problem that is sufficiently complex you need to know the acceleration vectors at the beginning and end of a time step. Since the trajectory is not a simple ellipse, and therefore since accelerations are subject to complex changes with position, the integral of acceleration or velocity cannot be simplified to a simple equation.

Rephrasing what you said, the core problem with computers can be simplified to the very abstract, its not in the size of the numeric processors, the problem lies in that for a computer, time is also digital, it runs on clock cycles (square waves). By definition it will always be looking at the problem from a stepping of time. This is compounded by the fact that producing results is not intuitive to a computer, and typically processes in nature that are happening in parallel must be serially processed in a computer. it requires brute fore processing and complexity compounds the problem astronomically.

Solving the N-body problem is quite simple if one has a computer made of a core processor, assumes that no masses interact outside the system, and that masses cannot collide. Each mass has its own processor, and each body-body interaction has its own processor. At the level of 3ghz you can calculate the interactions, feed then to the body 'processor' with sums the accelerations and produces new Position and Velocity vectors, which the feeds them to the core processors that then updates.

The number of needed processors in best equipped N-body computer can be estimatated as 1 + number of bodies + Σ_i=2toN (i-1) with the time constraint being the amount of time it takes to accumulate all the vectors for each body. So therefore the minimum time with multiple bodies mean accumulation time * (number of bodies -1). As you can see as the number of bodies increase, the length of time between time intervals also has to increase.

So in trying to create an accretion problem on my PC, the amount of local space 'looks' (simple for solving the local space problem) becomes smaller if one shrinks the time. It therefore places more mass in the point-mass GM part of the equation. There have been other attempts to solve the problem, however, from what I have read they solved the problem by sweeping the computation problem by making assumptions about what happens in local space, and following the interactions of local spaces (sort of like kerbals SOI).

I don't know if this helps the OP or not, I would simply make the following simpler point

If only SOI is producing gravity then you can use elliptical/orbital equations and the solution should not vary at the origin (assign it Pe). I would not try to solve every N-body problem, It is not necessary to solve earths orbit given pluto's orbit. So what you could do with say earth is to spin it around the sun for say 20 orbits and see which 10 or 20 bodies produce the most accelerations. The same would be true fore the entire system.

From the equation I give above you can estimate how fast you want to update, because the amount of time to process a N-body problem/2 body problem is going to roughly be that number/2. E.g. a 20 body problem would be 210/2 = 105, therefore if were looking at your system every second, you would then want to look every 105 seconds and at some point you would find that there is a tradeoff (compromise the result for the most accelerative body) and therefore a need to reduce the number of bodies examined in each accumulation.

[Kertherina]

@PB666:

I can't follow you completely, but i plan to consider every major body for n-body physics for now. If the performance gets too low i can still change that, Pluto is a good example for a body that won't have a significant impact most of the time.

Here's how i currently plan to deal with n-body gravity:

I have only started to look into this (so many things to learn for a simulation like this...) so this is all subject to change, but so far Encke's method seems to be still relatively simple, yet pretty accurate. Basically you predict the unperturbed orbit which you can do exactly, and then introduce the perturbations that change the orbit. The error being proportional to the perturbations means it will also be small since perturbations are small in general so it sounds pretty neat so far.

[PB666]

@PB666:

I can't follow you completely, but i plan to consider every major body for n-body physics for now. If the performance gets too low i can still change that, Pluto is a good example for a body that won't have a significant impact most of the time.

Here's how i currently plan to deal with n-body gravity:

I have only started to look into this (so many things to learn for a simulation like this...) so this is all subject to change, but so far Encke's method seems to be still relatively simple, yet pretty accurate. Basically you predict the unperturbed orbit which you can do exactly, and then introduce the perturbations that change the orbit. The error being proportional to the perturbations means it will also be small since perturbations are small in general so it sounds pretty neat so far.

This is long so if you are TL:DR type go to the next post or thread:D

I guess if I was trying to over simplify the problem, lets say we are going to study Jupiters orbit, and then we placed all the objects in the solar system inside jupiters orbit, the closer the objects were to the sun, the easier it is to solve the orbital problem (except for the tidal issues, which might actually get worse, lol). When they are in the sun its a two body problem. By the same token, if we move all the masses in the solar system to the periphery of Sun's sphere of influence, the problem is also now very easy. Jupiter goes flying forward until it approaches some mass. Assuming that we start with a 'good' position and velocity, what we really want to know is what the acceleration (or some cozy function thereof) is going to be over the next time interval of simulation.

I think I left something out of the computation on the other page, the calculations for the interaction vectors for each interaction are not as simple as an accumulation (explained at the end of the post). Lets give every object a P (im going to use matrix but not matrix math)

Lets say object zero at any moment can be defined as ⁰O_t = {⁰M, ⁰P_t, ⁰V_t, ⁰A_t, ⁰E_t} Mass, Position, Velocity, Acceleration, Elliptical (moment). This sounds easy but it really is not. Lets say ⁰O is the central body of the system a massive star. The elliptical will remain an elliptical until we apply A*t, then both the PV&E will change.

is the central body

⁰P_t=0 = {⁰X_0,⁰Y₀⁰Z₀}

⁰V_t=0 = {^0XV_0,^0YV₀^0ZV₀}

⁰E_t=0 = {⁰Pe_0,⁰I₀, ⁰a₀, and the various Theta (arguments)}

These three mean that first you have to define a reference system, which most logically means it has to be the systems center, which is going to be different from ⁰P₀. This also means that ⁰V <> 0.

This means you need to calculate the systems mass,^easy and the center of gravity. Since Y will be perpendicular to the plane one needs to decide on a plane (either a massive satellite or the weighted plane between all satellites). So in the Solar system who do we have to include to be accurate.

1. Jupiter-yep (and its satellites)-yep [hmm maybe we need a separate processor to integrate these]

Let suppose they were all in a line this would be the point in the suns 'orbit' that its center would be a maximum distance from the systems center its apoapsis.

2. Saturn-Yep (and its satellites)-yep

3. Neptune-Yep

4. Uranus-Yep

5. Venus-Yep

6. Earth-Yep

7. Mars-Yep

So the question then becomes, how many more masses need to be added before the system center stops moving. This not the hypothetical, its how much more mass needs to be added before your floating point digit is no longer relevant. For example the asteroid belt objects mass are sort of balanced, in addition also objects that are roughly pluto's that are so numerous that, in essense pluto's imbalancing effect is negligable. The other non-functional limitation is that because of the limitations in our understanding of mass and gravity in our solar system, the suns mass is only accurate to the 4th decimal place. For the sake of argument lets extend it to the fifth. It may be true that this is because we don't know G past that many places, but as we move further from the sun or earth, these calculations become less certain. Since we are placing all our planets in a line we can use the semi-major axis (a) as distance to ⁱX, 0 and the center will be ⁰O. Jupiter moves suns orbit some 743,000 km (this is an example so don't fault me if these calculations are not perfect, the problem I worked on was a star undergoing accretion in its planetary disk) To do this right you need to accumulate the total weight of each system (I didn't). I don't think it matters were you place the satellites around the planets.

So if we do this we see that each of the 4 largest planets listed move the system center away from the sun on the magnitude 8 scale. Adding Earth drops the outward motion down to magnitude 5. Mars magnitude 4. Mercury 3.

We can see here that the distance between relative distance added by mars between the system center and neptune is now below relative error of mass determination for Neptune or the Sun. We can sweep that under the rug for a moment because we can argue that mu is precise even if G and ⁱM is not. Even there, mu for bodies is based on alot of assumptions and so if we carry this more than 2 to 3 more magnitudes we really are having a Don Quixote moment. If we look at the attractive affect of Pluto on the sun, its 1/10⁷ that of jupiters, it about 1/1000th the error in Jupiters mass determination. Yeah sure mu can extend this, but only two digits or so.

[TABLE=width: 360]

[TR]

[TD][/TD]

[TD=colspan: 2][/TD]

[TD][/TD]

[/TR]

[TR]

[TD][/TD]

[TD]Center-O distance[/TD]

[TD]A-effect

[/TD]

[TD][/TD]

[/TR]

[TR]

[TD][/TD]

[TD]Max Dev (m)[/TD]

[TD]on sun[/TD]

[TD]dDev/O (m)[/TD]

[/TR]

[TR]

[TD]Jupiter[/TD]

[TD=align: right]742621394

[/TD]

[TD=align: right]2.09E-07[/TD]

[TD=align: center] 742621394

[/TD]

[/TR]

[TR]

[TD]Saturn[/TD]

[TD=align: right]1151676894[/TD]

[TD=align: right]1.85E-08[/TD]

[TD=align: right]409055500[/TD]

[/TR]

[TR]

[TD]Neptune[/TD]

[TD=align: right]1383038072[/TD]

[TD=align: right]3.38E-10[/TD]

[TD=align: right]231361179[/TD]

[/TR]

[TR]

[TD]Uranus[/TD]

[TD=align: right]1508129542[/TD]

[TD=align: right]7.03E-10[/TD]

[TD=align: right]125091470[/TD]

[/TR]

[TR]

[TD]Earth[/TD]

[TD=align: right]1508573705[/TD]

[TD=align: right]1.78E-08[/TD]

[TD=align: right]444163[/TD]

[/TR]

[TR]

[TD]Venus[/TD]

[TD=align: right]1508834535[/TD]

[TD=align: right]2.77E-08[/TD]

[TD=align: right]260830[/TD]

[/TR]

[TR]

[TD]Mars[/TD]

[TD=align: right]1508907523[/TD]

[TD=align: right]8.24E-10[/TD]

[TD=align: right]72988[/TD]

[/TR]

[TR]

[TD]Mercury[/TD]

[TD=align: right]1508916876[/TD]

[TD=align: right]6.57E-09[/TD]

[TD=align: right]9353[/TD]

[/TR]

[TR]

[TD]Pluto[/TD]

[TD=align: right]1508955363[/TD]

[TD=align: right]2.52E-14[/TD]

[TD=align: right]38487[/TD]

[/TR]

[TR]

[TD]Ceres[/TD]

[TD=align: right]1508955559[/TD]

[TD=align: right]3.67E-13[/TD]

[TD=align: right]195[/TD]

[/TR]

[/TABLE]

We can see from this that Ceres is not going to significantly alter the position of the sun in any way that we can measure. Its mass effect on the suns acceleration is at the border line of our knowledge of the sun's or Jupiters Mu the primary determinants where the Sun is being pulled and how far its center is from the system center. So the issue is whether one should waste processor time on anything smaller than pluto when calculating the orbits of distal objects.

So now we have a system center. We can then get the position information for the Sun, but now we need to assign velocity vectors, yikes, which we do not have. We can assign acceleration but until we run our simulation for a time we wont have an dynamic equilibrium representation of velocity. If we take our sun at apogee of its orbit we could use the elliptical equation to derive a velocity but the effects of inner planets is likely only to muddle the argument, you would need to focus on the big four. The velocity vector however is idealized, there are inclination vectors that also need to be considered (Imagine me handwavign vigorously- substitute some intense calculations). We therefore need to know what the Suns orbital motion was like the last time it reached Apogee as a good starting point. It is assumed that the motion of the Sun defines the XZ plane thus Y is perpindicular to the line that bisects the velocity vector at P creating a plane that Y is perpindicular X and Z can be assigned relative to the periapsis of Jupiter. From this map we can select N=10 or so objects to detemine the next value A objects

That being done you can move on.

⁰A₀ = G * Sum_{i = 1 to N} M_i * (ⁱP₀ - ⁰P₀)² That is the initial acceleration vector, if time is large enough and A is also large enough, you will need to go back and calculate ⁰d₁ and detemine the change in acceleration and calculate the average acceleration vector at over dt. This is why is best to keep the number of effectors low. Because you are going to have to do the same calculation on all objects in the system, but each object will have different effectors depending on is position in space (e.g. Pluto does no have to worry about say ceres when its on the other side of the Sun, but it does have to worry about the Sun and Jupiter everywhere).

I did underestimate the processing power needed, because as you can seen the critical calculation A is not a simple accumulation, as acceleration increases, there will be a need to increase the number of passes (improvements on A) . The ideal computer needs interaction subprocessors (multiple) to assess each dimension and parse distance into X,Y,Z, it then needs to square the values and send it to the interaction processor which then accumulates the values and takes the square root and (assuming that the other ⁱP₁ (i>0) do not change significantly from tempⁱP₁), acceleration would be recalculated also at ⁰Pt=1u (dP = V₀u + 1/2au²)

and then some formula approximating acceleration between 0 and t = u would be applied.

Then there is one final caveot, the t * a metric will be a good allocator of clock-cycle resources except when two objects approach each other (ruling out gravitational waves some collision might send to each other). As this occurs you may have to go into the process with a much more refined unit of time. Because error>meter changes of orbit (e.g. Pe) could have profound changes if the two objects approach, because neither object can be treated like point sources anymore. In the event that two object distances are below the maximum added altitude of terrain one runs into uncertainty outcomes about the result, more calculable interactions would be a atmospheric interaction. Unless you are also following the rotation of the planet and its wobble in its system you have Shrodinger's cat. I would also consider you lucky, my simulation never got this far, I estimated I would need a million processors running a year to get a nice collision. If you have alot (billions) of objects you simply cannot calculate all the close calls, you have to use probabilistics.

I

[K^2]

@PB666:

I can't follow you completely, but i plan to consider every major body for n-body physics for now. If the performance gets too low i can still change that, Pluto is a good example for a body that won't have a significant impact most of the time.

Here's how i currently plan to deal with n-body gravity:

I have only started to look into this (so many things to learn for a simulation like this...) so this is all subject to change, but so far Encke's method seems to be still relatively simple, yet pretty accurate. Basically you predict the unperturbed orbit which you can do exactly, and then introduce the perturbations that change the orbit. The error being proportional to the perturbations means it will also be small since perturbations are small in general so it sounds pretty neat so far.

Good plan. But even then, doing proper numerical integration on plaents' paths is a pain. I would recommend starting with planets on rails, and seeing how well you can simulate trajectory of the spacecraft. You can grab data for planets in the Solar system from Horizons. Make sure to set coordinate system to Sun's barycenter. Otherwise, you'll have some fictitious forces due to accelerated frame of reference. And you can compare your simulation to trajectories of some asteroids and comets. If you can get decent simulation for an asteroid, you'll know that your integration method is good. Then I'd give a crack to integrating motion of planets using Encke's.

[Kertherina]

@PB666:

I don't have much to say about all this since i have yet to dig into this topic, but thank you so much for writing that out, that was really helpful.

@K^2:

Currently i have set the origin of the coordinate system to the starting point of the Sun. I don't see how this leads to problems though, i have a fixed frame of reference and calculate everything from there. That shouldn't be a problem, or is it?

Edited March 16, 2015 by Kertherina

[K^2]

Currently i have set the origin of the coordinate system to the starting point of the Sun. I don't see how this leads to problems though, i have a fixed frame of reference and calculate everything from there. That shouldn't be a problem, or is it?

So long as the Sun's starting velocity isn't zero, it should be fine. But choosing center of mass of the Solar System as origin makes a lot of the computations easier. It makes it easy to place initial position and velocity of the Sun to counterbalance all of the other objects in the System.

What you want to make absolutely certain of is that whatever coordinate system you've chosen, center of mass either remains at rest or moves at uniform velocity. This should take care of itself if you simulate the system and all of the planets. But if you plan to use pre-determined positions of planets for testing, you will have to make sure that the coordinate system choice is appropriate and consistent.

[Kertherina]

Just wanted to let you know, it works like a charm. I can easily afford to reorthonormalize with Gram-Schmidt every couple rotations (perks of having a very computation heavy simulation ^^), doing it every 10 rotations atm, where the floating point errors are obviously pretty much nonexistent.

[K^2]

Gram-Schmidt sounds like overkill for a 3-basis. Cross products can be computed with a handful of SSE2 instructions. What are you using for rotations? If you aren't scared of low level code, you might be better off writing rotation code that does orthonormalization during rotation in a single pass, saving you trouble of doing any other checks. Branching, in particular, can be avoided all together, and that's probably costing you more than the actual operations.

[Kertherina]

What are you using for rotations? If you aren't scared of low level code, you might be better off writing rotation code that does orthonormalization during rotation in a single pass, saving you trouble of doing any other checks. Branching, in particular, can be avoided all together, and that's probably costing you more than the actual operations.

I'm using quaternions. And since i'm still pretty new to programming (this is my first project bigger than "calculate some prime numbers"), low level coding does indeed sound scary.

Besides, if performance were a problem i could just reorthonormalize only every 100 or 1000 rotations. The error shouldn't be noticeable then either, the only reason i am doing it every 10 rotations is because there is no noticeable difference in total simulation performance.

And branching...well, again, as long as i don't notice a drop in performance it should be fine.

[K^2]

Oh, if you are doing all of your transforms with quats, you should store the transformation of your space ship as a quaternion as well. In that case, all you have to do is renormalize the quaternion once per tick. To get the body axes of the ship, just take basis axes and transform them with ship's quaternion. The error from a single transformation is minimal, and you don't need to worry about enforcing orthogonality.

By the way, so long as your quaternions are properly normalized, there is a good shortcut for computing the qvq* that you can use. By using the commutator operator defined as [a,b] = ab - ba.

qvq* = qq*v + q[v,q*] = v - q[v,q]

The scalar part of the [v,q] is zero and the vector part is twice the cross product v x q.

[Kertherina]

Thank you, that sounds all very helpful. Now i just need to think about how exactly to store the ship's orientation in a quaternion.

[K^2]

I want to apologize in advance if I'm reiterating a bunch of things you know already. It's a little hard for me to gauge your precise background. But hopefully, at least some of this is helpful.

Think about how you are storing the orientation of the ship now. You have the XYZ vectors of the ship written in world's coordinates. If you want to take some vector (x, y, z) in ship's coordinates to world coordinates, you compute a new vector xX+yY+zZ. But that's exactly the same thing as writing a matrix with columns (XYZ) and multiplying it by a column vector (x, y, z)^T. Storing the body axes is equivalent to storing a rotation matrix. Instead of storing a rotation matrix, you can just as easily store rotation quaternion.

Lets break it up a bit more. Suppose you have quaternion s that tells you how to transform some vector v in ship's coordinates to world coordinates. Write u = svs*. Now, suppose, you want to apply a rotation quaternion q to that vector u. That q might be precisely the change in rotation of the ship that you are currently applying to the axes. The fully transformed vector is now quq* = qsvs*q*. And because quaternions follow standard conjugation rules, s*q* = (qs)*. In other words, you can first compute a new quatrnion r = qs, and then apply that to original vector v as rvr*.

Starting from the beginning, and looking at transformations as they come in at each time tick q₁, q₂, q₃..., you start with, say, vector X = (1, 0, 0) as one of your body vectors. You then update it to q₁Xq₁*. Then you update that vector on the next step, and you get q₂q₁Xq₁*q₂*. And so on. Well, instead of keeping track of that X on each tick, all you really need to do is keep track of that product on the left, and call that your ship's transform quaternion. s = q_n...q₂q₁. At each time step, you just need to multiply by one more quaternion on the left. The numerical errors will manifest in normalization of that quaternion drifting, and it's a good idea to just renormalize it every once in a while.

I'm more used to single precision quaternion math, and that tends to drift annoyingly fast. So with single precision, I'd recommend renormalizing on each time step. Especially, since that's really cheap to do. With double precision, it might not matter as much, and you can get away with renormalizing much less frequently, but it's really inexpensive on modern CPUs, so you might as well just do renormalization on each tick.

[Kertherina]

I want to apologize in advance if I'm reiterating a bunch of things you know already. It's a little hard for me to gauge your precise background. But hopefully, at least some of this is helpful.

I am new to all of this and learning as i go, so thanks again, you are really, really helpful and this cleared things up yet again

What i wanted to say was i need to think about what exactly the real part of the quaternion would stand for since that would be the roll information i currently store with the separate starboard vector. So basically, how is the craft rolled when the real part is zero since that would be the basis for everything and i need to know how to initialize it at the beginning of the simulation. I just don't have a (visual) grasp of it just yet.

As for errors, there is so much stuff going on that i can easily just normalize each tick, it really makes no noticeable difference.