Fix for Issue #175. Use a larger number to test equality of fin point y

positions to reduce numerical stability in computation of chord lengths.

core/src/net/sf/openrocket/aerodynamics/barrowman/FinSetCalc.java

@@ -17,8 +17,14 @@ import net.sf.openrocket.util.LinearInterpolator;
 import net.sf.openrocket.util.MathUtil;
 import net.sf.openrocket.util.PolyInterpolator;
 
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+
 
 public class FinSetCalc extends RocketComponentCalc {
+	
+	private final static Logger logger = LoggerFactory.getLogger(FinSetCalc.class);
+	
 	private static final double STALL_ANGLE = (20 * Math.PI / 180);
 	
 	/** Number of divisions in the fin chords. */
@@ -112,7 +118,7 @@ public class FinSetCalc extends RocketComponentCalc {
 		// One fin without interference (both sub- and supersonic):
 		double cna1 = calculateFinCNa1(conditions);
 		
-
+		//		logger.debug("Component cna1 = {}", cna1);
 		
 		// Multiple fins with fin-fin interference
 		double cna;
@@ -134,6 +140,7 @@ public class FinSetCalc extends RocketComponentCalc {
 			cna = cna1 * finCount / 2.0;
 		}
 		
+		//		logger.debug("Component cna = {}", cna);
 		
 		// Take into account fin-fin interference effects
 		switch (interferenceFinCount) {
@@ -223,6 +230,7 @@ public class FinSetCalc extends RocketComponentCalc {
 			tau = 0;
 		cna *= 1 + tau; // Classical Barrowman
 		//		cna *= pow2(1 + tau);	// Barrowman thesis (too optimistic??)
+		//		logger.debug("Component cna = {}", cna);
 		
 		
 		
@@ -234,7 +242,11 @@ public class FinSetCalc extends RocketComponentCalc {
 		
 		
 		// Calculate CP position
+		//		logger.debug("Component macLead = {}", macLead);
+		//		logger.debug("Component macLength = {}", macLength);
+		//FIXME - macLength is incorrect!
 		double x = macLead + calculateCPPos(conditions) * macLength;
+		//		logger.debug("Component x = {}", x);
 		
 		
 		
@@ -342,7 +354,9 @@ public class FinSetCalc extends RocketComponentCalc {
 			double x2 = points[point].x;
 			double y2 = points[point].y;
 			
-			if (MathUtil.equals(y1, y2))
+			// Don't use the default EPSILON since it is too small
+			// and causes too much numerical instability in the computation of x below
+			if (MathUtil.equals(y1, y2, 0.001))
 				continue;
 			
 			int i1 = (int) (y1 * 1.0001 / span * (DIVISIONS - 1));
@@ -412,6 +426,7 @@ public class FinSetCalc extends RocketComponentCalc {
 			double y = i * dy;
 			
 			macLength += length * length;
+			logger.debug("macLength = {}, length = {}, i = {}", macLength, length, i);
 			macSpan += y * length;
 			macLead += chordLead[i] * length;
 			area += length;
@@ -427,6 +442,7 @@ public class FinSetCalc extends RocketComponentCalc {
 		}
 		
 		macLength *= dy;
+		logger.debug("macLength = {}", macLength);
 		macSpan *= dy;
 		macLead *= dy;
 		area *= dy;
@@ -606,6 +622,7 @@ public class FinSetCalc extends RocketComponentCalc {
 	 */
 	private double calculateCPPos(FlightConditions cond) {
 		double m = cond.getMach();
+		//		logger.debug("m = {} ", m);
 		if (m <= 0.5) {
 			// At subsonic speeds CP at quarter chord
 			return 0.25;
@@ -624,7 +641,7 @@ public class FinSetCalc extends RocketComponentCalc {
 			val += poly[i] * x;
 			x *= m;
 		}
-		
+		//		logger.debug("val = {}", val);
 		return val;
 	}
 	

core/src/net/sf/openrocket/util/MathUtil.java

@@ -227,17 +227,19 @@ public class MathUtil {
 	}
 	
 	
-
+	public static boolean equals(double a, double b, double epsilon) {
-	public static boolean equals(double a, double b) {
 		double absb = Math.abs(b);
 		
-		if (absb < EPSILON / 2) {
+		if (absb < epsilon / 2) {
 			// Near zero
-			return Math.abs(a) < EPSILON / 2;
+			return Math.abs(a) < epsilon / 2;
 		}
-		return Math.abs(a - b) < EPSILON * absb;
+		return Math.abs(a - b) < epsilon * absb;
 	}
 	
+	public static boolean equals(double a, double b) {
+		return equals(a, b, EPSILON);
+	}
 	
 	/**
 	 * Return the sign of the number.  This corresponds to Math.signum, but ignores
@@ -322,23 +324,23 @@ public class MathUtil {
 	 * @param t domain value at which to interpolate
 	 * @return returns Double.NaN if either list is null or empty or different size, or if t is outsize the domain.
 	 */
-	public static double interpolate( List<Double> domain, List<Double> range, double t ) {
+	public static double interpolate(List<Double> domain, List<Double> range, double t) {
 		
-		if ( domain == null || range == null || domain.size() != range.size() ) {
+		if (domain == null || range == null || domain.size() != range.size()) {
 			return Double.NaN;
 		}
 		
 		int length = domain.size();
-		if ( length <= 1 || t < domain.get(0) || t > domain.get( length-1 ) ) {
+		if (length <= 1 || t < domain.get(0) || t > domain.get(length - 1)) {
 			return Double.NaN;
 		}
 		
 		// Look for the index of the right end point.
 		int right = 1;
-		while( t > domain.get(right) ) {
+		while (t > domain.get(right)) {
-			right ++;
+			right++;
 		}
-		int left = right -1;
+		int left = right - 1;
 		
 		// Points are:
 		
@@ -346,12 +348,12 @@ public class MathUtil {
 		double deltay = range.get(right) - range.get(left);
 		
 		// For numerical stability, if deltax is small,
-		if ( Math.abs(deltax) < EPSILON ) {
+		if (Math.abs(deltax) < EPSILON) {
-			if ( deltay < -1.0 * EPSILON ) {
+			if (deltay < -1.0 * EPSILON) {
 				// return neg infinity if deltay is negative
 				return Double.NEGATIVE_INFINITY;
 			}
-			else if ( deltay > EPSILON ) {
+			else if (deltay > EPSILON) {
 				// return infinity if deltay is large
 				return Double.POSITIVE_INFINITY;
 			} else {
@@ -360,7 +362,7 @@ public class MathUtil {
 			}
 		}
 		
-		return  range.get(left) + ( t - domain.get(left) ) * deltay / deltax;
+		return range.get(left) + (t - domain.get(left)) * deltay / deltax;
 		
 	}
 	

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

@@ -3,6 +3,10 @@ package net.sf.openrocket.rocketcomponent;
 import static org.junit.Assert.assertEquals;
 import static org.junit.Assert.assertFalse;
 import static org.junit.Assert.assertTrue;
+import net.sf.openrocket.aerodynamics.AerodynamicForces;
+import net.sf.openrocket.aerodynamics.FlightConditions;
+import net.sf.openrocket.aerodynamics.WarningSet;
+import net.sf.openrocket.aerodynamics.barrowman.FinSetCalc;
 import net.sf.openrocket.material.Material;
 import net.sf.openrocket.material.Material.Type;
 import net.sf.openrocket.rocketcomponent.ExternalComponent.Finish;
@@ -124,6 +128,45 @@ public class FinSetTest extends BaseTestCase {
 		
 	}
 	
+	@Test
+	public void testWildmanVindicatorShape() throws Exception {
+		// This fin shape is similar to the aft fins on the Wildman Vindicator.
+		// A user noticed that if the y values are similar but not equal,
+		// the compuation of CP was incorrect because of numerical instability.
+		//
+		//     +-----------------+
+		//      \                 \
+		//       \                 \
+		//        +                 \
+		//       /                   \
+		//      +---------------------+
+		//
+		FreeformFinSet fins = new FreeformFinSet();
+		fins.setFinCount(1);
+		Coordinate[] points = new Coordinate[] {
+				new Coordinate(0, 0),
+				new Coordinate(0.02143125, 0.01143),
+				new Coordinate(0.009524999999999999, 0.032543749999999996),
+				new Coordinate(0.041275, 0.032537399999999994),
+				new Coordinate(0.066675, 0)
+		};
+		fins.setPoints(points);
+		Coordinate coords = fins.getCG();
+		assertEquals(0.00130, fins.getFinArea(), 0.00001);
+		assertEquals(0.03423, coords.x, 0.00001);
+		assertEquals(0.01427, coords.y, 0.00001);
+		
+		BodyTube bt = new BodyTube();
+		bt.addChild(fins);
+		FinSetCalc calc = new FinSetCalc(fins);
+		FlightConditions conditions = new FlightConditions(null);
+		AerodynamicForces forces = new AerodynamicForces();
+		WarningSet warnings = new WarningSet();
+		calc.calculateNonaxialForces(conditions, forces, warnings);
+		System.out.println(forces);
+		assertEquals(0.023409, forces.getCP().x, 0.0001);
+	}
+	
 	@Test
 	public void testFreeFormCGWithNegativeY() throws Exception {
 		// This particular fin shape is currently not allowed in OR since the y values are negative