Make atmospheric pressure editable
[sw/motorsim] / src / com / billkuker / rocketry / motorsim / Burn.java
index 96fd0ed60011b877b3d98dfc2dd626cd906a2572..0416d1ec0892a087710627ceb61555854534f99c 100644 (file)
@@ -2,6 +2,9 @@ package com.billkuker.rocketry.motorsim;
 \r
 \r
 \r
+import java.util.Date;\r
+import java.util.HashSet;\r
+import java.util.Set;\r
 import java.util.SortedMap;\r
 import java.util.TreeMap;\r
 \r
@@ -13,6 +16,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
@@ -23,18 +27,85 @@ import org.apache.log4j.Logger;
 import org.jscience.physics.amount.Amount;\r
 import org.jscience.physics.amount.Constants;\r
 \r
+import com.billkuker.rocketry.motorsim.Validating.ValidationException;\r
+\r
 public class Burn {\r
-       \r
        private static Logger log = Logger.getLogger(Burn.class);\r
-       protected final Motor motor;\r
        \r
-       private static final Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
+       private static final BurnSettings settings = new BurnSettings();\r
+       public static final BurnSettings getBurnSettings(){\r
+               return settings;\r
+       }\r
        \r
+       /**\r
+        * A class representing all the settigns one can change on a burn\r
+        * @author bkuker\r
+        */\r
+       public static class BurnSettings {\r
+               private BurnSettings(){};\r
+\r
+               public enum BurnVolumeMethod {\r
+                       DeltaVolume,\r
+                       SurfaceTimesRegression;\r
+               }\r
+               \r
+               private BurnVolumeMethod volumeMethod = BurnVolumeMethod.SurfaceTimesRegression;\r
+               private double regStepIncreaseFactor = 1.01;\r
+               private double regStepDecreaseFactor = .5;\r
+               private Amount<Pressure> chamberPressureMaxDelta = Amount.valueOf(.5, SI.MEGA(SI.PASCAL));\r
+               private Amount<Pressure> endPressure = Amount.valueOf(.1, RocketScience.PSI);\r
+               private Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
+               \r
+               public void setVolumeMethod(BurnVolumeMethod volumeMethod) {\r
+                       this.volumeMethod = volumeMethod;\r
+               }\r
+               public BurnVolumeMethod getVolumeMethod() {\r
+                       return volumeMethod;\r
+               }\r
+               public double getRegStepIncreaseFactor() {\r
+                       return regStepIncreaseFactor;\r
+               }\r
+               public void setRegStepIncreaseFactor(double regStepIncreaseFactor) {\r
+                       this.regStepIncreaseFactor = regStepIncreaseFactor;\r
+               }\r
+               public double getRegStepDecreaseFactor() {\r
+                       return regStepDecreaseFactor;\r
+               }\r
+               public void setRegStepDecreaseFactor(double regStepDecreaseFactor) {\r
+                       this.regStepDecreaseFactor = regStepDecreaseFactor;\r
+               }\r
+               public Amount<Pressure> getChamberPressureMaxDelta() {\r
+                       return chamberPressureMaxDelta;\r
+               }\r
+               public void setChamberPressureMaxDelta(Amount<Pressure> chamberPressureMaxDelta) {\r
+                       this.chamberPressureMaxDelta = chamberPressureMaxDelta;\r
+               }\r
+               public Amount<Pressure> getEndPressure() {\r
+                       return endPressure;\r
+               }\r
+               public void setEndPressure(Amount<Pressure> endPressure) {\r
+                       this.endPressure = endPressure;\r
+               }\r
+               public Amount<Pressure> getAtmosphereicPressure() {\r
+                       return atmosphereicPressure;\r
+               }\r
+               public void setAtmosphereicPressure(Amount<Pressure> atmosphereicPressure) {\r
+                       this.atmosphereicPressure = atmosphereicPressure;\r
+               }\r
+       }\r
        \r
-       private static double combustionEfficency = 0.97;\r
+       protected final Motor motor;\r
        \r
-       private static double densityRatio = 0.96;\r
+       private boolean burning = false;\r
+       private boolean done = false;\r
        \r
+       public interface BurnProgressListener{\r
+               public void setProgress(float p);\r
+               public void burnComplete();\r
+       }\r
+               \r
+       private Set<BurnProgressListener> bpls = new HashSet<Burn.BurnProgressListener>();\r
+\r
        public class Interval{\r
                public Amount<Duration> time;\r
                public Amount<Duration> dt;\r
@@ -51,6 +122,8 @@ public class Burn {
        protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();\r
        \r
        public SortedMap<Amount<Duration>,Interval> getData(){\r
+               if ( !done )\r
+                       throw new IllegalStateException("Burn not complete!");\r
                return data;\r
        }\r
        \r
@@ -59,69 +132,107 @@ public class Burn {
        }\r
 \r
        public Amount<Duration> burnTime(){\r
-               return data.lastKey();\r
+               return getData().lastKey();\r
        }\r
        \r
        public Burn(Motor m){\r
+               try {\r
+                       m.validate();\r
+               } catch (ValidationException e) {\r
+                       throw new IllegalArgumentException("Invalid Motor: " + e.getMessage());\r
+               }\r
                motor = m;\r
-               burn();\r
        }\r
        \r
-       private void burn(){\r
+       public void addBurnProgressListener( BurnProgressListener bpl ){\r
+               bpls.add(bpl);\r
+               if ( done )\r
+                       bpl.burnComplete();\r
+       }\r
+       \r
+       public void burn(){\r
+               synchronized(this){\r
+                       if ( burning )\r
+                               throw new IllegalStateException("Already burning!");\r
+                       burning = true;\r
+               }\r
                log.info("Starting burn...");\r
+               for (BurnProgressListener bpl : bpls ){\r
+                       bpl.setProgress(0);\r
+               }\r
                \r
-               Amount<Length> regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER);\r
-\r
-               //if ( motor.getGrain() instanceof Grain.DiscreteRegression )\r
-                       //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep();\r
+               int endPressureSteps = 0;\r
+               long start = new Date().getTime();\r
                \r
+               Amount<Length> regStep = Amount.valueOf(0.01, SI.MILLIMETER);\r
+\r
                Interval initial = new Interval();\r
                initial.time = Amount.valueOf(0, SI.SECOND);\r
                initial.dt = Amount.valueOf(0, SI.SECOND);\r
                initial.regression = Amount.valueOf(0, SI.MILLIMETER);\r
-               initial.chamberPressure = atmosphereicPressure;\r
+               initial.chamberPressure = settings.getAtmosphereicPressure();\r
                initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
                initial.thrust = Amount.valueOf(0, SI.NEWTON);\r
                \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
+                       assert(positive(regStep));\r
+                       regStep = regStep.times(settings.getRegStepIncreaseFactor());\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
+                       //Update BurnProgressListeners\r
+                       Amount<Dimensionless> a = next.regression.divide(motor.getGrain().webThickness()).to(Dimensionless.UNIT);\r
+                       for (BurnProgressListener bpl : bpls ){\r
+                               bpl.setProgress((float)a.doubleValue(Dimensionless.UNIT));\r
+                       }\r
+\r
                        \r
                        next.time = prev.time.plus(dt);\r
                        \r
-                       log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
+                       //log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
                        \r
-                       log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
+                       //log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
                        \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
-                       log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
+                       Amount<Volume> volumeBurnt;\r
+                       if ( settings.getVolumeMethod() == BurnSettings.BurnVolumeMethod.DeltaVolume ){\r
+                               volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
+                       } else {\r
+                               volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
+                       }\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
-                       log.debug("Mass Gen Rate: " + mGenRate);\r
+                       assert(positive(mGenRate));\r
+                       \r
+                       //log.debug("Mass Gen Rate: " + mGenRate);\r
                        \r
                        //Calculate specific gas constant\r
-                       Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());\r
+                       Amount<?> specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());\r
                        //This unit conversion helps JScience to convert nozzle flow rate to\r
                        //kg/s a little later on I verified the conversion by hand and\r
                        //JScience checks it too.\r
@@ -132,46 +243,72 @@ public class Burn {
                        \r
                        Amount<MassFlowRate> mNozzle;\r
                        {\r
-                               Amount<Pressure> pDiff = prev.chamberPressure.minus(atmosphereicPressure);\r
-                               log.debug("Pdiff: " + pDiff);\r
+                               Amount<Pressure> pDiff = prev.chamberPressure.minus(settings.getAtmosphereicPressure());\r
+                               //log.debug("Pdiff: " + pDiff);\r
                                Amount<Area> aStar = motor.getNozzle().throatArea();\r
                                double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats();\r
                                double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1)));\r
-                               Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
+                               Amount<?> sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
                                mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
-                               log.debug("Mass Exit Rate: " + mNozzle.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
-                       log.debug("Mass Storage Rate: " + massStorageRate);\r
+                       //log.debug("Mass Storage Rate: " + massStorageRate);\r
 \r
                        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.times(settings.getRegStepDecreaseFactor());\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
-                       log.debug("Chamber Product: " + next.chamberProduct);\r
+                       //log.debug("Chamber Product: " + next.chamberProduct);\r
                        \r
                        Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);\r
                        \r
                        log.debug("Product Density: " + combustionProductDensity);\r
                        \r
-                       next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);\r
+                       next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(settings.getAtmosphereicPressure()).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
-                       next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());\r
+                       Amount<Pressure> dp = next.chamberPressure.minus(prev.chamberPressure);\r
+                       if ( dp.abs().isGreaterThan(settings.getChamberPressureMaxDelta())){\r
+                               log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!");\r
+                               regStep = regStep.times(settings.getRegStepDecreaseFactor());\r
+                               continue step;\r
+                       }\r
+                       \r
+                       next.thrust = motor.getNozzle().thrust(next.chamberPressure, settings.getAtmosphereicPressure(), settings.getAtmosphereicPressure(), 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
+                       if ( i > 100 && next.chamberPressure.minus(settings.getAtmosphereicPressure()).abs().isLessThan(settings.getEndPressure())){\r
+                               log.info("Pressure at ~Patm on step " + i);\r
+                               endPressureSteps++;\r
+                               if ( endPressureSteps > 5 )\r
+                                       break;\r
                        }\r
+                       \r
+                       data.put(data.lastKey().plus(dt), next);\r
                }\r
 \r
+               long time = new Date().getTime() - start;\r
+               log.info("Burn took " + time + " millis.");\r
+               done = true;\r
+               for (BurnProgressListener bpl : bpls ){\r
+                       bpl.burnComplete();\r
+               }\r
        }\r
        \r
        @SuppressWarnings("unchecked")\r
@@ -186,15 +323,20 @@ public class Burn {
        }\r
        \r
        public Amount<Pressure> pressure(Amount<Duration> time){\r
-               return data.get(time).chamberPressure;\r
+               return getData().get(time).chamberPressure;\r
        }\r
        \r
        public Amount<Force> thrust(Amount<Duration> time){\r
-               return data.get(time).thrust;\r
+               return getData().get(time).thrust;\r
        }\r
        \r
        public Amount<Dimensionless> kn(Amount<Length> regression){\r
                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