From b39cd8a016a12ab6ecb08f07c2e298d12f8e2162 Mon Sep 17 00:00:00 2001 From: kruland Date: Fri, 17 Jan 2014 20:36:46 -0600 Subject: [PATCH] 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. --- .../aerodynamics/barrowman/FinSetCalc.java | 103 +++++++++------ core/src/net/sf/openrocket/util/MathUtil.java | 124 +++++++++--------- .../rocketcomponent/FinSetTest.java | 43 ++++++ 3 files changed, 166 insertions(+), 104 deletions(-) diff --git a/core/src/net/sf/openrocket/aerodynamics/barrowman/FinSetCalc.java b/core/src/net/sf/openrocket/aerodynamics/barrowman/FinSetCalc.java index 96abf7a62..293a04dd8 100644 --- a/core/src/net/sf/openrocket/aerodynamics/barrowman/FinSetCalc.java +++ b/core/src/net/sf/openrocket/aerodynamics/barrowman/FinSetCalc.java @@ -17,15 +17,21 @@ 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. */ protected static final int DIVISIONS = 48; - - + + protected double macLength = Double.NaN; // MAC length protected double macLead = Double.NaN; // MAC leading edge position protected double macSpan = Double.NaN; // MAC spanwise position @@ -42,7 +48,7 @@ public class FinSetCalc extends RocketComponentCalc { protected double[] chordTrail = new double[DIVISIONS]; protected double[] chordLength = new double[DIVISIONS]; - + protected final WarningSet geometryWarnings = new WarningSet(); private double[] poly = new double[6]; @@ -84,7 +90,7 @@ public class FinSetCalc extends RocketComponentCalc { public void calculateNonaxialForces(FlightConditions conditions, AerodynamicForces forces, WarningSet warnings) { - + if (span < 0.001) { forces.setCm(0); forces.setCN(0); @@ -98,28 +104,28 @@ public class FinSetCalc extends RocketComponentCalc { return; } - + // Add warnings (radius/2 == diameter/4) if (thickness > bodyRadius / 2) { warnings.add(Warning.THICK_FIN); } warnings.addAll(geometryWarnings); - - + + //////// Calculate CNa. ///////// // 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; double theta = conditions.getTheta(); double angle = baseRotation; - + // Compute basic CNa without interference effects if (finCount == 1 || finCount == 2) { // Basic CNa from geometry @@ -134,7 +140,8 @@ public class FinSetCalc extends RocketComponentCalc { cna = cna1 * finCount / 2.0; } - + // logger.debug("Component cna = {}", cna); + // Take into account fin-fin interference effects switch (interferenceFinCount) { case 1: @@ -215,7 +222,7 @@ public class FinSetCalc extends RocketComponentCalc { break; } */ - + // Body-fin interference effect double r = bodyRadius; double tau = r / (span + r); @@ -223,21 +230,26 @@ 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); + + - - // TODO: LOW: check for fin tip mach cone interference // (Barrowman thesis pdf-page 40) // TODO: LOW: fin-fin mach cone effect, MIL-HDBK page 5-25 - - - // Calculate CP position - double x = macLead + calculateCPPos(conditions) * macLength; - - + + // 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); + + + // Calculate roll forces, reduce forcing above stall angle // Without body-fin interference effect: @@ -246,9 +258,9 @@ public class FinSetCalc extends RocketComponentCalc { // With body-fin interference effect: forces.setCrollForce(finCount * (macSpan + r) * cna1 * (1 + tau) * cantAngle / conditions.getRefLength()); - - - + + + if (conditions.getAOA() > STALL_ANGLE) { // System.out.println("Fin stalling in roll"); forces.setCrollForce(forces.getCrollForce() * MathUtil.clamp( @@ -257,8 +269,8 @@ public class FinSetCalc extends RocketComponentCalc { forces.setCrollDamp(calculateDampingMoment(conditions)); forces.setCroll(forces.getCrollForce() - forces.getCrollDamp()); - - + + // System.out.printf(component.getName() + ": roll rate:%.3f force:%.3f damp:%.3f " + // "total:%.3f\n", // conditions.getRollRate(), forces.CrollForce, forces.CrollDamp, forces.Croll); @@ -285,7 +297,7 @@ public class FinSetCalc extends RocketComponentCalc { forces.setCside(0); forces.setCyaw(0); - + } @@ -304,7 +316,7 @@ public class FinSetCalc extends RocketComponentCalc { } - + /** * Pre-calculates the fin geometry values. */ @@ -329,7 +341,7 @@ public class FinSetCalc extends RocketComponentCalc { } } - + // Calculate the chord lead and trail positions and length Arrays.fill(chordLead, Double.POSITIVE_INFINITY); @@ -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)); @@ -388,7 +402,7 @@ public class FinSetCalc extends RocketComponentCalc { } } - + /* Calculate fin properties: * * macLength // MAC length @@ -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; @@ -516,8 +532,8 @@ public class FinSetCalc extends RocketComponentCalc { } - - + + private double calculateDampingMoment(FlightConditions conditions) { double rollRate = conditions.getRollRate(); @@ -527,7 +543,7 @@ public class FinSetCalc extends RocketComponentCalc { double mach = conditions.getMach(); double absRate = Math.abs(rollRate); - + /* * At low speeds and relatively large roll rates (i.e. near apogee) the * fin tips rotate well above stall angle. In this case sum the chords @@ -548,8 +564,8 @@ public class FinSetCalc extends RocketComponentCalc { return MathUtil.sign(rollRate) * finCount * sum; } - - + + if (mach <= CNA_SUBSONIC) { // System.out.println("BASIC: "+ // (component.getFinCount() * 2*Math.PI * rollRate * rollSum / @@ -594,8 +610,8 @@ public class FinSetCalc extends RocketComponentCalc { } - - + + /** * Return the relative position of the CP along the mean aerodynamic chord. * Below mach 0.5 it is at the quarter chord, above mach 2 calculated using an @@ -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; } @@ -689,7 +706,7 @@ public class FinSetCalc extends RocketComponentCalc { // // } - + @Override public double calculatePressureDragForce(FlightConditions conditions, double stagnationCD, double baseCD, WarningSet warnings) { @@ -727,7 +744,7 @@ public class FinSetCalc extends RocketComponentCalc { } // Airfoil assumed to have zero base drag - + // Scale to correct reference area drag *= finCount * span * thickness / conditions.getRefArea(); @@ -743,7 +760,7 @@ public class FinSetCalc extends RocketComponentCalc { double lead = component.toRelative(Coordinate.NUL, parent)[0].x; double trail = component.toRelative(new Coordinate(component.getLength()), - parent)[0].x; + parent)[0].x; /* * The counting fails if the fin root chord is very small, in that case assume @@ -768,8 +785,8 @@ public class FinSetCalc extends RocketComponentCalc { } if (interferenceFinCount < component.getFinCount()) { throw new BugException("Counted " + interferenceFinCount + " parallel fins, " + - "when component itself has " + component.getFinCount() + - ", fin points=" + Arrays.toString(component.getFinPoints())); + "when component itself has " + component.getFinCount() + + ", fin points=" + Arrays.toString(component.getFinPoints())); } } diff --git a/core/src/net/sf/openrocket/util/MathUtil.java b/core/src/net/sf/openrocket/util/MathUtil.java index 982d1b86f..434792a84 100644 --- a/core/src/net/sf/openrocket/util/MathUtil.java +++ b/core/src/net/sf/openrocket/util/MathUtil.java @@ -11,9 +11,9 @@ import org.slf4j.LoggerFactory; public class MathUtil { private static final Logger log = LoggerFactory.getLogger(MathUtil.class); - + public static final double EPSILON = 0.00000001; // 10mm^3 in m^3 - + /** * The square of x (x^2). On Sun's JRE using this method is as fast as typing x*x. * @param x x @@ -22,7 +22,7 @@ public class MathUtil { public static double pow2(double x) { return x * x; } - + /** * The cube of x (x^3). * @param x x @@ -31,11 +31,11 @@ public class MathUtil { public static double pow3(double x) { return x * x * x; } - + public static double pow4(double x) { return (x * x) * (x * x); } - + /** * Clamps the value x to the range min - max. * @param x Original value. @@ -50,7 +50,7 @@ public class MathUtil { return max; return x; } - + public static float clamp(float x, float min, float max) { if (x < min) return min; @@ -58,7 +58,7 @@ public class MathUtil { return max; return x; } - + public static int clamp(int x, int min, int max) { if (x < min) return min; @@ -66,8 +66,8 @@ public class MathUtil { return max; return x; } - - + + /** * Maps a value from one value range to another. * @@ -90,8 +90,8 @@ public class MathUtil { } return (value - fromMin) / (fromMax - fromMin) * (toMax - toMin) + toMin; } - - + + /** * Maps a coordinate from one value range to another. * @@ -115,8 +115,8 @@ public class MathUtil { double a = (value - fromMin) / (fromMax - fromMin); return toMax.multiply(a).add(toMin.multiply(1 - a)); } - - + + /** * Compute the minimum of two values. This is performed by direct comparison. * However, if one of the values is NaN and the other is not, the non-NaN value is @@ -127,7 +127,7 @@ public class MathUtil { return x; return (x < y) ? x : y; } - + /** * Compute the maximum of two values. This is performed by direct comparison. * However, if one of the values is NaN and the other is not, the non-NaN value is @@ -138,7 +138,7 @@ public class MathUtil { return y; return (x < y) ? y : x; } - + /** * Compute the minimum of three values. This is performed by direct comparison. * However, if one of the values is NaN and the other is not, the non-NaN value is @@ -151,9 +151,9 @@ public class MathUtil { return min(y, z); } } - - - + + + /** * Compute the minimum of three values. This is performed by direct comparison. * However, if one of the values is NaN and the other is not, the non-NaN value is @@ -162,8 +162,8 @@ public class MathUtil { public static double min(double w, double x, double y, double z) { return min(min(w, x), min(y, z)); } - - + + /** * Compute the maximum of three values. This is performed by direct comparison. * However, if one of the values is NaN and the other is not, the non-NaN value is @@ -176,7 +176,7 @@ public class MathUtil { return max(y, z); } } - + /** * Calculates the hypotenuse sqrt(x^2+y^2). This method is SIGNIFICANTLY * faster than Math.hypot(x,y). @@ -184,7 +184,7 @@ public class MathUtil { public static double hypot(double x, double y) { return Math.sqrt(x * x + y * y); } - + /** * Reduce the angle x to the range 0 - 2*PI. * @param x Original angle. @@ -194,7 +194,7 @@ public class MathUtil { double d = Math.floor(x / (2 * Math.PI)); return x - d * 2 * Math.PI; } - + /** * Reduce the angle x to the range -PI - PI. * @@ -207,8 +207,8 @@ public class MathUtil { double d = Math.rint(x / (2 * Math.PI)); return x - d * 2 * Math.PI; } - - + + /** * Return the square root of a value. If the value is negative, zero is returned. * This is safer in cases where rounding errors might make a value slightly negative. @@ -225,20 +225,22 @@ public class MathUtil { } return Math.sqrt(d); } - - - - public static boolean equals(double a, double b) { + + + public static boolean equals(double a, double b, double epsilon) { 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 * the special cases of zero and NaN. The value returned for those is arbitrary. @@ -251,20 +253,20 @@ public class MathUtil { public static double sign(double x) { return (x < 0) ? -1.0 : 1.0; } - + /* Math.abs() is about 3x as fast as this: public static double abs(double x) { return (x<0) ? -x : x; } */ - - + + public static double average(Collection values) { if (values.isEmpty()) { return Double.NaN; } - + double avg = 0.0; int count = 0; for (Number n : values) { @@ -273,12 +275,12 @@ public class MathUtil { } return avg / count; } - + public static double stddev(Collection values) { if (values.size() < 2) { return Double.NaN; } - + double avg = average(values); double stddev = 0.0; int count = 0; @@ -289,12 +291,12 @@ public class MathUtil { stddev = Math.sqrt(stddev / (count - 1)); return stddev; } - + public static double median(Collection values) { if (values.isEmpty()) { return Double.NaN; } - + List sorted = new ArrayList(values); Collections.sort(sorted, new Comparator() { @Override @@ -302,7 +304,7 @@ public class MathUtil { return Double.compare(o1.doubleValue(), o2.doubleValue()); } }); - + int n = sorted.size(); if (n % 2 == 0) { return (sorted.get(n / 2).doubleValue() + sorted.get(n / 2 - 1).doubleValue()) / 2; @@ -322,36 +324,36 @@ 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 domain, List range, double t ) { - - if ( domain == null || range == null || domain.size() != range.size() ) { + public static double interpolate(List domain, List range, double t) { + + 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) ) { - right ++; + while (t > domain.get(right)) { + right++; } - int left = right -1; - + int left = right - 1; + // Points are: double deltax = domain.get(right) - domain.get(left); double deltay = range.get(right) - range.get(left); - + // For numerical stability, if deltax is small, - if ( Math.abs(deltax) < EPSILON ) { - if ( deltay < -1.0 * EPSILON ) { + if (Math.abs(deltax) < 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,8 +362,8 @@ public class MathUtil { } } - return range.get(left) + ( t - domain.get(left) ) * deltay / deltax; - + return range.get(left) + (t - domain.get(left)) * deltay / deltax; + } - + } diff --git a/core/test/net/sf/openrocket/rocketcomponent/FinSetTest.java b/core/test/net/sf/openrocket/rocketcomponent/FinSetTest.java index 9febe3256..58dc0a440 100644 --- a/core/test/net/sf/openrocket/rocketcomponent/FinSetTest.java +++ b/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