From ac2339b8b44ae0325a689b635a510144a519e86c Mon Sep 17 00:00:00 2001 From: Joe Pfeiffer Date: Tue, 21 Aug 2018 17:29:07 -0600 Subject: [PATCH 1/4] Fix Average Thrust Calculation (fixes issue #441) Remove test for short time interval before first data point in thrust curve. Comment said it was for numerical stability; multiplying by a small number and then adding doesn't introduce any instabilities I'm aware of in this code. Add parentheses to clarify that values are being multiplied by time intervals, not divided. --- .../net/sf/openrocket/motor/ThrustCurveMotor.java | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java index ac4b37a6b..3cf858a9c 100644 --- a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java +++ b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java @@ -321,23 +321,20 @@ public class ThrustCurveMotor implements Motor, Comparable, Se // portion from startTime through time[timeIndex] double avgImpulse = 0.0; - // For numeric stability. - if( time[timeIndex+1] - startTime > 0.001 ) { - avgImpulse = (MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]) + thrust[timeIndex+1]) - / 2.0 * (time[timeIndex+1] - startTime); - } + avgImpulse = ((MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]) + thrust[timeIndex+1]) + / 2.0) * (time[timeIndex+1] - startTime); // Now add the whole steps; timeIndex++; - while( timeIndex < time.length -1 && endTime >= time[timeIndex+1] ) { - avgImpulse += (thrust[timeIndex] + thrust[timeIndex+1]) / 2.0 * (time[timeIndex+1]-time[timeIndex]); + while ( timeIndex < time.length -1 && endTime >= time[timeIndex+1] ) { + avgImpulse += ((thrust[timeIndex] + thrust[timeIndex+1]) / 2.0) * (time[timeIndex+1]-time[timeIndex]); timeIndex++; } // Now add the bit after the last time index if ( timeIndex < time.length -1 ) { double endInstImpulse = MathUtil.map( endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - avgImpulse += (thrust[timeIndex] + endInstImpulse) / 2.0 * (endTime - time[timeIndex]); + avgImpulse += ((thrust[timeIndex] + endInstImpulse) / 2.0) * (endTime - time[timeIndex]); } return avgImpulse / (endTime - startTime); From c089e3ecb777b47fe7f2da839eafd234a6ad9e02 Mon Sep 17 00:00:00 2001 From: Joe Pfeiffer Date: Sun, 26 Aug 2018 14:53:00 -0600 Subject: [PATCH 2/4] Code clarification; should make no difference to results --- .../sf/openrocket/motor/ThrustCurveMotor.java | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java index 3cf858a9c..8bb8afddd 100644 --- a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java +++ b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java @@ -313,28 +313,28 @@ public class ThrustCurveMotor implements Motor, Comparable, Se if ( endTime <= time[timeIndex+1] ) { // we are completely within this time slice so the computation of the average is pretty easy: - double avgImpulse = MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - avgImpulse += MathUtil.map(endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - avgImpulse /= 2.0; - return avgImpulse; + double startThrust = MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); + double endThrust = MathUtil.map(endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); + return (startThrust + endThrust) / 2.0; } - - // portion from startTime through time[timeIndex] + double avgImpulse = 0.0; - avgImpulse = ((MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]) + thrust[timeIndex+1]) - / 2.0) * (time[timeIndex+1] - startTime); + + // portion from startTime through time[timeIndex+1] + double startThrust = MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); + avgImpulse = (time[timeIndex+1] - startTime) * (startThrust + thrust[timeIndex+1]) / 2.0; // Now add the whole steps; timeIndex++; while ( timeIndex < time.length -1 && endTime >= time[timeIndex+1] ) { - avgImpulse += ((thrust[timeIndex] + thrust[timeIndex+1]) / 2.0) * (time[timeIndex+1]-time[timeIndex]); + avgImpulse += (time[timeIndex+1] - time[timeIndex]) * (thrust[timeIndex] + thrust[timeIndex+1]) / 2.0; timeIndex++; } // Now add the bit after the last time index if ( timeIndex < time.length -1 ) { - double endInstImpulse = MathUtil.map( endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - avgImpulse += ((thrust[timeIndex] + endInstImpulse) / 2.0) * (endTime - time[timeIndex]); + double endThrust = MathUtil.map( endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); + avgImpulse += ((thrust[timeIndex] + endThrust) / 2.0) * (endTime - time[timeIndex]); } return avgImpulse / (endTime - startTime); From d2cdea21131c85d2550e9e29dc71f0411d224bdd Mon Sep 17 00:00:00 2001 From: Joe Pfeiffer Date: Sun, 26 Aug 2018 16:32:44 -0600 Subject: [PATCH 3/4] Little bit more massaging for clarity (replace avgImpulse with impulse) --- core/src/net/sf/openrocket/motor/ThrustCurveMotor.java | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java index 8bb8afddd..ec5f82605 100644 --- a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java +++ b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java @@ -318,26 +318,26 @@ public class ThrustCurveMotor implements Motor, Comparable, Se return (startThrust + endThrust) / 2.0; } - double avgImpulse = 0.0; + double impulse = 0.0; // portion from startTime through time[timeIndex+1] double startThrust = MathUtil.map(startTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - avgImpulse = (time[timeIndex+1] - startTime) * (startThrust + thrust[timeIndex+1]) / 2.0; + impulse = (time[timeIndex+1] - startTime) * (startThrust + thrust[timeIndex+1]) / 2.0; // Now add the whole steps; timeIndex++; while ( timeIndex < time.length -1 && endTime >= time[timeIndex+1] ) { - avgImpulse += (time[timeIndex+1] - time[timeIndex]) * (thrust[timeIndex] + thrust[timeIndex+1]) / 2.0; + impulse += (time[timeIndex+1] - time[timeIndex]) * (thrust[timeIndex] + thrust[timeIndex+1]) / 2.0; timeIndex++; } // Now add the bit after the last time index if ( timeIndex < time.length -1 ) { double endThrust = MathUtil.map( endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - avgImpulse += ((thrust[timeIndex] + endThrust) / 2.0) * (endTime - time[timeIndex]); + impulse += ((thrust[timeIndex] + endThrust) / 2.0) * (endTime - time[timeIndex]); } - return avgImpulse / (endTime - startTime); + return impulse / (endTime - startTime); } @Override From 7a04bd567c9a37ad825b52bbe21da47a8faf667e Mon Sep 17 00:00:00 2001 From: Joe Pfeiffer Date: Mon, 27 Aug 2018 09:16:16 -0600 Subject: [PATCH 4/4] missed reversing the operands in the calculation of last bit of impulse --- core/src/net/sf/openrocket/motor/ThrustCurveMotor.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java index ec5f82605..3e851f260 100644 --- a/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java +++ b/core/src/net/sf/openrocket/motor/ThrustCurveMotor.java @@ -334,7 +334,7 @@ public class ThrustCurveMotor implements Motor, Comparable, Se // Now add the bit after the last time index if ( timeIndex < time.length -1 ) { double endThrust = MathUtil.map( endTime, time[timeIndex], time[timeIndex+1], thrust[timeIndex], thrust[timeIndex+1]); - impulse += ((thrust[timeIndex] + endThrust) / 2.0) * (endTime - time[timeIndex]); + impulse += (endTime - time[timeIndex]) * (thrust[timeIndex] + endThrust) / 2.0; } return impulse / (endTime - startTime);