Correct physics constants - center of mass and moments of inertia.

[Swamp-Ig]

It would be better if the correct physics constants were used for parts as created by procedural parts - including the moment of intertia and the center of mass.

This bug is to track development.

[Swamp-Ig]

For the bezier cones I've used pappus's theorem to work out the volume

V = 2pi * Area of lamina * geometric centroid y axis geometric centroid y axis = Moment about y / Area of lamina

so:

V = 2pi * moment about y

moment about y = integrate y^2 dx dy (which is a surface integral)

Can then use green's theorem to turn it into a line integral - equation 13:

= 1/2 integral_0^1 x^2 y' dt

I know what x(t) and y(t) are from the definition of the bezier, so can do that calculus. I'll spare you the workings but it works out as:

float M_y = 1f/4f* ((p1.y-p0.y)*(1f/3f * p0.x*p0.x + 1f/4f * p0.x*p1.x + 1f/14f* p0.x*p2.x + 3f/28f* p1.x*p1.x+ 3f/28f * p1.x*p2.x + 1f/84f * p0.x*p3.x + 3f/70f * p2.x*p2.x + 1f/35f* p1.x*p3.x + 1f/28f* p2.x*p3.x + 1f/84f* p3.x*p3.x) +
                    (p2.y-p1.y)*(1f/12f* p0.x*p0.x + 1f/7f * p0.x*p1.x + 1f/14f* p0.x*p2.x + 3f/28f* p1.x*p1.x+ 6f/35f * p1.x*p2.x + 2f/105f* p0.x*p3.x + 3f/28f * p2.x*p2.x + 1f/14f* p1.x*p3.x + 1f/7f * p2.x*p3.x + 1f/12f* p3.x*p3.x) +
                    (p3.y-p2.y)*(1f/84f* p0.x*p0.x + 1f/28f* p0.x*p1.x + 1f/35f* p0.x*p2.x + 3f/70f* p1.x*p1.x+ 3f/28f * p1.x*p2.x + 1f/84f * p0.x*p3.x + 3f/28f * p2.x*p2.x + 1f/14f* p1.x*p3.x + 1f/4f * p2.x*p3.x + 1f/3f * p3.x*p3.x));

where p0, p1, p2, and p3 are the control points of the bezier. (In the highly unlikely event you do the maths to check it, I'm aware that there's a factor of 1/2 included).

But you can't do the same trick for the moment of intertia, because there's no pappus centroid theorem.

Thinking about it, you can work this out with some calculus. You have some function that defines the boundary of the volume that will be revolved:

f(s) = ( f_x(s) , f_y(s) ) where s = 0..1

so: x = f_x(s) y = f_y(s)

Pick the origin as the point to compute the moment about. I'll switch now to using D for complete and d for partial integrals

so for each point on the curve, each Ds, you have a thin cylinder that will be revolved around the y axis.

We know the formula for moments of inertia of cylinders

Ix = Iz = 1/12 m ( 3 r^2 + h^2 ) Iy = m ( r^2 / 2 + y^2 ) (using the parallel axis theorem to find about the origin )

We need to do some changes of variables.

So r = x obviously.

Now... This is a bit I'm not so sure about:

h = df_y/ds = y'

You can use the partial derivative here since you're only interested in the variance with changing y, right? If not, how can I do the change of variables?

m = rho pi r^2 h = rho pi x^2 y' (rho = density)

Then integrate running your point around the curve so:

Ix = integrate_0^1 1/12 rho pi x^2 y' ( 3 x^2 + y'^2 ) Ds = rho pi/4 integrate_0^1 x^4 y' + 4 x^2 y'^3 Ds

Iz = integrate_0^1 rho pi x^2 y' ( x^2 / 2 + y^2 ) Ds = rho pi/2 integrate_0^1 x^4 y' + 2 x^2 y^2 y' Ds

Now this makes sense, because of the dependence on y' and all odd powers then if the curve starts running down, you get a negative value and subtract off from your moment.

I don't want to go ahead and do the integral with all those forth powers until I feel the maths is right. If you can tell me if the bold section is right then I'm happy with the rest of it.

[eggrobin]

Ok, so I haven't read the math in detail, it's 01:30 here and I'm a tad tired, but if it's just a 4-point Bézier (I somehow assumed it was a spline, which is more of the same, but makes things somewhat tedious) which does not go back it should be completely trivial, since then the function f(x) is just a single cubic polynomial (the Bézier form gives you additional information as the parametrisation of its graph, but you can obviously go back to f; a Bézier curve is just a pair of cubics (x(t), y(t)), if it does not go back---x is monotonous---, you can change variables to (x, y(x))). In that case it's just a question of writing f down and applying the formulae you need.

[Swamp-Ig]

You found my bug obviously. Get some sleep :)

Beziers aren't that easy to change so that y = f(x). Besides I'm feeling like I really want to solve the original problem now. :)

If you can take a look at that one highlighted line and ensure it is correct, then the rest of it falls into place.

[eggrobin]

Bézier curves are easy to turn into a cubic under the assumptions you make. Let P_k = (ξ_k, η_k) for k = 0, 1, 2, 3 be the control points of your cubic Bézier. the tangents are m₀ = (η₁ - η₀) / (ξ₁ - ξ₀) and m₁ = (η₂ - η₃) / (ξ₂ - ξ₃).

The cubic is the unique cubic through P₀ and P₃ with those tangents. Let x₀ = ξ₀, x₁ = ξ₃, p₀ = η₀, p₁ = η₃, Wikipedia gives you the formula you want (p(x) = f(x)).

EDIT: I gave you the wrong link, fixed. EDIT2: Fixed wrong notation.

[eggrobin]

Actually, representing your curves as (p₀, m₀, p₁, m₁) might make more sense, since that exactly embodies subset of splines you allow.

[Swamp-Ig]

You're right, this is actually a better approach. I wasn't aware of Cubic Hermite splines until you pointed them out to me.

I may well reinterpret this using a monotone cubic approach, for a more general case when there's a bunch of points.

[eggrobin]

Have you thought about the convex collider issue? In order to force your shape to be convex you want the function to be concave, i.e., m₀ > m₁.

I'm not sure you should care about monotony (and actually that would disallow stuff like ellipsoids), but maybe there's something I haven't thought of.

[Swamp-Ig]

Yeh I know about the convex collider issue.

But... is my above maths right? Can you just use the partial, rather than the complete integral?

Using the same method seems to work for mass (ends up with the same formula as far above)

[eggrobin]

I still haven't looked very carefully at it, but the bold line looks fishy. It helps to take a physical look at things here and to consider the units (you can formalise that, and I use it in Principia, but let's not get into that): Let's say our parametrisation is as a function of time,

h          = df_y/ds             = y'
[Distance] = [Distance] / [Time] = [Speed]

[eggrobin]

The integrals seem correct though. You can formalise things differently, and avoid having to reason with awkward infinitesimals, but that leads to a lot of useless abstraction.

[Swamp-Ig]

Go get some sleep :smile:

I think h is actually dy, (speed) because you're deliberately integrating infinitesimal slices.

[Swamp-Ig]

You can ignore the dx part of it. Ideally the thin cylinders would be thin truncated cones, with base radius of x and top radius of x + dx. But since dx is infinitesimal you end up with the cylinder case anyhow.

... and height dy obviously.

[eggrobin]

Well, actually, h = dy is a distance, but dy = df_y/ds * ds, which is what you integrate (and that's where this notation gets so tautological it's confusing, unless you really formalise nonstandard analysis which is an interesting rabbit hole; differentials are the way to go). Note that when your curve is just the cubic in x described above, this simplifies to s = x, df_y/ds = df/dx = f'.

[Swamp-Ig]

dy = ∂f_y/∂s * ds

... does it?

I thought it was:

dy = df/ds * ds = (∂f_x/∂s dx + ∂f_y/∂s dy) ds

which really doesn't help.

although really it's h=∂y isn't it? (since we're only interested in the variance in Y)

I'm really confused :confused:

---

Anyhow, assuming you are correct.

m(s) = ρ π x² dy

I_x(s) = I_z(s) = 1/12 m(s) ( 3 x² + dy² ) = 1/4 ρ π ( x⁴ dy + 4 x² dy³ ) = 1/4 ρ π ( x⁴ y' ds + 4 x² (y' ds)³ )

I don't know how to handle the (y' ds)³. Can't one just neglect the third order differential?

I_x = 1/4 ρ π integrate x⁴ y' ds s=0..1

I_y(s) = m ( x² / 2 + y² ) I_y = ρ π/2 integral x⁴y' + 2 x² y² y' ds

---

Incidentally I don't think the bernstein polynomial above is unique. Your tangent functions are scalars, but there's a magnitude involved:

[eggrobin]

"dy = df/ds * ds = (∂f_x/∂s dx + ∂f_y/∂s dy) ds" has the wrong units again. Moreover, since f_y is a function R[time] -> R[length] of s, and f is a function R[time] - > R²[length], (univariate) there are no partial derivatives here. We're just taking the y-component of the derivative, dy = df_y/ds * ds. Well, that notation is just confusing. Let's be clearer. The local linear approximation is h := Δy = df_y/ds * Δs = f_y' * Δs

You're right about splines, the inverse of a cubic obviously isn't an affine map, so Hermite splines do not represent all Béziers under your constraints. mea culpa, mea maxima culpa. They represent a subset which is easy to work with, so if you were to switch to splines rather than single segments (which would be a good idea even if you stick to Béziers), you could make your life simpler by using Hermite splines. There would be a loss of generality though and writing checks so that the spline does not 'go backwards' is easy.

EDITED FOR CLARIFICATIONS AND BETTER NOTATIONS

[Swamp-Ig]

Ok cool...

ermmm... does that mean my equations are right?

Incidentally the I_y can be simplified to:

I_y = ρ π/2 integral x⁴ y' ds

since we're only interested in the moment about the center of mass so the second term drops out.

Interestingly, using the same dy = df_y/ds * ds, then you can use the results from this paper, and get the same result, which is vindicating :smile:

[Swamp-Ig]

I checked the integrals on spheres and tori - seems to get the right result. :smile:

Going to run with it.

Here's the final formulas:

f(s) = ( f_x(s), f_y(s) ) x = f_x(s) , x' = ∂f_x(s)/∂s y = f_y(s) , y' =∂f_y(s)/∂s

M = ρ π integral x² y' ds, s=0..1

ȳ = π integral y x² y' ds, s=0..1 x̄ = z̄ = 0

The moments are about the centre of mass: (0, ȳ, 0)

I_y = ρ π/2 integral x⁴ y' ds, s=0..1 I_x = I_z = I_y / 2

[Swamp-Ig]

Just spent a few hours doing the maths. Here's the results:

a = x₀, b=x₁, c=x₂, d=x₃ e = (y₁-y₀), f = (y₂-y₁), g = (y₃-y₂),

i_y = pi/2 (
  1/5 a^4 e
+ 1/5 d^4 g 
+ 1/35 (6 a^3 b e + a^4 f)
+ 1/35 (6 c d^3 g + d^4 f)
+ 1/455 (54 a^2 b^2 e + 24 a^3 b f + 12 a^3 c e + a^4 g)
+ 1/455 (54 c^2 d^2 g + 24 c d^3 f + 12 b d^3 g + d^4 e)
+ 1/455 (27 a^2 b^2 f + 27 a b^3 e + 27 a^2 b c e + 6 a^3 c f + 3 a^3 b g + a^3 d e)
+ 1/455 (27 c^2 d^2 f + 27 c^3 d g + 27 b c d^2 g + 6 b d^3 f + 3 c d^3 e + a d^3 g)
+ 1/5005 (324 a b^2 c e + 216 a b^3 f + 216 a^2 b c f + 81 b^4 e + 54 a^2 c^2 e + 54 a^2 b^2 g + 36 a^2 b d e + 12 a^3 c g + 8 a^3 d f)
+ 1/5005 (324 b c^2 d g + 216 c^3 d f + 216 b c d^2 f + 81 c^4 g + 54 c^2 d^2 e + 54 b^2 d^2 g + 36 a c d^2 g + 12 b d^3 e + 8 a d^3 f)
+ 1/5005 (324 a b^2 c f + 162 b^3 c e + 162 a b c^2 e + 81 b^4 f + 54 a b^3 g + 54 a^2 c^2 f + 54 a b^2 d e + 54 a^2 b c g + 36 a^2 b d f + 18 a^2 c d e + 2 a^3 d g)
+ 1/5005 (324 b c^2 d f + 162 b c^3 g + 162 b^2 c d g + 81 c^4 f + 54 c^3 d e + 54 b^2 d^2 f + 54 a c^2 d g + 54 b c d^2 e + 36 a c d^2 f + 18 a b d^2 g + 2 a d^3 e)
+ 1/5005 (216 b^3 c f + 216 a b c^2 f + 162 b^2 c^2 e + 108 a b^2 c g + 72 a b c d e + 72 a b^2 d f + 36 a c^3 e + 36 b^3 d e + 27 b^4 g + 24 a^2 c d f + 18 a^2 c^2 g + 12 a^2 b d g + 2 a^2 d^2 e) 
+ 1/5005 (216 b c^3 f + 216 b^2 c d f + 162 b^2 c^2 g + 108 b c^2 d e + 72 a b c d g + 72 a c^2 d f + 36 a c^3 g + 36 b^3 d g + 27 c^4 e + 24 a b d^2 f + 18 b^2 d^2 e + 12 a c d^2 e + 2 a^2 d^2 g)
+ 1/1430 (81 b^2 c^2 f + 36 a b c d f + 27 b c^3 e + 27 b^3 c g + 27 b^2 c d e + 27 a b c^2 g + 18 a c^3 f + 18 b^3 d f + 9 a c^2 d e + 9 a b^2 d g + 3 a b d^2 e + 3 a^2 c d g + a^2 d^2 f)
)

[eggrobin]

I think it might be useful to recapitulate things, including a complete derivation, since we both said a lot of stupid things (by the way, dimensional analysis on the integrals in the second post shows they make no sense). Moreover, it's the only way to check your result above.

Notation: ℝ_[D] is the real 1-dimensional vector space of real scalar physical quantities with dimensions D, e. g., 1 m ∈ ℝ_[Length]. ℝ_[D]^k is the real k-dimensional vector space of real k-dimensional vector quantities with dimensions D, e. g., (-54 m, 42 m) ∈ ℝ_[Length]². We consider the parametrisation f of the curve γ as a function of time (the choice of time is arbitrary, but it is useful to make that dimensionful to catch errors; time is a pretty intuitive choice). The curve is then: f : ℝ_[Time] ∋ [0, 1 s] → ℝ_[Length]². f = (x, y), where x, y : [0, 1 s] →ℝ_[Length]. γ = f([0, 1 s]).

First and foremost, there are no partial derivatives here, since f : ℝ_[Time] →ℝ_[Length]² is univariate. Please do not write ∂. the derivative of f, f' : ℝ_[Time] →ℝ_[Speed]², is f' = (x', y'). Using Leibniz notation: df/dt = (dx/dt, dy/dt).

We want to compute the moment of inertia of the volume whose boundary is the surface of revolution generated by γ by rotation around the y axis, and which is full on the left side when following the parametrisation f of γ, under the assumption of uniform density ρ.

We will without loss of generality (since our calculation will do implicit prolongation horizontally and on the y axis of an open curve) assume that γ is a closed curve, and that f is injective.

We then simply need to sum small cylinders around the y axis along γ, which will be counted positively when f goes up and negatively when f goes down (thus handling the case of curves that go backwards if you're interested in that. For colliders, some smart partitioning should be done, but it's another problem).

We will use dirty physicist infinitesimals for the sake of brevity.

We consider the cylinder whose side is bounded by (an approximation of) the surface generated by f([t, t + dt]). Using Wikipedia's inertia tensor for a solid cylinder of radius r, height h and mass m (their cylinder is around the z axis, ours is around the y axis) with r = x(t), h = y(t + dt) - y(t) = y' dt, m = π x(t)² y' dt ρ, we get dI_x = dI_z = 1/12 π x(t)² y'(t) dt ρ (3 x(t)² + dt²), dI_y = 1/2 π x(t)² y'(t) dt ρ x(t)². As the infinitesimal is a first-order approximation, dI_x = dI_z = 1/4 π y'(t) dt ρ x(t)⁴, dI_y = 1/2 π y'(t) dt ρ x(t)⁴.

dI_x = dI_z = dI_y / 2, so we'll only look at I_y.

We now integrate that, I_y = y = ∫_{[0, 1 s]} dI_y(t) = 1/2 ∫_{[0, 1 s]} π y'(t) ρ x(t)⁴ dt.

This matches your latest result.

We now consider the special case where f is a cubic Bézier with control points P_k = (x_k, y_k) ∈ ℝ_[Length]². f(t) = (1 - t)³P₀ + 3(1 - t)²t P₁ + 3(1 - t)t²P₂ + t³P₃. x(t) = (1 - t)³x₀ + 3(1 - t)²t x₁ + 3(1 - t)t²x₂ + t³x₃. Let a = y₁ - y₀, b = y₂ - y₁, c = y₃ - y₂, we have y'(t) = 3(1 - t)²a + 6(1 - t)t b + 3t²c. Performing the integration with Mathematica yields a polynomial that equals yours (except you forgot the ρ).

There remains the question of writing it down and calculating it. Note that you do not want to use Math.Pow for integer exponents as it is quite slow. Some experiments with Mathematica yielded the following "simplest" form for the polynomial. It should be both quick and reasonably accurate. Note that I omitted the factor of ρπ/2.

(b*(286*x0*x0*x0*x0+162*x1*x1*x1*x1+162*x2*x2*x2*x2+432*x2*x2*x2*x3+4*x0*x0*x0*(
132*x1+33*x2+4*x3)+18*x1*x1*x1*(24*x2+7*x3)+594*x2*x2*x3*x3+27*x1*x1*(21*x2*x2+1
6*x2*x3+4*x3*x3)+x0*x0*(594*x1*x1+108*x2*x2+48*x2*x3+72*x1*(6*x2+x3)+7*x3*x3)+52
8*x2*x3*x3*x3+2*x0*(216*x1*x1*x1+63*x2*x2*x2+72*x2*x2*x3+36*x1*x1*(9*x2+2*x3)+36
*x2*x3*x3+6*x1*(36*x2*x2+21*x2*x3+4*x3*x3)+8*x3*x3*x3)+12*x1*(36*x2*x2*x2+54*x2*
x2*x3+36*x2*x3*x3+11*x3*x3*x3)+286*x3*x3*x3*x3)+a*(2002*x0*x0*x0*x0+162*x1*x1*x1
*x1+22*x0*x0*x0*(78*x1+12*x2+x3)+36*x1*x1*x1*(9*x2+2*x3)+9*x1*x1*(36*x2*x2+21*x2
*x3+4*x3*x3)+2*x0*x0*(594*x1*x1+9*x1*(33*x2+4*x3)+2*(27*x2*x2+9*x2*x3+x3*x3))+x0
*(594*x1*x1*x1+72*x2*x2*x2+63*x2*x2*x3+108*x1*x1*(6*x2+x3)+24*x2*x3*x3+3*x1*(108
*x2*x2+48*x2*x3+7*x3*x3)+4*x3*x3*x3)+3*x1*(63*x2*x2*x2+72*x2*x2*x3+36*x2*x3*x3+8
*x3*x3*x3)+2*(27*x2*x2*x2*x2+54*x2*x2*x2*x3+54*x2*x2*x3*x3+33*x2*x3*x3*x3+11*x3*
x3*x3*x3))+c*(22*x0*x0*x0*x0+54*x1*x1*x1*x1+9*x1*x1*x1*(21*x2+8*x3)+x0*x0*x0*(66
*x1+4*(6*x2+x3))+108*x1*x1*(3*x2*x2+3*x2*x3+x3*x3)+x0*x0*(108*x1*x1+36*x2*x2+21*
x2*x3+12*x1*(9*x2+2*x3)+4*x3*x3)+x0*(108*x1*x1*x1+72*x2*x2*x2+108*x2*x2*x3+9*x1*
x1*(24*x2+7*x3)+72*x2*x3*x3+9*x1*(21*x2*x2+16*x2*x3+4*x3*x3)+22*x3*x3*x3)+6*x1*(
54*x2*x2*x2+108*x2*x2*x3+99*x2*x3*x3+44*x3*x3*x3)+2*(81*x2*x2*x2*x2+297*x2*x2*x2
*x3+594*x2*x2*x3*x3+858*x2*x3*x3*x3+1001*x3*x3*x3*x3)))/10010

EDIT: Of course, you'll want to declare variables for all those reused powers of x_k.

[Swamp-Ig]

Oh you have mathmatica. Awesome. I'll use your form, since I'll half bet that I've introduced some subtle bug somewhere in the calculation. Incidentally, I've done a slight consistency change to your above notation. The x(t) was written as a function, where as y'(t) was not.

I was also aware of missing out the ρ ,

I've been using alpha, which works well but doesn't allow you very much space to do expansions, so I had to do a lot of factorization by hand.

I was just about to embark on calculating the center of mass, which looks like:

COM = ( x̄ , ȳ ) / M

x̄ = 0 ȳ = π ∫_{[0, 1 s]} y(t) x(t)² y'(t) dt

Do you think you could shove that through mathmatica?

I think I've learned more about calculus in the last few weeks than in an entire year as an undergrad back in 1994 :smile:

[eggrobin]

Mathematica is indeed very often useful for numerical optimisation. It already helped for atmospheric splines in RSS.

Haven't looked at the centre of mass yet, but I have some refinements regarding moments. I fooled with Mathematica a little more, here are the fully optimised calculations (using ExperimentalOptimizeExpression` and some tweaks to make it optimize exponents). The temporary variables have pretty unenlightened names, but they should be legal identifiers. The code is unpalatable in any form, so we might as well overoptimise it. :smile:

// local variables to reduce the number of calculations.
double temp_511 = x0*x0;
double temp_517 = x0*temp_511;
double temp_544 = x1*x1;
double temp_545 = x1*temp_544;
double temp_556 = x2*x2;
double temp_557 = x2*temp_556;
double temp_562 = x3*x3;
double temp_566 = 4*temp_562;
double temp_578 = x3*temp_562;
double temp_589 = 36*x2*temp_562;
double temp_527 = x0*temp_517;
double temp_580 = 9*x2;
double temp_582 = 2*x3;
double temp_583 = temp_580+temp_582;
double temp_554 = x1*temp_545;
double temp_555 = 162*temp_554;
double temp_590 = 21*x2*x3;
double temp_591 = 36*temp_556;
double temp_592 = temp_590+temp_591+temp_566;
double temp_450 = 33*x2;
double temp_478 = 4*x3;
double temp_573 = 594*temp_544;
double temp_570 = 6*x2;
double temp_571 = temp_570+x3;
double temp_569 = 48*x2*x3;
double temp_574 = 108*temp_556;
double temp_575 = 7*temp_562;
double temp_587 = 72*x3*temp_556;
double temp_588 = 63*temp_557;
double temp_594 = 8*temp_578;
double temp_559 = x2*temp_557;
double temp_602 = x3*temp_578;
double temp_535 = 24*x2;
double temp_536 = 7*x3;
double temp_538 = temp_535+temp_536;
double temp_625 = 72*temp_557;
double temp_564 = 16*x2*x3;
double temp_565 = 21*temp_556;
double temp_567 = temp_564+temp_565+temp_566;
double temp_663 = 108*x3*temp_556;
double temp_563 = 594*temp_556*temp_562;
// I_y / (ρ π/2).
double Density2PithsMomentY =
(b*(4*(132*x1+temp_450+temp_478)*temp_517+286*temp_527+18*temp_538*temp_545
+temp_555+432*x3*temp_557+162*temp_559+temp_563+27*temp_544*temp_567+temp_511*
(temp_569+72*x1*temp_571+temp_573+temp_574+temp_575)+528*x2*temp_578+12*x1*(54*x3
*temp_556+36*temp_557+11*temp_578+temp_589)+2*x0*(216*temp_545+36*temp_544*temp_583
+temp_587+temp_588+temp_589+6*x1*temp_592+temp_594)+286*temp_602)+a*(22*(78*x1+12
*x2+x3)*temp_517+2002*temp_527+temp_555+2*temp_511*(9*x1*(temp_450+temp_478)+2*(
9*x2*x3+27*temp_556+temp_562)+temp_573)+36*temp_545*temp_583+9*temp_544*temp_592
+3*x1*(temp_587+temp_588+temp_589+temp_594)+2*(54*x3*temp_557+27*temp_559+54
*temp_556*temp_562+33*x2*temp_578+11*temp_602)+x0*(594*temp_545+63*x3*temp_556+24*x2
*temp_562+108*temp_544*temp_571+3*x1*(temp_569+temp_574+temp_575)+4*temp_578
+temp_625))+c*(22*temp_527+9*(21*x2+8*x3)*temp_545+54*temp_554+108*temp_544*(3*x2*x3
+3*temp_556+temp_562)+temp_517*(66*x1+4*temp_571)+temp_511*(108*temp_544+
temp_566+12*x1*temp_583+temp_590+temp_591)+2*(297*x3*temp_557+81*temp_559+temp_563+858
*x2*temp_578+1001*temp_602)+6*x1*(54*temp_557+99*x2*temp_562+44*temp_578+temp_663)
+x0*(9*temp_538*temp_544+108*temp_545+72*x2*temp_562+9*x1*temp_567+22*temp_578
+temp_625+temp_663)))/10010;

// Result.
Iy = Density2PithsMomentY * density * 0.5 * Math.Pi;
Ix = Iy / 2;
Iz = Ix;

EDIT declared the local variables; if this were C++ they should be const. Silly C# doesn't allow readonly local variables.

[eggrobin]

Since optimised code is illegible anyway (and the unoptimised version is hardly understandable in the first place) you may just want it to take as little room as possible. Here is an equivalent but more compact version (reduced the identifiers, folded lines accordingly):

// local variables to reduce the number of calculations.
double t511=x0*x0; double t517=x0*t511; double t544=x1*x1; double t545=x1*t544;
double t556=x2*x2; double t557=x2*t556; double t562=x3*x3; double t566=4*t562; 
double t578=x3*t562; double t589=36*x2*t562; double t527=x0*t517; 
double t580=9*x2; double t582=2*x3; double t583=t580+t582; double t554=x1*t545; 
double t555=162*t554; double t590=21*x2*x3; double t591=36*t556; double t592=t590+t591+t566; 
double t450=33*x2; double t478=4*x3; double t573=594*t544; double t570=6*x2; 
double t571=t570+x3; double t569=48*x2*x3; double t574=108*t556; double t575=7*t562;
double t587=72*x3*t556; double t588=63*t557; double t594=8*t578; 
double t559=x2*t557; double t602=x3*t578; double t535=24*x2; double t536=7*x3; 
double t538=t535+t536; double t625=72*t557; double t564=16*x2*x3; double t565=21*t556; 
double t567=t564+t565+t566; double t663=108*x3*t556; double t563=594*t556*t562;

// I_y / (ρ π/2).
double Density2PithsMomentY =
(b*(4*(132*x1+t450+t478)*t517+286*t527+18*t538*t545+t555+432*x3*t557+162*t559
+t563+27*t544*t567+t511*(t569+72*x1*t571+t573+t574+t575)+528*x2*t578+12*x1*(54*x3
*t556+36*t557+11*t578+t589)+2*x0*(216*t545+36*t544*t583+t587+t588+t589+6*x1*t592
+t594)+286*t602)+a*(22*(78*x1+12*x2+x3)*t517+2002*t527+t555+2*t511*(9*x1*(t450
+t478)+2*(9*x2*x3+27*t556+t562)+t573)+36*t545*t583+9*t544*t592+3*x1*(t587+t588+
t589+t594)+2*(54*x3*t557+27*t559+54*t556*t562+33*x2*t578+11*t602)+x0*(594*t545+63*
x3*t556+24*x2*t562+108*t544*t571+3*x1*(t569+t574+t575)+4*t578+t625))+c*(22*t527+9*
(21*x2+8*x3)*t545+54*t554+108*t544*(3*x2*x3+3*t556+t562)+t517*(66*x1+4*t571)
+t511*(108*t544+t566+12*x1*t583+t590+t591)+2*(297*x3*t557+81*t559+t563+858*x2*t578
+1001*t602)+6*x1*(54*t557+99*x2*t562+44*t578+t663)+x0*(9*t538*t544+108*t545+72*x2
*t562+9*x1*t567+22*t578+t625+t663)))/10010;

// Result.
Iy = Density2PithsMomentY * density * 0.5 * Math.Pi;
Ix = Iy / 2;
Iz = Ix;

EDIT: same as above.

[eggrobin]

Random thought: the assumption of uniform density is really bad for empty tanks. Perhaps the mass of the structure of the tank should be modelled as a thin shell? That makes the integration more interesting because we need to use conic frusta, but then I'll just dump it into Mathematica anyway, so it shouldn't be significantly harder.

[Swamp-Ig]

Yeh I had thought of that. What you realy want is a hollow shape for the tank, plus a filled shape inside that empties from top to bottom. The COM would also shift as it empties.

Then you could get even trickier and use the current acceleration vector. But then what if there's rotation and fuel sloshes up the sides?

And in zero G... assume uniform fuel density I guess...

:smile:

[Swamp-Ig]

More to say on this point but am using a phone at the moment so is a bit slow to type...

[eggrobin]

Something consistent with HoneyFox's EngineIgnitor fuel 'bubble' modelling would be nice actually. In zero G, fuel eventually settles on the tank walls due to surface tension. He considers rotation and acceleration, so the work has already been done. If sloshing had an actual effect for the mechanics of the rocket things would get far more interesting. There should be an option to add baffles to alleviate (but not suppress) that. Anyway, time to sleep.

[Swamp-Ig]

Non uniform density, and including the ullage bubble - now that gets tricky.

SO we have three components: The wall of the vessel W which is a thin walled shape with constant density ρ_w The vessel contents C model this as a solid shape that conforms to the vessel walls with density ρ_c . The ullage bubble E - a 'negative space' in the contents. density -ρ_c.

C is done already.

---

W - for simplification, lets consider the thickness to be negligible as it makes things easier. Finding a suitable dV_w is a bit trickier. Consider each dV - it will be a cylindrical tube, with height h, outer radius r, and thickness k. (k can be +ve or -ve, won't matter)

V_{cylindrical tube} = π (r²-(r-k)²) h

so doing the substitution of variables:

r=x(t) h=y(t+dt) - y(t) = y'(t) dt k=x(t+dt) - x(t) = x'(t) dt

dV = π (x(t)² - (x(t)²-2x(t)x'(t)dt+x'(t)²dt²) ) y'(t) dt

= π 2x(t)x'(t)y'(t) dt² - x'(t)²y'(t) dt³ = 0

Hmm... problem :smile:

It's pretty clear what the issue is - a zero width, zero height cylindrical tube is going to have no volume naturally.

So.. we need some thickness. We can't just use a constant thickness in x, because as y' approaches zero then the thickness vanishes vertically, so your ends don't count if you have a cylinder, and any slope inward or outward ends up counting less.

I had some thoughts about using a thickness vector perpendicular to f', but that would 'double count' when you went around corners.

Hmm.. a bit stuck :confounded:

---

E - that's reasonably easy to calculate, especially if we don't allow horizontal sloshing since things end up being surfaces of revolution again.

You'd assume the bubble is some kind of conic section in shape, and the volumes, centre of masses, and moments are all either well described or fairly straightforward to calculate.

[eggrobin]

Ullage simulation is complicated, so I put that into another issue, #43. Let's just go with an infinitely thin structural shell and uniformly distributed fuel in the meantime. The reason you're getting confused is that you need to integrate thin cone frusta, not thin cylinders (otherwise, as you say, the ends don't count and bad things happen generally).

[eggrobin]

An API question: do we know how to change the centre of mass in-flight? Sarbian strongly advises against using CoMOffset, since half the engine ignores it. I'm sure we can find some sort of workaround, but we should keep that in mind.

[Swamp-Ig]

Yeh, you can just update the Rigidbody directly. Set CoMOffset directly also since that's what RCS build aid looks for

[Swamp-Ig]

For W: moved from #43

Use surface area * some constant (which we'll incorporate into the density) as the approximation of volume

V_{conic frustra} = π (2r+k) sqrt(k²+h²)

so doing the substitution of variables:

r=x(t) h=y(t+dt) - y(t) = y'(t) dt k=x(t+dt) - x(t) = x'(t) dt

dV = π (2x(t) - x'(t)dt) sqrt(x'(t)²dx²+y'(t)²dt²)

(I had an inkling arc length was going to end up in there)

I'm not sure how to do the integral however.

[eggrobin]

The moment of inertia of a thin cone frustum. It is easy to compute the moments of a solid homogeneous cone frustum with mass m, radii a and b, and height h (using Wolfram|Alpha or integrating cylinders). One gets, if the symmetry is along the z axis (not the same convention as above),

iFrustum[a_, b_, h_, m_] := m {
  {(3 (a^4 + a^3 b + a^2 b^2 + a b^3 + b^4) + 2 (a^2 + 3 a b + 6 b^2) h^2)/(20 (a^2 + a b + b^2)), 0, 0},
  {0, (3 (a^4 + a^3 b + a^2 b^2 + a b^3 + b^4) + 2 (a^2 + 3 a b + 6 b^2) h^2)/(20 (a^2 + a b + b^2)), 0},
  {0, 0, (3 (a^4 + a^3 b + a^2 b^2 + a b^3 + b^4))/(10 (a^2 + a b + b^2))}
}

The volume of that frustum is

vFrustum[a_, b_, h_] := 1/3 π h (a^2 + a b + b^2)

One then removes from a frustum (a, b, h) a frustum (a-ε, b-ε, h), choosing the masses so that both frusta have the same density and the resulting shell has mass m. We thus have only thick sides remaining, the mass of those sides being m. We let ε tend to 0,

Limit[
  iFrustum[a, b, h, m / (1 - vFrustum[a - ε, b - ε, h] / vFrustum[a, b, h])]
    - iFrustum[a - ε, b - ε, h, m / (vFrustum[a, b, h] / vFrustum[a - ε, b - ε, h] - 1)],
  ε -> 0]

The result is

{
  {((3 (a + b) (a^2 + b^2) + 2 (a + 3 b) h^2) m)/(12 (a + b)), 0, 0},
  {0, ((3 (a + b) (a^2 + b^2) + 2 (a + 3 b) h^2) m)/(12 (a + b)), 0},
  {0, 0, 1/2 (a^2 + b^2) m}
}

We can check this result with a thin cylinder (a = b = r), a disk (a = r, b = 0, h = 0), or a thin hoop (a = b = r, h = 0).

Next step, integrating those frusta!

[Swamp-Ig]

just thinking - will we need to include the change of axis in the z-term? I don't think so but if you're already integrating using mathematica, you might as well put it in.

[eggrobin]

What change of axis?

[Swamp-Ig]

Parallel axis theorem

The x and y axies are always colinear, but the z axis changes as you do the integration.

Therefore the I_zz term should be 1/2 (a² + b² + 2 z²) m

I think it will just cancel out however, you'll end up with a term allowing that does the change of origin to centered about the CoM

[eggrobin]

Echoes from IRC So, let's go back to the earlier convention, y is the symmetry axis of the tank. The moments that are affected are Ix and Iz, and this cannot be neglected. This also applies to the previous calculation, and Ix and Iz will have to be redone there. Note that the perpendicular axis theorem does not apply here, since the tank is not planar.