Consider shoulders in moment of inertia calculations

core/src/net/sf/openrocket/rocketcomponent/RocketComponent.java

@@ -2732,6 +2732,7 @@ public abstract class RocketComponent implements ChangeSource, Cloneable, Iterab
 	
 	protected static final double ringLongitudinalUnitInertia(double outerRadius,
 			double innerRadius, double length) {
+		// axis is through center of mass
 		// 1/12 * (3 * (r1^2 + r2^2) + h^2)
 		return (3 * (MathUtil.pow2(innerRadius) + MathUtil.pow2(outerRadius)) + MathUtil.pow2(length)) / 12;
 	}

core/src/net/sf/openrocket/rocketcomponent/SymmetricComponent.java

@@ -34,8 +34,8 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 	private double planCenter = Double.NaN;
 	protected double volume = Double.NaN;
 	private double fullVolume = Double.NaN;
-	private double longitudinalInertia = Double.NaN;
+	protected double longitudinalUnitInertia = Double.NaN;
-	private double rotationalInertia = Double.NaN;
+	protected double rotationalUnitInertia = Double.NaN;
 	protected Coordinate cg = null;
 
 	public SymmetricComponent() {
@@ -284,29 +284,29 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 	}
 
 	/*
-	 * Obtain longitudinal unit inertia Ixx
+	 * Obtain longitudinal unit inertia Iyy = Izz
 	 *
-	 * @return longitudinal unit inertia Ixx
+	 * @return longitudinal unit inertia Iyy
 	 */
 	@Override
 	public double getLongitudinalUnitInertia() {
-		if (Double.isNaN(longitudinalInertia)) {
+		if (Double.isNaN(longitudinalUnitInertia)) {
 			calculateProperties();
 		}
-		return longitudinalInertia;
+		return longitudinalUnitInertia;
 	}
 
 	/*
-	 * Obtain rotational unit inertia Iyy = Izz
+	 * Obtain rotational unit inertia Ixx
 	 *
 	 * @return rotational unit inertia
 	 */
 	@Override
 	public double getRotationalUnitInertia() {
-		if (Double.isNaN(rotationalInertia)) {
+		if (Double.isNaN(rotationalUnitInertia)) {
 			calculateProperties();
 		}
-		return rotationalInertia;
+		return rotationalUnitInertia;
 	}
 
 	/**
@@ -335,20 +335,22 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 
 	/**
 	 * helper method for calculateProperties()
-	 * returns the rotational moment of inertia of a solid frustum
+	 * returns the unit rotational moment of inertia of a solid frustum
 	 * handles special case of cylinder correctly
-	 * @param l length (height) of frustum
+	 * see https://web.mit.edu/8.13/8.13c/references-fall/aip/aip-handbook-section2c.pdf
+	 * page 2-41, table 2c-2
+	 * Caution!  Returns 10/3 times the correct answer.  Will need to be corrected at end
 	 * @param r1 radius of fore end of frustum
 	 * @param r2 radius of aft end of frustum
 	 * @return rotational moment of inertia
 	 */
-	private double calculateRotMOI(double l, double r1, double r2) {
+	private double calculateUnitRotMOI(double r1, double r2) {
 		// check for cylinder special case
 		if (Math.abs(r1 - r2) < MathUtil.EPSILON) {
-			return pow2(r1)/2.0;
+			return 10.0*pow2(r1)/6.0;
 		}
 		
-		return 3.0 * (Math.pow(r2, 5) - Math.pow(r1, 5)) / (10.0 * (Math.pow(r2, 3) - Math.pow(r1, 3)));
+		return (Math.pow(r2, 5) - Math.pow(r1, 5)) / (Math.pow(r2, 3) - Math.pow(r1, 3));
 	}
 
 	/**
@@ -423,8 +425,8 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 		planCenter = 0;
 		fullVolume = 0;
 		volume = 0;
-		longitudinalInertia = 0;
+		longitudinalUnitInertia = 0;
-		rotationalInertia = 0;
+		rotationalUnitInertia = 0;
 		
 		double cgx = 0;
 
@@ -480,22 +482,22 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 			final double dCGx = dV * (x1 + dCG);
 			
 			// rotational moment of inertia
-			final double Izzo = calculateRotMOI(l, r1o, r2o);
+			final double Ixxo = calculateUnitRotMOI(r1o, r2o);
-			final double Izzi = calculateRotMOI(l, r1i, r2i);
+			final double Ixxi = calculateUnitRotMOI(r1i, r2i);
-			final double Izz = Izzo * fullCG.weight - Izzi * innerCG.weight;
+			final double Ixx = Ixxo * fullCG.weight - Ixxi * innerCG.weight;
 
 			// longitudinal moment of inertia -- axis through CG of division
-			double Ixx = calculateLongMOI(l, r1o, r2o, fullCG) - calculateLongMOI(l, r1i, r2i, innerCG);
+			double Iyy = calculateLongMOI(l, r1o, r2o, fullCG) - calculateLongMOI(l, r1i, r2i, innerCG);
 
 			// move to axis through forward end of component
-			Ixx += dV * pow2(x1 + dCG);
+			Iyy += dV * pow2(x1 + dCG);
 			
 			// Add to the volume-related components
 			volume += dV;
 			fullVolume += dFullV;
 			cgx += dCGx;
-			rotationalInertia += Izz;
+			rotationalUnitInertia += Ixx;
-			longitudinalInertia += Ixx;
+			longitudinalUnitInertia += Iyy;
 
 			// Wetted area ( * PI at the end)
 			wetArea += (r1o + r2o) * Math.sqrt(pow2(r1o - r2o) + pow2(l));
@@ -506,19 +508,20 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 			final double planMoment = dA*x1 + 2.0*pow2(l)*(r1o/6.0 + r2o/3.0);
 			planCenter += planMoment;
 		}
-
-		// get unit moments of inertia
-		rotationalInertia /= volume;
-		longitudinalInertia /= volume;
 		
 		if (planArea > 0)
 			planCenter /= planArea;
 
+		// get unit moments of inertia
+		rotationalUnitInertia /= volume;
+		longitudinalUnitInertia /= volume;
+
 		// Correct for deferred constant factors
 		volume *= Math.PI / 3.0;
 		fullVolume *= Math.PI / 3.0;
 		cgx *= Math.PI / 3.0;
 		wetArea *= Math.PI;
+		rotationalUnitInertia *= 3.0 / 10.0;
 		
 		if (volume < 0.0000000001) { // 0.1 mm^3
 			volume = 0;
@@ -532,13 +535,13 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 		
 		// a component so small it has no volume can't contribute to moment of inertia
 		if (MathUtil.equals(volume, 0)) {
-			rotationalInertia = 0;
+			rotationalUnitInertia = 0;
-			longitudinalInertia = 0;
+			longitudinalUnitInertia = 0;
 			return;
 		}
 		
 		// Shift longitudinal inertia to CG
-		longitudinalInertia = longitudinalInertia - pow2(cg.x);
+		longitudinalUnitInertia = longitudinalUnitInertia - pow2(cg.x);
 		
 	}
 
@@ -554,8 +557,8 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
 			planCenter = Double.NaN;
 			volume = Double.NaN;
 			fullVolume = Double.NaN;
-			longitudinalInertia = Double.NaN;
+			longitudinalUnitInertia = Double.NaN;
-			rotationalInertia = Double.NaN;
+			rotationalUnitInertia = Double.NaN;
 			cg = null;
 		}
 	}

core/src/net/sf/openrocket/rocketcomponent/Transition.java

@@ -13,10 +13,10 @@ import net.sf.openrocket.startup.Application;
 import net.sf.openrocket.util.Coordinate;
 import net.sf.openrocket.util.MathUtil;
 
-
 public class Transition extends SymmetricComponent implements InsideColorComponent {
 	private static final Translator trans = Application.getTranslator();
 	private static final double CLIP_PRECISION = 0.0001;
+	final static double MINFEATURE = 0.001;
 
 
 	private Shape type;
@@ -663,9 +663,9 @@ public class Transition extends SymmetricComponent implements InsideColorCompone
 	@Override
 	public Collection<Coordinate> getComponentBounds() {
 		Collection<Coordinate> bounds = super.getComponentBounds();
-		if (foreShoulderLength > 0.001)
+		if (foreShoulderLength > MINFEATURE)
 			addBound(bounds, -foreShoulderLength, foreShoulderRadius);
-		if (aftShoulderLength > 0.001)
+		if (aftShoulderLength > MINFEATURE)
 			addBound(bounds, getLength() + aftShoulderLength, aftShoulderRadius);
 		return bounds;
 	}
@@ -673,54 +673,119 @@ public class Transition extends SymmetricComponent implements InsideColorCompone
 	@Override
 	protected void calculateProperties() {
 		super.calculateProperties();
-		if (getForeShoulderLength() > 0.001) {
-			final double or = getForeShoulderRadius();
-			final double ir = Math.max(getForeShoulderRadius() - getForeShoulderThickness(), 0);
-			
-			cg = cg.average(ringCG(getForeShoulderRadius(), ir, -getForeShoulderLength(), 0,
-					getMaterial().getDensity()));
 
-			volume += ringVolume( or, ir, getForeShoulderLength() );
+		// only adjust properties if there is in fact at least one shoulder
+		if ((getForeShoulderLength() > MINFEATURE) || (getAftShoulderLength() > MINFEATURE)) {
+			// we'll work with volumes and not masses because density is uniform and it'll save some
+			// multiplications and divisions
+
+			// accumulate data for later calculation of properties with shoulders added
+			double transVolume = volume;
+			double transLongMOI = longitudinalUnitInertia * transVolume;
+			double transRotMOI = rotationalUnitInertia * transVolume;
+			final Coordinate transCG = cg;
+
+			double foreCapVolume = 0.0;
+			Coordinate foreCapCG = Coordinate.ZERO;
+			double foreCapLongMOI = 0.0;
+			double foreCapRotMOI = 0.0;
+			if (isForeShoulderCapped()) {
+				final double ir = Math.max(getForeShoulderRadius() - getForeShoulderThickness(), 0);
+
+				foreCapCG = ringCG(ir, 0, -getForeShoulderLength(),
+								   getForeShoulderThickness() - getForeShoulderLength(),
+								   getMaterial().getDensity());
+				
+				foreCapVolume = ringVolume(ir, 0, getForeShoulderThickness() );
+
+				foreCapLongMOI = ringLongitudinalUnitInertia(ir, 0, getForeShoulderThickness())*foreCapVolume;
+				
+				foreCapRotMOI += ringRotationalUnitInertia(ir, 0.0) * foreCapVolume;
+			}
 			
+			double foreShoulderVolume = 0.0;
+			Coordinate foreShoulderCG = Coordinate.ZERO;
+			double foreShoulderLongMOI = 0.0;
+			double foreShoulderRotMOI = 0.0;
+			if (getForeShoulderLength() > MINFEATURE) {
+				final double or = getForeShoulderRadius();
+				final double ir = Math.max(getForeShoulderRadius() - getForeShoulderThickness(), 0);
+
+				foreShoulderCG = ringCG(getForeShoulderRadius(), ir, -getForeShoulderLength(), 0,
+										getMaterial().getDensity());
+				
+				foreShoulderVolume = ringVolume( or, ir, getForeShoulderLength() );
+
+				foreShoulderLongMOI = ringLongitudinalUnitInertia(or, ir, getForeShoulderLength())*foreShoulderVolume;
+				
+				foreShoulderRotMOI = ringRotationalUnitInertia(or, ir) * foreShoulderVolume;
+			}
+			
+
+			double aftShoulderVolume = 0.0;
+			Coordinate aftShoulderCG = Coordinate.ZERO;
+			double aftShoulderLongMOI = 0.0;
+			double aftShoulderRotMOI = 0.0;
+			if (getAftShoulderLength() > MINFEATURE) {
+				final double or = getAftShoulderRadius();
+				final double ir = Math.max(getAftShoulderRadius() - getAftShoulderThickness(), 0);
+
+				aftShoulderCG = ringCG(getAftShoulderRadius(), ir, getLength(),
+									   getLength() + getAftShoulderLength(),
+									   getMaterial().getDensity());
+
+				aftShoulderVolume = ringVolume(or, ir, getAftShoulderLength());
+
+				aftShoulderLongMOI = ringLongitudinalUnitInertia(or, ir, getAftShoulderLength())*aftShoulderVolume;
+				
+				aftShoulderRotMOI = ringRotationalUnitInertia(or, ir) * aftShoulderVolume;
+			}
+
+			double aftCapVolume = 0.0;
+			Coordinate aftCapCG = Coordinate.ZERO;
+			double aftCapLongMOI = 0.0;
+			double aftCapRotMOI = 0.0;
+			if (isAftShoulderCapped()) {
+				final double ir = Math.max(getAftShoulderRadius() - getAftShoulderThickness(), 0);
+
+				aftCapCG = ringCG(ir, 0,
+								  getLength() + getAftShoulderLength() - getAftShoulderThickness(),
+								  getLength() + getAftShoulderLength(), getMaterial().getDensity());
+
+				aftCapVolume = ringVolume(ir, 0, getAftShoulderThickness() );
+
+				aftCapLongMOI = ringLongitudinalUnitInertia(ir, 0, getForeShoulderThickness())*aftCapVolume;
+				
+				aftCapRotMOI = ringRotationalUnitInertia(ir, 0.0) * aftCapVolume;			
+			}
+
+			// Combine results
+			volume = foreCapVolume + foreShoulderVolume + transVolume + aftShoulderVolume + aftCapVolume;
+
+			final double cgx = foreCapCG.x * foreCapCG.weight +
+				foreShoulderCG.x * foreShoulderCG.weight +
+				transCG.x * transCG.weight +
+				aftShoulderCG.x * aftShoulderCG.weight +
+				aftCapCG.x * aftCapCG.weight;
+			
+			final double mass = foreCapCG.weight + foreShoulderCG.weight + transCG.weight + aftShoulderCG.weight + aftCapCG.weight;
+
+			cg = new Coordinate(cgx / mass, 0, 0, mass);
+
+			// need to use parallel axis theorem to move longitudinal MOI to CG of component
+			foreCapLongMOI += pow2(cg.x - foreCapCG.x) * foreCapVolume;
+			foreShoulderLongMOI += pow2(cg.x - foreShoulderCG.x) * foreShoulderVolume;
+			transLongMOI += pow2(cg.x - transCG.x) * transVolume;
+			aftShoulderLongMOI += pow2(cg.x - aftShoulderCG.x) * aftShoulderVolume;
+			aftCapLongMOI += pow2(cg.x - aftCapCG.x) * aftCapVolume;
+
+			final double longMOI = foreCapLongMOI + foreShoulderLongMOI + transLongMOI + aftShoulderLongMOI + aftCapLongMOI;
+			longitudinalUnitInertia = longMOI/volume;
+
+			final double rotMOI = foreCapRotMOI + foreShoulderRotMOI + transRotMOI + aftShoulderRotMOI + aftCapRotMOI;
+			rotationalUnitInertia = rotMOI/volume;
 		}
-		if (isForeShoulderCapped()) {
+ 	}
-			final double ir = Math.max(getForeShoulderRadius() - getForeShoulderThickness(), 0);
-			cg = cg.average(ringCG(ir, 0, -getForeShoulderLength(),
-					getForeShoulderThickness() - getForeShoulderLength(),
-					getMaterial().getDensity()));
-
-			volume += ringVolume(ir, 0, getForeShoulderThickness() );
-		}
-
-		if (getAftShoulderLength() > 0.001) {
-			final double or = getAftShoulderRadius();
-			final double ir = Math.max(getAftShoulderRadius() - getAftShoulderThickness(), 0);
-			
-			cg = cg.average(ringCG(getAftShoulderRadius(), ir, getLength(),
-					getLength() + getAftShoulderLength(), getMaterial().getDensity()));
-			
-			volume += ringVolume(or, ir, getAftShoulderLength() );			
-		}
-		if (isAftShoulderCapped()) {
-			final double or = getAftShoulderRadius();
-			final double ir = Math.max(getAftShoulderRadius() - getAftShoulderThickness(), 0);
-			
-			cg = cg.average(ringCG(ir, 0,
-					getLength() + getAftShoulderLength() - getAftShoulderThickness(),
-					getLength() + getAftShoulderLength(), getMaterial().getDensity()));
-
-			volume += ringVolume(ir, 0, getAftShoulderThickness() );			
-		}
-	}
-
-
-	/*
-	 * The moments of inertia are not explicitly corrected for the shoulders.
-	 * However, since the mass is corrected, the inertia is automatically corrected
-	 * to very nearly the correct value.
-	 */
-
-
 
 	/**
 	 * Returns the name of the component ("Transition").
@@ -884,7 +949,7 @@ public class Transition extends SymmetricComponent implements InsideColorCompone
 					length = radius;
 				}
 
-				if (param < 0.001)
+				if (param < MINFEATURE)
 					return CONICAL.getRadius(x, radius, length, param);
 
 				// Radius of circle is:

core/test/net/sf/openrocket/rocketcomponent/SymmetricComponentVolumeTest.java

@@ -337,6 +337,7 @@ public class SymmetricComponentVolumeTest extends BaseTestCase {
 
 	@Test
 	public void testTransitionVsTubeFilled() {
+		System.err.println(">>>> tube filled");
 		// BodyTubes use closed form solutions for mass properties, while Transitions use
 		// numerical integration from SymmetricComponent. Properties should agree.
 		final double radius = 1.0;
@@ -354,18 +355,21 @@ public class SymmetricComponentVolumeTest extends BaseTestCase {
 		assertEquals("Length is incorrect", bt1.getLength(), trans1.getLength(), EPSILON);
 		assertEquals("Forward radius is incorrect", bt1.getRadius(0), trans1.getRadius(0), EPSILON);
 		assertEquals("Aft radius is incorrect", bt1.getRadius(bt1.getLength()), trans1.getRadius(trans1.getLength()), EPSILON);
+		
 		assertEquals("Volume is incorrect", bt1.getComponentVolume(), trans1.getComponentVolume(), EPSILON);
 		assertEquals("CG is incorrect", bt1.getComponentCG().x, trans1.getComponentCG().x, EPSILON);
+
 		assertEquals("Longitudinal moment of inertia is incorrect", bt1.getLongitudinalUnitInertia(), trans1.getLongitudinalUnitInertia(), EPSILON);
-		assertEquals("Rotational moment of inertia is incorrect", bt1.getRotationalUnitInertia(), trans1.getRotationalUnitInertia(), EPSILON);
+		assertEquals("Rotational unit moment of inertia is incorrect", bt1.getRotationalUnitInertia(), trans1.getRotationalUnitInertia(), EPSILON);
+
 		assertEquals("Wetted area is incorrect", bt1.getComponentWetArea(), trans1.getComponentWetArea(), EPSILON);
 		assertEquals("Planform area is incorrect", bt1.getComponentPlanformArea(), trans1.getComponentPlanformArea(), EPSILON);
 		assertEquals("Planform centroid is incorrect", bt1.getComponentPlanformCenter(), trans1.getComponentPlanformCenter(), EPSILON);
-
 	}
 	
 	@Test
 	public void testTransitionVsTubeHollow() {
+		System.err.println(">>>> tube hollow");
 		final double radius = 1.0;
 		final double innerRadius = 0.1;
 		final double length = 10.0;
@@ -386,6 +390,7 @@ public class SymmetricComponentVolumeTest extends BaseTestCase {
 		assertEquals("Longitudinal unit moment of inertia is incorrect", bt2.getLongitudinalUnitInertia(), trans2.getLongitudinalUnitInertia(), EPSILON);
 		assertEquals("Rotational unit moment of inertia is incorrect", bt2.getRotationalUnitInertia(), trans2.getRotationalUnitInertia(), EPSILON);
 		assertEquals("Wetted area is incorrect", bt2.getComponentWetArea(), trans2.getComponentWetArea(), EPSILON);
+
 		assertEquals("Planform area is incorrect", bt2.getComponentPlanformArea(), trans2.getComponentPlanformArea(), EPSILON);
 		assertEquals("Planform centroid is incorrect", bt2.getComponentPlanformCenter(), trans2.getComponentPlanformCenter(), EPSILON);
 	}