Maybe we should just standardize on vectors of rotation for everything, and kiss quaternions good-bye. The point of using quaternions was that they were supposedly easy to interpolate; but by now we know better.

Well, here is the thing. I don't exactly understand why people use unit quaternions when the same information can be contained in a vector. I think it has to do with this though: A vector has all the information you need for the rotation, but there is no vector operation that will apply the rotations correctly. You have to use a rotation matrix. Consider the un-rotated default position of an object.

In quaternion notation, that is (1,0,0,0). In vector notation, it is (0,0,0). Now say you want to rotate Pi radians about the x-axis. In quaternion notation, this rotation is (0,1,0,0). In vector notation, this would be (Pi, 0, 0). This seems ok so far. We have the axis and the angle through which we rotated. Now however, consider an additional rotation about the y-axis of Pi radians. The combination of these two rotations is equivalent to a single rotation of Pi radians about the z-axis. Just visualize rolling a plane(airplane) upside down, and then pitching it over 180 degrees. It is now facing the opposite direction from where it started, but it is facing upright again. Just like if it had yawed 180 degrees.

You can see this by doing quaternion multiplication of (0,1,0,0)*(0,0,1,0)=(0,0,0,-1). However, if we do the same thing we did last time with the vector, we get (Pi,Pi,0). This is a rotation, but it does not take the object where it is supposed to be. Cumulative rotations don't work with the vector operations. This is why you need the rotation matrices.

You can "store" orientation as a 3 element vector, but to actually perform the rotation, you need to either construct quaternions of both the current orientation and the rotation, or construct a rotation matrix that can act on the 3-element vector.

According to wiki-pedia, quaternions are preferred because they don't suffer from something called gimbal lock that can happen with rotations. However, I have not been able to find any place on the internet that explains wtf gimbal lock is.

Likewise, you can compute the delta_theta_vector as a vector, but to actually apply it as a rotation, you need to construct a quaternion from it. This is why angular_velocity and angular_acceleration can be (and should be) vectors instead of quaternions. They are instantaneous derivatives. They aren't an actual finite rotation.

Sadly, quaternions are necessary for rotations to work correctly, but I currently can only do angular displacement calculations using vectors. So, unless I can figure out how to do it with quaternions(not sure its even possible), it is inevitable that we need to switch between the two... calculating the 3-element vector representing the rotation and then turning that into quaternion that we apply to the current orientation.

First we need to answer how many gravitational bodies we intend to simulate,

This is a HUGE question. The I we just do objects large enough to have assumed a mostly spherical shape, we already have hundreds of objects. If we start doing anything large enough to impact your space ship in a collision, then we basically need everything larger than a medicine ball. That is a number I wouldn't even like to guess at, but if I had to, I would say at least millions, perhaps billions of objects. You have the rings of the gas giants, the asteroid belt, the ort cloud (though it may be too far out for us to care about), comets.

I am actually starting to wonder if it might make more sense to treat each planetary system as its own "system". We would obviously need some complicated sector grid. With differentially rotating sectors within the system such that planets always stay in the same grid as they go around the sun. Aside from the problem of gridding, this would allow resources to be focused entirely on the system the player is in allowing us to model it in detail. This would mean a "loading" screen in-system, but if we are smart, the loading screens will be far away from anything of interest, and players can use jump-gates for in-system travel between planets and, in most cases, never see the loading screen as it is covered up by a fancy wormhole graphic...just like I presume inter-system travel loading screens are handled.

Can't find it.

1st option: linear interpolation. rk4 needs to give you the rotation <angle,axis>that gets you from the current state to the future state. If you just have the current state <angle,axis> and future state <angle,axis>, you could, in principle extract the rotation that got you there, but it is computationally expensive and complicated. It is more efficient to just carry it around in the state. So, the rotation <angle,axis> object is given. The value of angle is just the average angular velocity (av) times the time step(ts). So, I write interp_angle=(angle/ts)*t where t is the elapsed time since current state. You can see that when t=ts, we are at the future state and the value of angle is what it should be. Then, I construct the interp_rot_quaternion=quaternion(cos(interp_angle/2),axis*sin(interp_angle/x)). I then say interp_quaternion=interp_rot_quaternion*current_state_orientation_quaternion. This gives you the orientation at a time in between the current state and future state. And, you are guaranteed to end up where you are supposed to at t=ts.

It isn't pseudo-code, but that is the algorithm.

The gravitation of asteroids is too feeble to have any effect perceivable to the player.

Well, yes and no. Ceres is a large asteroid with .27m/s^2 acceleration at its surface, but there are few of those that exert significant gravity at any large distance. That is how I sized the spheres of influence.

a=GM/r^2--> r=sqrt(GM)/sqrt(a)...in my sphere of influence calculations before, I took my 'a' threshold to be the sun's gravitational pull at 100 AU which is around 6*10^-7 m/s^2. If you make the gravitational threshold higher, you will need to make exceptions for the sun and probably some of the gas giants with lots of moons to make sure that these objects are influencing the objects that are bound to them gravitationally. Otherwise, pluto just goes floating off into space. I also chose the low threshold so that even if something is traveling very fast, it will likely get picked up by the sphere of influence well before the gravitational pull is noticeable to the player.

Also, without micro-gravity, it will be impossible to keep the asteroid field from dispersing. The asteroids will all have small random velocities in the radial direction. In reality, the gravitational pull of the asteroids in the belt pulls back any errant asteroids that start to wander out. Without this, we will need to create some sort of artificial acceleration that keeps these objects in a belt.

While the effect of gravity on the player may be negligible. The cumulative effects of gravity on objects in the system(or the lack thereof IS perceivable to the player). If the asteroid belt disperses, or some of Jupiter's moons going flying off...or Saturn's rings disperse, these things will be noticed.

Another thing is that pluto holds charon in orbit via micro-gravity. It may be appropriate to have each object in the system have a list of "attractors", but have them be a fixed list for most objects. Only players would have dynamic lists. Objects like asteroids that could be moved from their normal orbits by player intervention may also become partially dynamic once they have been removed with their fixed list continuing to influence them all the time, but having the option to add other influences temporarily if it wanders into the sphere of influence of another object.

Damn! Tensor of inertia is a matrix, not a vector. Couldn't it be expressed more succinctly as three quantities?, --namely, the moment of inertia around each axis? I would envision us reading just 3 quantities from units.csv, then expanding that to a matrix, if necessary (in RigidBody's constructor); and maybe the inverse tensor of inertia as well, if it speeds up calculations.

It is not as simple as a moment of inertia about each axis. I can't really explain it that well. It is best seen through demonstration, but first, we need to calculate the moi tensor of an actual non-symmetrical object. Let me try and find a suitable object and get back to you. You will definitely want to carry around the inverse tensor of inertia. In fact, that is really the only thing that we need. We don't care about the actual moment of inertia. We are always doing torque*inverse_moi. Maybe for purposes of computing changes in the moment of inertia tensor, we will need it, but the computational expense of inverting matrices is such that I would definitely recommend carrying around the inverse.