Added a bunch of assert positives.

Found RegStep too big would make them fail
Added some adaptive regStep changing, expiermental

diff --git a/src/com/billkuker/rocketry/motorsim/Burn.java b/src/com/billkuker/rocketry/motorsim/Burn.java
index 03bf47c04cf0be1abb21b719b7a0d39b404c2e7d..49d9ef20f0071901f032586fedd0eb1be46c9e96 100644 (file)
--- a/src/com/billkuker/rocketry/motorsim/Burn.java
+++ b/src/com/billkuker/rocketry/motorsim/Burn.java
@@ -2,6 +2,7 @@ package com.billkuker.rocketry.motorsim;
 \r
 \r
 \r
+import java.util.Date;\r
 import java.util.SortedMap;\r
 import java.util.TreeMap;\r
 \r
@@ -13,6 +14,7 @@ import javax.measure.quantity.Length;
 import javax.measure.quantity.Mass;\r
 import javax.measure.quantity.MassFlowRate;\r
 import javax.measure.quantity.Pressure;\r
+import javax.measure.quantity.Quantity;\r
 import javax.measure.quantity.Temperature;\r
 import javax.measure.quantity.Velocity;\r
 import javax.measure.quantity.Volume;\r
@@ -64,8 +66,9 @@ public class Burn {
        \r
        private void burn(){\r
                log.info("Starting burn...");\r
+               long start = new Date().getTime();\r
                \r
-               Amount<Length> regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER);\r
+               Amount<Length> regStep = Amount.valueOf(0.1, SI.MILLIMETER);\r
 \r
                //if ( motor.getGrain() instanceof Grain.DiscreteRegression )\r
                        //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep();\r
@@ -80,27 +83,32 @@ public class Burn {
                \r
                data.put(Amount.valueOf(0, SI.SECOND), initial);\r
                \r
-               for ( int i = 0; i < 5000; i++ ){\r
-\r
+               step:\r
+               for ( int i = 0; i < 5000; i++ ) {\r
+                       regStep = regStep.times(1.01);\r
+                       \r
                        Interval prev = data.get(data.lastKey());\r
                        log.debug(prev);\r
                        log.debug("Step " + i + " ==============================");\r
                        Interval next = new Interval();\r
                        \r
                        Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);\r
+                       assert(positive(burnRate));\r
                        \r
                        log.debug("Burn Rate: " + burnRate);\r
                        \r
                        Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);\r
+                       assert(positive(dt));\r
                        next.dt = dt;\r
                        \r
-                       data.put(data.lastKey().plus(dt), next);\r
+\r
                        \r
                        log.debug("Dt: " + dt);\r
                        \r
                        next.regression = prev.regression.plus(regStep);\r
+                       assert(positive(next.regression));\r
                        \r
-                       log.info("Regression: " + next.regression);\r
+                       log.debug("Regression: " + next.regression);\r
                        \r
                        next.time = prev.time.plus(dt);\r
                        \r
@@ -110,9 +118,12 @@ public class Burn {
                        \r
                        //TODO Amount<Volume> volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
                        Amount<Volume> volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
+                       assert(positive(volumeBurnt));\r
                        //log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
                        \r
                        Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).divide(dt).to(MassFlowRate.UNIT);\r
+                       assert(positive(mGenRate));\r
+                       \r
                        //log.debug("Mass Gen Rate: " + mGenRate);\r
                        \r
                        //Calculate specific gas constant\r
@@ -136,6 +147,7 @@ public class Burn {
                                mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
                                //log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT));                \r
                        }\r
+                       assert(positive(mNozzle));\r
                        \r
                        Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);\r
                        \r
@@ -144,6 +156,12 @@ public class Burn {
                        next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));\r
                        \r
                        //Product can not go negative!\r
+                       if ( !positive(next.chamberProduct) ){\r
+                               log.warn("ChamberProduct Negative on step " + i + "!, Adjusting regstep down and repeating step!");\r
+                               regStep = regStep.divide(2);\r
+                               continue step;\r
+                       }\r
+                       assert(positive(next.chamberProduct));\r
                        if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) )\r
                                next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
                        \r
@@ -154,19 +172,34 @@ public class Burn {
                        log.debug("Product Density: " + combustionProductDensity);\r
                        \r
                        next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);\r
+                       assert(positive(next.chamberPressure));\r
                        \r
                        next.chamberPressure = Amount.valueOf(\r
                                        next.chamberPressure.doubleValue(SI.PASCAL),\r
                                        SI.PASCAL);\r
                        \r
+                       Amount<Pressure> dp = next.chamberPressure.minus(prev.chamberPressure);\r
+                       if ( dp.abs().isGreaterThan(Amount.valueOf(.5, SI.MEGA(SI.PASCAL)))){\r
+                               log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!");\r
+                               regStep = regStep.divide(2);\r
+                               continue step;\r
+                       }\r
+                       \r
                        next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());\r
+                       assert(positive(next.thrust));\r
                        \r
                        if ( i > 100 && next.chamberPressure.approximates(atmosphereicPressure)){\r
                                log.info("Pressure at Patm on step " + i);\r
                                break;\r
                        }\r
+                       \r
+                       data.put(data.lastKey().plus(dt), next);\r
+\r
+                       assert(positive(regStep));\r
                }\r
 \r
+               long time = new Date().getTime() - start;\r
+               log.info("Burn took " + time + " millis.");\r
        }\r
        \r
        @SuppressWarnings("unchecked")\r
@@ -192,4 +225,9 @@ public class Burn {
                return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);\r
        }\r
        \r
+       \r
+       private <Q extends Quantity> boolean positive(Amount<Q> a){\r
+               return ( a.isGreaterThan(a.minus(a)) || a.equals(a.minus(a)));\r
+       }\r
+\r
 }\r