Physics simulations in python

[wumpus]

On 11/7/2018 at 7:32 PM, K^2 said:

*shrug* No argument in general, but the smallest servos I've worked with just happen to be 9g each as well, and you'll want at least two of these, plus a power pack to keep them running. So you're probably looking at something like 30g and $20+ of hardware before you consider MCU. That still doesn't make cost and weight of Pi Zero trivial, but certainly blunts the impact of going with something smaller.

Anyways, without more details, "It depends," is probably the most that can be said, yeah.

I suspect that the actuators (or whatever actually moves the fins) and the interface it needs will determine what computer it needs.  Hobbyist computers may seem low-volume, but nothing compared to computer-controlled model rocket fins.  Presumably the thing will be based on RC-aircraft parts, and check there before solidifying the rest of the electronics.

[Gargamel]

11 hours ago, wumpus said:

I suspect that the actuators (or whatever actually moves the fins) and the interface it needs will determine what computer it needs.  Hobbyist computers may seem low-volume, but nothing compared to computer-controlled model rocket fins.  Presumably the thing will be based on RC-aircraft parts, and check there before solidifying the rest of the electronics.

You'd be surprised, the software libraries that drive the servos controls are very flexible for micro controllers, but most servos run on industry standards for communication.    They have to, as the parts are all modular.  The transmitter sends a signal out, which is on the 2.4ghz spectrum (which is the modern standard used for RC), and is encoded to only be interpreted by the specific receiver it is paired with.  That way you can have many many transmitters working on the same frequency (No, I don't know the details of how it works).   The receiver then interprets the data sent to it, and sends out the appropriate signals via the various ports it has.  You can then stick pretty much any servo you want into any of the ports on a reciever.  The biggest issue is if you have the right amount of power being supplied to the servos to run them.    You really have to step into industrial application servos before you start getting compatibility issues, and even then, a simple software library on the micro controller should fix that problem.  You just can't use industrial style servos on standard off the shelf RC systems.   TLDR: So as long as you have the right power supplied for what the servo demands, a programmable micro controller will be able to run just about any servo with the right library. 

I use the same/similar servos on my RC racing yachts as I do on the electronics projects I mess around with.   We've gotten real fancy and slaved multiple boats to one transmitter, it was fun to watch them sail in formation, but it never ended well, LOL. 

[Mad Rocket Scientist]

24 minutes ago, Gargamel said:

[...]

I use the same/similar servos on my RC racing yachts as I do on the electronics projects I mess around with.   We've gotten real fancy and slaved multiple boats to one transmitter, it was fun to watch them sail in formation, but it never ended well, LOL. 

Just got nerd-sniped trying to think about how to make that work.

I personally don't have much experience with RC servos, but my grandfather flew RC planes and I was always amazed how well they worked. They never needed configuration or calibration and I have never seen one fail. Not only that, modding them for continuous rotation is (relatively) easy.

On making a smooth version of the thrust curves:

This is the kind of data I have (actually, there are 4 formats, but this is the basics):

   0.0070 23.0
   0.018 25.0
   0.027 20.25
   0.066 20.25
   0.073 18.5
   0.094 20.25
   0.112 20.75
   0.137 19.75
   0.163 21.5
   0.202 20.75
   0.231 20.75
   0.254 22.75
   0.27 20.75
   0.504 20.0
   0.536 18.25
   0.607 17.0
   0.687 14.75
   0.751 14.25
   0.84 11.25
   0.998 8.25
   1.024 8.25
   1.248 2.5
   1.385 0.0

The number on the left is time, the one on the right is thrust. This data has relatively few points.

I've had very little trouble turning this into a bunch of lists, but making a curve function for it is where the challenge is. This is what the data looks like when graphed:

First, I pursued bezier curves, with each point as a control point for the bezier curve:

As you can see, it smooths the curve, but loses much of the detail. More importantly, I'm not sure if the bezier package I'm using actually supports n-degree curves, since it seems to look that same whether I input any manual degree or automatically calculate the number of degrees from the number of control points (1 per control point). It also seems to break with more than about 50 control points.

Then I want for cubic splines, making the function pass through all the points on the original data (here seen overlaid on the original data):

This is the most successful, since it preserves all of the original points, and provides an easy to use function for all points on the spline. However, for data that has large areas between points or is noisy, it produces points that are quite far away. Many of the other data sets work much better than this one.

In the long run, probably none of this matters too much since this data probably won't match what we actually get from the engines, but should I continue to try to refine the interpolation?

[Kerbart]

2 hours ago, Mad Rocket Scientist said:

(...)

I've had very little trouble turning this into a bunch of lists, but making a curve function for it is where the challenge is. This is what the data looks like when graphed: (...)

Not to derail the discussion, but there's a savefig(<filename>) method on your chart plots (assuming you're usig matplotlib or something similar; it certainly has the look of it). Spares you the effort of saving screenshots and they look prettier here, too. You can even save them as SVG files.

Edited November 11, 2018 by Kerbart

[K^2]

 

5 hours ago, Gargamel said:

You'd be surprised, the software libraries that drive the servos controls are very flexible for micro controllers, but most servos run on industry standards for communication.

Absolutely. Some of the modern digital servos can get a little fancy, but Pulse Width Modulation, which is standard for analog servos, is almost universally supported even by these. PWM is almost trivial to produce, as you just have to turn output pin on and off with the right timing. Duration of the pulse determines the angle the servo will hold. That's it. No fancy code or hardware drivers necessary.

Edited November 11, 2018 by K^2

[Gargamel]

12 hours ago, K^2 said:

PWM is almost trivial to produce, as you just have to turn output pin on and off with the right timing.

Most arduino based chips already have at least a couple PWM pins native to the board/chip.

16 hours ago, Mad Rocket Scientist said:

and I have never seen one fail.

Oh trust me, they do fail.   Our main sail servos have to hold 20lbs of force, with some wind gusts putting 50lbs on them momentarily.  They do burn out eventually.    Whats worse is when your rudder servo fails, and it just spins in the middle of the lake for hours. 

13 hours ago, Kerbart said:

Not to derail the discussion,

Don't worry, this discussion already has more tangents than a geometry class. 

[wumpus]

On 11/10/2018 at 8:19 PM, Gargamel said:

You'd be surprised, the software libraries that drive the servos controls are very flexible for micro controllers, but most servos run on industry standards for communication.    They have to, as the parts are all modular.

The other issue is the accelerometer (or whatever you are using to determine landing).  You probably also need some sort of laser/radar distance measurement to the ground.  Depending on the TWR of the rocket, these will likely have some hefty speed/latency issues: try to make sure that I²C and/or SPI is good enough for both of these before making your final selection (it might be worth it to get a ardino or raspberry pi that includes an accelerometer on some sort of fast internal bus).

Looking at the ATTiny 8 pin chip, it looks like pins 2&3 are tied to the crystal/ceramic resonator.  That leaves one set of x-directional fins, one set of y directional fins, and one input pin (which might be an I2C bus to both the accelerometer and rangefinder, that just sounds too slow for a hoverslam).  I'd recommend at least stretching it out to a 16 bit package to use the included extra two ports for something.

Once the physics simulation is done in python (remember the point of the thread?), you should know just how much latency you can eat on your accelerometer and range finder.  That will likely drive the selection of  those parts, an in turn find a microcontroller than can interface with them (to a sufficiently low latency).  I'd expect that a ATTiny can react faster than an actuator can move, but it needs the input signals fast enough as well.

[Gargamel]

9 hours ago, wumpus said:

I'd recommend at least stretching it out to a 16 bit package to use the included extra two ports for something.

Yes, the 16 pin ATTiny isn't that much bigger, and I don't think I can tell the difference in mass in my hand.   Obviously there is one, but it's pretty negligble for our needs. 

9 hours ago, wumpus said:

Depending on the TWR of the rocket, these will likely have some hefty speed/latency issues:

The accelerometers I've used for Arduino don't have much latency from what I can tell.   I think the one I was playing with was sampling around 100hz, and was responsive enough to make a steady cam. 

As for the range finding, the ultrasonic sensors I mentioned might have enough latency to cause an issue, but they should be able to source one that is fast enough.

9 hours ago, wumpus said:

That leaves one set of x-directional fins, one set of y directional fins, and one input pin (which might be an I2C bus to both the accelerometer and rangefinder, that just sounds too slow for a hoverslam).

I've never played with I2C, but given how many components we're using, I don't think it's necessary, and most everything can be wired directly to the chip. 

Redacted comment about fins, cause I'm dumb. 

[K^2]

I can think of a bunch of potential problems with using a sonar on a rocket. There are optical range-finders with 3m+ range. That ought to be better. Lower latency, fewer problems with noise, etc.

[Mad Rocket Scientist]

If I model the rocket engine as an object with a location relative to the rocket's origin, where does the CoT go? At the origin of the engine object? Or does it really matter?

EDIT: How should I calculate the CoT from the position and rotation of multiple engines?

Edited November 19, 2018 by Mad Rocket Scientist

[K^2]

1 hour ago, Mad Rocket Scientist said:

If I model the rocket engine as an object with a location relative to the rocket's origin, where does the CoT go? At the origin of the engine object? Or does it really matter?

It matters a little. Possibly not enough, but you might as well try to get it as close as possible. For a real rocket engine, this is a complicated question, because the pressure from exhaust is distributed between combustion chamber and nozzle bell. Fortunately, most hobby rocket engines don't produce high enough exhaust velocity to warrant a bell, so you can basically ignore the nozzle. The CoT for these should be at the top of the combustion cavity. That's not going to be quite the top of the motor, since even if you're using one without ejection charge, there will still be a plug at the top. It's going to vary from motor to motor, but about a centimeter down from the top of the plug is going to put you in the ball park. CoT should, ideally, be centered on the axis of symmetry. But this is another thing you might want to add a slight variation to in your models just to make sure that the real world doesn't ruin your model. After all, hobby motors aren't made to that high of a spec. You might want to vary direction of thrust very slightly as well for the same reason.

1 hour ago, Mad Rocket Scientist said:

How should I calculate the CoT from the position and rotation of multiple engines?

You probably don't need to model an overall effective CoT. You could, but it's probably going to be easier just to apply forces and torques from each one individually.

To compute net force, you don't need to worry about CoT at all. You just add thrust vectors from all engines and you have net thrust. For net torque, compute torque vectors from each engine individually, then add them together just like you did with thrust. For each engine, you need to compute the arm, which is CoT coordinate minus rocket's CoM coordinate, call that r. Then torque = r x F for each engine.

In terms of recomputing CoT positions as the rocket rotates, you can either apply the transformation matrix to each arm before multiplying, or you can use the fact that net torque transforms the same way. So you can do your torque computations in the rocket's coordinate system, and then rotate it to world coordinates using rocket's rotation matrix. The angular momentum EoM: dL/dt = net torque, works in world coordinates, so you definitely want to have net torque in world coordinates before you do angular momentum update. Once you have updated angular momentum, you need to compute angular velocity by using rocket's moment of inertia tensor, which is in rocket's coordinate system! So you'll have to either rotate angular momentum to rocket's or rotate moment of inertia tensor to world (keeping in mind that you can rotate principal axes like you do any other vector) and use angular velocity in world coordinates to update rotation matrix.

Yeah, the rotation bit is quite convoluted, but if you can keep it straight which quantities are preserved in which coordinate system, you should be fine.

P.S. If trying to derive Verlet or RK for rotation is too much, you can probably get away with just using Forward Euler method for the rotation part. I've done that in some simulations before, and I can't say I've noticed much difference. Fortunately, most of the types of situations which make FE bad for CoM update just don't arise with rotations all that often.

[Mad Rocket Scientist]

On 11/18/2018 at 7:06 PM, K^2 said:

It matters a little. Possibly not enough, but you might as well try to get it as close as possible. For a real rocket engine, this is a complicated question, because the pressure from exhaust is distributed between combustion chamber and nozzle bell. Fortunately, most hobby rocket engines don't produce high enough exhaust velocity to warrant a bell, so you can basically ignore the nozzle. The CoT for these should be at the top of the combustion cavity. That's not going to be quite the top of the motor, since even if you're using one without ejection charge, there will still be a plug at the top. It's going to vary from motor to motor, but about a centimeter down from the top of the plug is going to put you in the ball park. CoT should, ideally, be centered on the axis of symmetry. But this is another thing you might want to add a slight variation to in your models just to make sure that the real world doesn't ruin your model. After all, hobby motors aren't made to that high of a spec. You might want to vary direction of thrust very slightly as well for the same reason.

That makes sense.

Quote

You probably don't need to model an overall effective CoT. You could, but it's probably going to be easier just to apply forces and torques from each one individually.

I'm not sure whether to be relieved that that is the easier solution or concerned that it is.

Quote

To compute net force, you don't need to worry about CoT at all. You just add thrust vectors from all engines and you have net thrust. For net torque, compute torque vectors from each engine individually, then add them together just like you did with thrust. For each engine, you need to compute the arm, which is CoT coordinate minus rocket's CoM coordinate, call that r. Then torque = r x F for each engine.

In terms of recomputing CoT positions as the rocket rotates, you can either apply the transformation matrix to each arm before multiplying, or you can use the fact that net torque transforms the same way. So you can do your torque computations in the rocket's coordinate system, and then rotate it to world coordinates using rocket's rotation matrix. The angular momentum EoM: dL/dt = net torque, works in world coordinates, so you definitely want to have net torque in world coordinates before you do angular momentum update. Once you have updated angular momentum, you need to compute angular velocity by using rocket's moment of inertia tensor, which is in rocket's coordinate system! So you'll have to either rotate angular momentum to rocket's or rotate moment of inertia tensor to world (keeping in mind that you can rotate principal axes like you do any other vector) and use angular velocity in world coordinates to update rotation matrix.

Yeah, the rotation bit is quite convoluted, but if you can keep it straight which quantities are preserved in which coordinate system, you should be fine.

Right, so I'm essentially dealing with the translational and rotational forces separately.

Quote

P.S. If trying to derive Verlet or RK for rotation is too much, you can probably get away with just using Forward Euler method for the rotation part. I've done that in some simulations before, and I can't say I've noticed much difference. Fortunately, most of the types of situations which make FE bad for CoM update just don't arise with rotations all that often.

My calculus is too weak to even derive the velocity Verlet, so I'll probably do the same thing I did for that: look up how it is used in simulations and try to understand the formula and code itself, not the derivation.

How should I deal with local rocket coordinate systems? The simulation takes in the center of mass of the entire vehicle, but the CoM's relative position within the rocket changes throughout the flight. I was thinking of defining the locations of all the parts relative to the arbitrary origin, then computing the relative position from that to the CoM, and moving the origin by that much every time step. So the CoM is globally simulated, then the origin is defined relative to the CoM, then all the parts are defined relative to the origin. But is that just over complicating things?

Edited November 20, 2018 by Mad Rocket Scientist

[K^2]

You can have the definitions of rocket parts given relative to arbitrary origin and include CoM relative to that origin. But for run-time, I would move all relevant coordinates to be relative to CoM, so that CoM is your origin. Having CoM at (0, 0, 0) makes a lot of math WAY easier. For starters, locations of your engine CoTs are your moment arms for torque.

If you do this, the entire state of your rocket is given by these 6 vectors, all in world coordinates:

1. Orientation of Rocket's X axis.
2. Orientation of Rocket's Y axis.
3. Orientation of Rocket's Z axis.
4. CoM location.
5. CoM velocity.
6. Angular momentum.

If you look up how transformation matrices are done in video games (there tends to be more good literature on that than simulations), 1.-4. basically give you the rows of the matrix that transforms vectors in rocket's coordinate system to world coordinate system. (The last column of the 4x4 matrix is always going to be [0, 0, 0, 1], and is sometimes omitted from being recorded for that reason.) Alternatively, you can split translation and rotation. In which case, 1.-3. give you the 3 rows of the 3x3 rotation matrix. Since you aren't writing for performance, whichever way makes more sense to you is good enough.

The update rules for the 3 axes are pretty straight forward. dX/dt = ω x X, and same for Y and Z. Of course, numerical errors will cause norms and relative orientations to drift. So after each step it's useful to renormalize all 3 axes and ensure that they are mutually orthogonal. (E.g, you could just set Z = X x Y, then X = Y x Z. There are cleaner ways, but this will probably do.)

The only thing that's missing is ability to transform back from world coordinate system back to rocket's. And you do that by inverting the transform matrix. In general, inverting a 4x4 matrix is a bit of work, but transform matrices are special. The rotational part is an SO(3) matrix, and they are Unitary. So to compute inverse of the rotation 3x3 matrix, just transpose it! The translation bit, 4th row, is only slightly trickier. You negate it, and then using the inverse rotation you've just derived, transform it from world space to rocket's space.

And just to save you a bit of confusion, I've used DirectX convention of writing transformation as x' = x M, where M is transformation matrix, x is vector before transformation ans x' is the same vector after. If you prefer to multiply vectors from the right, so that x' = M x, then switch the words "row' and 'column' everywhere in the above text, and the math will still work.

[Mad Rocket Scientist]

Alright, so school is over and I have some more time to work on this again. Currently I'm using this code for the position (velocity verlet):

rocket.position += rocket.velocity * dt + 0.5 * rocket.acceleration * dt * dt

But then I'm using forward Euler for velocity:

rocket.velocity += rocket.acceleration * dt 

This PDF on slide 53 shows another way of calculating velocity. Should I try to use that (I'm not sure how to compute F_{t + dt} from position_{t + dt} and velocity_{t + dt}, I'm currently calculation force as the sum of thrusts divided by mass) or is Euler close enough? The slide uses h where I'd use dt, X where I'd use position, and V where I'd use velocity, FYI.

And I'm afraid I'm going to have to cry uncle on trying to figure out Verlet for rotation. Does anyone have any sources for what Euler would look like for rotation? I think I understand most of it (torque ~ acceleration, angular momentum ~ velocity, and the three unit orientation vectors which make up the rotation matrix ~ position), however the way torque affects, say, a tube is different dependent on what direction it is applied in. I think that difference can be modeled with the moment of inertia tensor, but I'm not sure how to make that, or what it would even look like. Would it be a three dimensional array in 3D? 

[K^2]

23 minutes ago, Mad Rocket Scientist said:

Should I try to use that (I'm not sure how to compute F_{t + dt} from position_{t + dt} and velocity_{t + dt}, I'm currently calculation force as the sum of thrusts divided by mass) or is Euler close enough? The slide uses h where I'd use dt, X where I'd use position, and V where I'd use velocity, FYI.

Doing velocity portion properly is critical when you have position-dependent (or velocity-dependent) forces. For example, if you were building a KSP clone, and you wanted to add joints, which are simulated as springs, the forces on joints would depend on relative positions of your ship parts. And if you didn't properly account for F_t and F_t+dt, your ship will get Krakened almost instantly.

In your case, your main forces are gravity, which is constant, and thrust, which depends on time, but not position or velocity. So the worst you'll have is thrust effectively changing half a tick out of sync with your simulation, which is such a minor effect, that you really don't need to bother.

For future reference, if you did have forces that depend on position, there is a super simple hack to account for it in Velocity Verlet.

# Compute acceleration at time t
force_t = compute_force(time, position) # Compute F_t
accel_t = force_t / mass # Use rocket's current mass to compute accel.

# Update position and time.
position = position + velocity * dt + 0.5 * accel * dt * dt # Update position. position is now position_t+dt!!!
time = time + dt
mass = mass - compute_mass_flow(time) * dt # Simple forward-Euler on mass, if you're even keepint track of it.

# Now compute acceleration at time t+dt
force_tpdt = compute_force(time, position) # Compute F_t+dt, because time is now t+dt and position is position_t+dt
accel_tpdt = force_tpdt / mass # Mass has been updated as well, so this is accel_t+dt

# Update velocity
velocity = velocity + 0.5 * (accel_t + accel_tpdt)

Note that on the next iteration of the loop, accel_t will evaluate to exactly the same value as accel_tpdt from the previous iteration of the loop. So as long as you initialize accel_tpdt sensibly outside of the loop, then inside you can replace first two lines with accel_t = accel_tpdt.

 

Where things get worse is when your force depends on velocity. And technically, drag does. So your code will have compute_force(time, position, velocity), and because of that, the Velocity Verlet becomes an implicit method. So a little bit of jargon. Explicit methods are the ones where previous steps of integration provide you with all the variables you need to know to compute the next step. So in the loop above, we are using position_t and position_t+dt, but by the time we need the later, we have it computed. So Velocity Verlet is explicit when forces depend only on position and time. The method is called implicit if it relies on information from a future state, one that you're yet to compute, in order to advance. Because we need accel_t+dt to compute velocity_t+dt, if force_t+dt depends on velocity_t+dt, you got yourself a circular dependence.

The proper way of getting about it is by using an iterative methods that solves for accel_t+dt and velocity_t+dt simultaneously. For drag, it'd probably be sufficient to simply take current velocity to compute accel_t+dt, and use that to compute velocity_t+dt, then use this as your velocity to compute accel_t+dt again, and update your guess for velocity_t+dt, and so on until the two values stop changing (converge). It shouldn't take too many iterations.

That said, drag rarely changes in a way that produces resonance. You'd have to have some model that produces realistic turbulence to have that, and if you're dealing with that, it's a whole another class of problem. I would honestly not bother. Just use the current value for velocity when running force computations, and treat it as an explicit method.

As a general rule, if you're forced to use an implicit method, odds are, you can just reduce step size, and use an explicit method, and still get better precision out of fewer computations. The only hard exception I know are central potential problems, such as gravity. If you're trying to predict a path of a comet through Solar system with enough precision to distinguish between close fly-by and a planetary impact, you have to use implicit methods. For basically everything else, you should be able to scrape by with explicit methods.

23 minutes ago, Mad Rocket Scientist said:

And I'm afraid I'm going to have to cry uncle on trying to figure out Verlet for rotation. Does anyone have any sources for what Euler would look like for rotation? I think I understand most of it (torque ~ acceleration, angular momentum ~ velocity, and the three unit orientation vectors which make up the rotation matrix ~ position), however the way torque affects, say, a tube is different dependent on what direction it is applied in. I think that difference can be modeled with the moment of inertia tensor, but I'm not sure how to make that, or what it would even look like. Would it be a three dimensional array in 3D? 

Moment of inertia tensor is a 3x3 matrix. That matrix is diagonal in coordinate system spanned by its principal axes. For a rocket with typical symmetries, the minor axis will run along the center line of the rocket (so directly up when rocket stands on the ground). If your rocket has rotational symmetry around that axis with any angle less than 180° (So if you have 3 identical fins, 120°, or 4 identical fins, 90°) then the other two axes will be degenerate, and you can point them wherever you like. Along rocket's X and Y is usually convenient.

So if you set your coordinate system up like this, then you don't really need to track the full matrix. You'll simply have the angular_velocity_x = angular_monentum_x / inertia_x, etc for y, z. But this only works in coordinate system of the rocket. This is not the X, Y, Z of your world. So your update loop might look something like this:

# First, compute torques, just like with forces.
torque_t_w = compute_torques(time, rot_matrix, ...) # The parameter list might get a little crazy if you arne't doing this with objects.
torque_t_l = vector_to_local(torque_t_w, rot_matrix) # This function will take world coordinates vector and transform it to local (rocket) coordinates vector.

# Now we need to compute update to orientation (rot_matrix) of the rocket.
ang_momentum_l = vector_to_local(ang_momentum, rot_matrix)
ang_velocity_l = ang_momentum_l / inertia # This should be a component-wise division. If you aren't using numpy, might need to be a loop.
ang_delta_l = ang_velocity_l * dt + 0.5 * dt * dt * torque_t_l / inertia # Here, we use torques by analogy to acceleration in position update.
ang_delta_w = vector_to_world(ang_delta_w, rot_matrix)

# Finally, we have to apply this and udpate time.
rot_matrix = rot_matrix * rodriguez(ang_delta_w) # Double-check if you should be multiplying from left or right. rodriguez() is Rodriguez rotation matrix.
time = time + dt

# And now we do update for the things that happen at t+dt
torque_tpdt_w = compute_torques(time, rot_matrix, ...)
torque_tpdt_l = vector_to_local(torque_tpdt_w, rot_matrix)

ang_momentum_l = ang_momentum_l + 0.5 * (torque_t_l + torque_tpdt_l)
ang_momentum_w = vector_to_world(angl_momentum_l, rot_matrix)

Ok, so two things to note. The primary quantity we're keeping track of is angular momentum, not angular velocity. Because the former is conserved, ant the later isn't.

Second, note that I'm updating time = time + dt in the middle. This is because this block of code isn't meant to be separate from the position update above. The two are meant to work together. In fact, you'll probably want to compute torques and forces in one pass both at the top and at the bottom of the loop. Mostly, the order isn't critical, but you need to make sure that everything that happens before the time update happens before, and everything after happens after. 

I've also shorthanded a LOT of matrix math. I hope you can follow it, but that's going to be a big chunk of your effort. The one thing I want to call out specifically is that I'm treating inertia as a vector, rather than matrix. That's because everything's in principal axes coordinates, and the only numbers I care about are on the diagonal of that matrix (everything else is zero in this CS). So we can store it as a vector and not bother with the full tensor. Again, so long as your rocket's local Z runs along the symmetry axis.

[Mad Rocket Scientist]

Thanks, these posts are really helpful.

Edited December 17, 2018 by Mad Rocket Scientist