- * The dimensional factor is taken from the {@code problem}, while the discrete
- * pulse width (a multiplier of the {@code grid} parameter {@code tau} is
- * calculated using the {@code gridTime} method.
- *
- *
- * @param problem the problem, used to extract the dimensional time factor
- * @param grid a grid used to discretise the {@code pulse}
- */
-
- public DiscretePulse(Problem problem, Grid grid) {
- timeFactor = problem.getProperties().timeFactor();
- this.grid = grid;
- this.pulse = problem.getPulse();
-
- recalculate();
-
- pulse.getPulseShape().init(this);
- pulse.addListener(e -> {
- recalculate();
- pulse.getPulseShape().init(this);
-
- });
- }
-
- /**
- * Uses the {@code PulseTemporalShape} of the {@code Pulse} object to calculate
- * the laser power at the specified moment of {@code time}.
- *
- * @param time the time argument
- * @return the laser power at the specified moment of {@code time}
- */
-
- public double laserPowerAt(double time) {
- return pulse.evaluateAt(time);
- }
-
- /**
- * Recalculates the {@code discretePulseWidth} by calling {@code gridTime} on
- * the physical pulse width and {@code timeFactor}.
- *
- * @see pulse.problem.schemes.Grid.gridTime(double,double)
- */
-
- public void recalculate() {
- final double width = ((Number) pulse.getPulseWidth().getValue()).doubleValue();
- discretePulseWidth = grid.gridTime(width, timeFactor);
- }
-
- /**
- * Gets the discrete pulse width defined by {@code DiscretePulse}.
- *
- * @return a double, representing the discrete pulse width.
- */
-
- public double getDiscreteWidth() {
- return discretePulseWidth;
- }
-
- /**
- * Gets the physical {@code Pulse}
- *
- * @return the {@code Pulse} object
- */
-
- public Pulse getPulse() {
- return pulse;
- }
-
- /**
- * Gets the {@code Grid} object used to construct this {@code DiscretePulse}
- *
- * @return the {@code Grid} object.
- */
-
- public Grid getGrid() {
- return grid;
- }
-
-}
\ No newline at end of file
+public class DiscretePulse implements Serializable {
+
+ private static final long serialVersionUID = 5826506918603729615L;
+ private final Grid grid;
+ private final Pulse pulse;
+ private final ExperimentalData data;
+
+ private double widthOnGrid;
+ private double characteristicTime;
+ private double invTotalEnergy; //normalisation factor
+
+ /**
+ * This number shows how small the actual pulse may be compared to the
+ * half-time. If the pulse is shorter than
+ *
+ * The dimensional factor is taken from the {@code problem}, while the
+ * discrete pulse width (a multiplier of the {@code grid} parameter
+ * {@code tau} is calculated using the {@code gridTime} method.
+ *
+ *
+ * @param problem the problem, used to extract the dimensional time factor
+ * @param grid a grid used to discretise the {@code pulse}
+ */
+ public DiscretePulse(Problem problem, Grid grid) {
+ this.grid = grid;
+ characteristicTime = problem.getProperties().characteristicTime();
+ this.pulse = problem.getPulse();
+
+ Object ancestor
+ = Objects.requireNonNull(problem.specificAncestor(SearchTask.class),
+ "Problem has not been assigned to a SearchTask");
+
+ data = (ExperimentalData) (((SearchTask) ancestor).getInput());
+ init();
+
+ PropertyHolderListener phl = e -> {
+ characteristicTime = problem.getProperties().characteristicTime();
+ widthOnGrid = 0;
+ init();
+ };
+
+ pulse.addListener(e -> {
+ widthOnGrid = 0;
+ init();
+ });
+ problem.addListener(phl);
+
+ }
+
+ /**
+ * Uses the {@code PulseTemporalShape} of the {@code Pulse} object to
+ * calculate the laser power at the specified moment of {@code time}.
+ *
+ * @param time the time argument
+ * @return the laser power at the specified moment of {@code time}
+ */
+ public double laserPowerAt(double time) {
+ return invTotalEnergy * pulse.getPulseShape().evaluateAt(time);
+ }
+
+ /**
+ * Recalculates the {@code discretePulseWidth} by calling {@code gridTime}
+ * on the physical pulse width and {@code timeFactor}.
+ *
+ * @see pulse.problem.schemes.Grid.gridTime(double,double)
+ */
+ public final void init() {
+ final double nominalWidth = ((Number) pulse.getPulseWidth().getValue()).doubleValue();
+ final double resolvedWidth = resolvedPulseWidthSeconds();
+
+ final double EPS = 1E-10;
+
+ double oldValue = widthOnGrid;
+ this.widthOnGrid = pulseWidthGrid();
+
+ /**
+ * The pulse is too short, which makes calculations too expensive. Can
+ * we replace it with a rectangular pulse shape instead?
+ */
+ if (nominalWidth < resolvedWidth - EPS && oldValue < EPS) {
+ //change shape to rectangular
+ var shape = new RectangularPulse();
+ pulse.setPulseShape(shape);
+ shape.init(null, this);
+ } else {
+ pulse.getPulseShape().init(data, this);
+ }
+
+ invTotalEnergy = 1.0 / totalEnergy();
+ }
+
+ /**
+ * Optimises the {@code Grid} parameters so that the timestep is
+ * sufficiently small to enable accurate pulse correction.
+ *
+ * This can change the {@code tauFactor} and {@code tau} variables in the
+ * {@code Grid} object if {@code discretePulseWidth/(M - 1) < grid.tau},
+ * where M is the required number of pulse calculations.
+ *
+ *
+ * @see PulseTemporalShape.getRequiredDiscretisation()
+ */
+ public double pulseWidthGrid() {
+ //minimum number of points for pulse calculation
+ int reqPoints = pulse.getPulseShape().getRequiredDiscretisation();
+ //physical pulse width in time units
+ double experimentalWidth = (double) pulse.getPulseWidth().getValue();
+
+ //minimum resolved pulse width in time units for that specific problem
+ double resolvedWidth = resolvedPulseWidthSeconds();
+
+ double pWidth = Math.max(experimentalWidth, resolvedWidth);
+
+ final double EPS = 1E-10;
+
+ double newTau = pWidth / characteristicTime / reqPoints;
+
+ double result = 0;
+
+ if (newTau < grid.getTimeStep() - EPS) {
+ double newTauFactor = (double) grid.getTimeFactor().getValue() / 2.0;
+ grid.setTimeFactor(derive(TAU_FACTOR, newTauFactor));
+ result = pulseWidthGrid();
+ } else {
+ result = grid.gridTime(pWidth, characteristicTime);
+ }
+
+ return result;
+ }
+
+ /**
+ * Calculates the total pulse energy using a numerical integrator.The
+ * normalisation factor is then equal to the inverse total energy.
+ *
+ * @return the total pulse energy, assuming sample area fully covered by the
+ * beam
+ */
+ public final double totalEnergy() {
+ var pulseShape = pulse.getPulseShape();
+
+ var integrator = new MidpointIntegrator(new Segment(0, widthOnGrid)) {
+
+ @Override
+ public double integrand(double... vars) {
+ return pulseShape.evaluateAt(vars[0]);
+ }
+
+ };
+
+ return integrator.integrate();
+ }
+
+ /**
+ * Gets the discrete dimensionless pulse width, which is a multiplier of the
+ * current grid timestep. The pulse width is converted to the dimensionless
+ * pulse width by dividing the real value by
- * Calls the constructor of the superclass, after which calculates the
- * {@code discretePulseSpot} using the {@code gridRadialDistance} method of this
- * class. The dimension factor is defined as the sample diameter.
- *
- *
- * @param problem a two-dimensional problem
- * @param grid the two-dimensional grid
- */
+ /**
+ * The constructor for {@code DiscretePulse2D}.
+ *
+ * Calls the constructor of the superclass, after which calculates the
+ * {@code discretePulseSpot} using the {@code gridRadialDistance} method of
+ * this class. The dimension factor is defined as the sample diameter.
+ *
+ *
+ * @param problem a two-dimensional problem
+ * @param grid the two-dimensional grid
+ */
+ public DiscretePulse2D(ClassicalProblem2D problem, Grid2D grid) {
+ super(problem, grid);
+ var properties = (ExtendedThermalProperties) problem.getProperties();
+ calcPulseSpot(properties);
+ properties.addListener(e -> calcPulseSpot(properties));
+ }
- public DiscretePulse2D(ClassicalProblem2D problem, Grid2D grid) {
- super(problem, grid);
- var properties = (ExtendedThermalProperties)problem.getProperties();
- coordFactor = (double) properties.getSampleDiameter().getValue() / 2.0;
- var pulse = (Pulse2D)problem.getPulse();
- discretePulseSpot = grid.gridRadialDistance((double) pulse.getSpotDiameter().getValue() / 2.0, coordFactor);
+ /**
+ * This calculates the dimensionless, discretised pulse function at a
+ * dimensionless radial coordinate {@code coord}.
+ *
+ * It uses a Heaviside function to determine whether the {@code radialCoord}
+ * lies within the {@code 0 <= radialCoord <= discretePulseSpot} interval.
+ * It uses the {@code time} parameter to determine the discrete pulse
+ * function using {@code evaluateAt(time)}.
+ *
+ * @param time the time for calculation
+ * @param radialCoord - the radial coordinate [length dimension]
+ * @return the pulse function at {@code time} and {@code coord}, or 0 if
+ * {@code coord > spotDiameter}.
+ * @see pulse.problem.laser.PulseTemporalShape.laserPowerAt(double)
+ */
+ public double evaluateAt(double time, double radialCoord) {
+ return laserPowerAt(time)
+ * (0.5 + 0.5 * signum(discretePulseSpot - radialCoord));
+ }
- }
+ /**
+ * Calculates the laser power at a give moment in time. The total laser
+ * energy is normalised over a beam partially illuminating the sample
+ * surface.
+ *
+ * @param time a moment in time (in dimensionless units)
+ * @return the laser power in arbitrary units
+ */
+ @Override
+ public double laserPowerAt(double time) {
+ return normFactor * super.laserPowerAt(time);
+ }
- /**
- * This calculates the dimensionless, discretised pulse function at a
- * dimensionless radial coordinate {@code coord}.
- *
- * It uses a Heaviside function to determine whether the {@code radialCoord}
- * lies within the {@code 0 <= radialCoord <= discretePulseSpot} interval. It
- * uses the {@code time} parameter to determine the discrete pulse function
- * using {@code evaluateAt(time)}.
- *
- * @param time the time for calculation
- * @param radialCoord - the radial coordinate [length dimension]
- * @return the pulse function at {@code time} and {@code coord}, or 0 if
- * {@code coord > spotDiameter}.
- * @see pulse.problem.laser.PulseTemporalShape.laserPowerAt(double)
- */
+ private void calcPulseSpot(ExtendedThermalProperties properties) {
+ sampleRadius = (double) properties.getSampleDiameter().getValue() / 2.0;
+ evalPulseSpot();
+ }
- public double evaluateAt(double time, double radialCoord) {
- return laserPowerAt(time) * (0.5 + 0.5 * signum(discretePulseSpot - radialCoord));
- }
+ /**
+ * Calculates the {@code discretePulseSpot} using the
+ * {@code gridRadialDistance} method.
+ *
+ * @see pulse.problem.schemes.Grid2D.gridRadialDistance(double,double)
+ */
+ public final void evalPulseSpot() {
+ var pulse = (Pulse2D) getPhysicalPulse();
+ var grid2d = (Grid2D) getGrid();
+ final double spotRadius = (double) pulse.getSpotDiameter().getValue() / 2.0;
+ discretePulseSpot = grid2d.gridRadialDistance(spotRadius, sampleRadius);
+ grid2d.adjustStepSize(this);
+ normFactor = sampleRadius * sampleRadius / spotRadius / spotRadius;
+ }
- /**
- * Calls the superclass method, then calculates the {@code discretePulseSpot}
- * using the {@code gridRadialDistance} method.
- *
- * @see pulse.problem.schemes.Grid2D.gridRadialDistance(double,double)
- */
+ public final double getDiscretePulseSpot() {
+ return discretePulseSpot;
+ }
- @Override
- public void recalculate() {
- super.recalculate();
- final double radius = (double) ((Pulse2D) getPulse()).getSpotDiameter().getValue() / 2.0;
- discretePulseSpot = ((Grid2D) getGrid()).gridRadialDistance(radius, coordFactor);
- }
+ public final double getRadialConversionFactor() {
+ return sampleRadius;
+ }
- public double getDiscretePulseSpot() {
- return discretePulseSpot;
- }
+ /**
+ * A smaller tolerance factor is set for 2D calculations
+ */
+ @Override
+ public int getWidthToleranceFactor() {
+ return WIDTH_TOLERANCE_FACTOR;
+ }
-}
\ No newline at end of file
+}
diff --git a/src/main/java/pulse/problem/laser/ExponentiallyModifiedGaussian.java b/src/main/java/pulse/problem/laser/ExponentiallyModifiedGaussian.java
index 77f803b0..d4490d2b 100644
--- a/src/main/java/pulse/problem/laser/ExponentiallyModifiedGaussian.java
+++ b/src/main/java/pulse/problem/laser/ExponentiallyModifiedGaussian.java
@@ -6,185 +6,158 @@
import static pulse.properties.NumericProperties.def;
import static pulse.properties.NumericProperties.derive;
import static pulse.properties.NumericProperty.requireType;
-import static pulse.properties.NumericPropertyKeyword.INTEGRATION_SEGMENTS;
import static pulse.properties.NumericPropertyKeyword.SKEW_LAMBDA;
import static pulse.properties.NumericPropertyKeyword.SKEW_MU;
import static pulse.properties.NumericPropertyKeyword.SKEW_SIGMA;
-import java.util.List;
+import java.util.Set;
-import pulse.math.FixedIntervalIntegrator;
-import pulse.math.MidpointIntegrator;
-import pulse.math.Segment;
import pulse.properties.NumericProperty;
import pulse.properties.NumericPropertyKeyword;
-import pulse.properties.Property;
/**
* Represents the exponentially modified Gaussian function, which is given by
* three independent parameters (μ, σ and λ).
- *
+ *
* @see Wikipedia page
*
*/
-
public class ExponentiallyModifiedGaussian extends PulseTemporalShape {
- private double mu;
- private double sigma;
- private double lambda;
- private double norm;
-
- private final static int DEFAULT_POINTS = 100;
- private FixedIntervalIntegrator integrator;
-
- /**
- * Creates an exponentially modified Gaussian with the default parameter values.
- */
-
- public ExponentiallyModifiedGaussian() {
- mu = (double) def(SKEW_MU).getValue();
- lambda = (double) def(SKEW_LAMBDA).getValue();
- sigma = (double) def(SKEW_SIGMA).getValue();
- norm = 1.0;
- integrator = new MidpointIntegrator(new Segment(0.0, getPulseWidth()),
- derive(INTEGRATION_SEGMENTS, DEFAULT_POINTS)) {
-
- @Override
- public double integrand(double... vars) {
- return evaluateAt(vars[0]);
- }
-
- };
- }
-
- /**
- * This calls the superclass {@code init method} and sets the normalisation
- * factor to .
- */
-
- @Override
- public void init(DiscretePulse pulse) {
- super.init(pulse);
- norm = 1.0; // resets the normalisation factor to unity
- norm = 1.0 / area(); // calculates the area. The normalisation factor is then set to the inverse of
- // the area.
- }
-
- /**
- * Uses numeric integration (midpoint rule) to calculate the area of the pulse
- * shape corresponding to the selected parameters.
- *
- * @return the area
- */
-
- private double area() {
- integrator.setBounds(new Segment(0.0, getPulseWidth()));
- return integrator.integrate();
- }
-
- /**
- * Evaluates the laser power function. The error function is calculated using
- * the ApacheCommonsMath library tools.
- *
- * @see https://tinyurl.com/ExpModifiedGaussian
- * @param time is measured from the 'start' of laser pulse
- */
-
- @Override
- public double evaluateAt(double time) {
- final var reducedTime = time / getPulseWidth();
-
- final double lambdaHalf = 0.5 * lambda;
- final double sigmaSq = sigma * sigma;
-
- return norm * lambdaHalf * exp(lambdaHalf * (2.0 * mu + lambda * sigmaSq - 2.0 * reducedTime))
- * erfc((mu + lambda * sigmaSq - reducedTime) / (sqrt(2) * sigma));
-
- }
-
- @Override
- public List listedTypes() {
- var list = super.listedTypes();
- list.add(def(SKEW_MU));
- list.add(def(SKEW_LAMBDA));
- list.add(def(SKEW_SIGMA));
- return list;
- }
-
- /**
- * @return the μ parameter
- */
-
- public NumericProperty getMu() {
- return derive(SKEW_MU, mu);
- }
-
- /**
- * @return the σ parameter
- */
-
- public NumericProperty getSigma() {
- return derive(SKEW_SIGMA, sigma);
- }
-
- /**
- * @return the λ parameter
- */
-
- public NumericProperty getLambda() {
- return derive(SKEW_LAMBDA, lambda);
- }
-
- /**
- * Sets the {@code SKEW_LAMBDA} parameter
- *
- * @param p the λ parameter
- */
-
- public void setLambda(NumericProperty p) {
- requireType(p, SKEW_LAMBDA);
- this.lambda = (double) p.getValue();
- }
-
- /**
- * Sets the {@code SKEW_MU} parameter
- *
- * @param p the μ parameter
- */
-
- public void setMu(NumericProperty p) {
- requireType(p, SKEW_MU);
- this.mu = (double) p.getValue();
- }
-
- /**
- * Sets the {@code SKEW_SIGMA} parameter
- *
- * @param p the σ parameter
- */
-
- public void setSigma(NumericProperty p) {
- requireType(p, SKEW_SIGMA);
- this.sigma = (double) p.getValue();
- }
-
- @Override
- public void set(NumericPropertyKeyword type, NumericProperty property) {
- switch (type) {
- case SKEW_MU:
- setMu(property);
- break;
- case SKEW_LAMBDA:
- setLambda(property);
- break;
- case SKEW_SIGMA:
- setSigma(property);
- break;
- default:
- break;
- }
- firePropertyChanged(this, property);
- }
-
-}
\ No newline at end of file
+ private static final long serialVersionUID = -4437794069527301235L;
+ private double mu;
+ private double sigma;
+ private double lambda;
+
+ private final static int MIN_POINTS = 10;
+
+ /**
+ * Creates an exponentially modified Gaussian with the default parameter
+ * values.
+ */
+ public ExponentiallyModifiedGaussian() {
+ mu = (double) def(SKEW_MU).getValue();
+ lambda = (double) def(SKEW_LAMBDA).getValue();
+ sigma = (double) def(SKEW_SIGMA).getValue();
+ }
+
+ public ExponentiallyModifiedGaussian(ExponentiallyModifiedGaussian another) {
+ super(another);
+ this.mu = another.mu;
+ this.sigma = another.sigma;
+ this.lambda = another.lambda;
+ }
+
+ /**
+ * Evaluates the laser power function. The error function is calculated
+ * using the ApacheCommonsMath library tools.
+ *
+ * @see https://tinyurl.com/ExpModifiedGaussian
+ * @param time is measured from the 'start' of laser pulse
+ */
+ @Override
+ public double evaluateAt(double time) {
+ final var reducedTime = time / getPulseWidth();
+
+ final double lambdaHalf = 0.5 * lambda;
+ final double sigmaSq = sigma * sigma;
+
+ return lambdaHalf * exp(lambdaHalf * (2.0 * mu + lambda * sigmaSq - 2.0 * reducedTime))
+ * erfc((mu + lambda * sigmaSq - reducedTime) / (sqrt(2) * sigma));
+
+ }
+
+ /**
+ * @see pulse.properties.NumericPropertyKeyword.SKEW_MU
+ * @see pulse.properties.NumericPropertyKeyword.SKEW_LAMBDA
+ * @see pulse.properties.NumericPropertyKeyword.SKEW_SIGMA
+ */
+ @Override
+ public Set listedKeywords() {
+ var set = super.listedKeywords();
+ set.add(SKEW_MU);
+ set.add(SKEW_LAMBDA);
+ set.add(SKEW_SIGMA);
+ return set;
+ }
+
+ /**
+ * @return the μ parameter
+ */
+ public NumericProperty getMu() {
+ return derive(SKEW_MU, mu);
+ }
+
+ /**
+ * @return the σ parameter
+ */
+ public NumericProperty getSigma() {
+ return derive(SKEW_SIGMA, sigma);
+ }
+
+ /**
+ * @return the λ parameter
+ */
+ public NumericProperty getLambda() {
+ return derive(SKEW_LAMBDA, lambda);
+ }
+
+ /**
+ * Sets the {@code SKEW_LAMBDA} parameter
+ *
+ * @param p the λ parameter
+ */
+ public void setLambda(NumericProperty p) {
+ requireType(p, SKEW_LAMBDA);
+ this.lambda = (double) p.getValue();
+ }
+
+ /**
+ * Sets the {@code SKEW_MU} parameter
+ *
+ * @param p the μ parameter
+ */
+ public void setMu(NumericProperty p) {
+ requireType(p, SKEW_MU);
+ this.mu = (double) p.getValue();
+ }
+
+ /**
+ * Sets the {@code SKEW_SIGMA} parameter
+ *
+ * @param p the σ parameter
+ */
+ public void setSigma(NumericProperty p) {
+ requireType(p, SKEW_SIGMA);
+ this.sigma = (double) p.getValue();
+ }
+
+ @Override
+ public void set(NumericPropertyKeyword type, NumericProperty property) {
+ switch (type) {
+ case SKEW_MU:
+ setMu(property);
+ break;
+ case SKEW_LAMBDA:
+ setLambda(property);
+ break;
+ case SKEW_SIGMA:
+ setSigma(property);
+ break;
+ default:
+ break;
+ }
+ firePropertyChanged(this, property);
+ }
+
+ @Override
+ public PulseTemporalShape copy() {
+ return new ExponentiallyModifiedGaussian(this);
+ }
+
+ @Override
+ public int getRequiredDiscretisation() {
+ return MIN_POINTS;
+ }
+
+}
diff --git a/src/main/java/pulse/problem/laser/NumericPulse.java b/src/main/java/pulse/problem/laser/NumericPulse.java
new file mode 100644
index 00000000..a71d65fd
--- /dev/null
+++ b/src/main/java/pulse/problem/laser/NumericPulse.java
@@ -0,0 +1,176 @@
+package pulse.problem.laser;
+
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
+import static pulse.properties.NumericProperties.derive;
+import static pulse.properties.NumericPropertyKeyword.PULSE_WIDTH;
+
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.analysis.interpolation.AkimaSplineInterpolator;
+import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
+
+import pulse.input.ExperimentalData;
+import pulse.problem.statements.Problem;
+import pulse.properties.NumericProperty;
+import pulse.properties.NumericPropertyKeyword;
+import pulse.tasks.SearchTask;
+
+import pulse.baseline.FlatBaseline;
+import pulse.tasks.Calculation;
+import pulse.util.FunctionSerializer;
+
+/**
+ * A numeric pulse is given by a set of discrete {@code NumericPulseData}
+ * measured independently using a pulse diode.
+ *
+ * @see pulse.problem.laser.NumericPulseData
+ *
+ */
+public class NumericPulse extends PulseTemporalShape {
+
+ private static final long serialVersionUID = 6088261629992349844L;
+ private NumericPulseData pulseData;
+ private transient UnivariateFunction interpolation;
+
+ private final static int MIN_POINTS = 20;
+
+ public NumericPulse() {
+ //intentionally blank
+ }
+
+ /**
+ * Copy constructor
+ *
+ * @param pulse another numeric pulse, the data of which will be copied
+ */
+ public NumericPulse(NumericPulse pulse) {
+ super(pulse);
+ this.pulseData = new NumericPulseData(pulseData);
+ }
+
+ /**
+ * Defines the pulse width as the last element of the time sequence
+ * contained in {@code NumericPulseData}. Calls {@code super.init}, then
+ * interpolates the input pulse using spline functions and normalises the
+ * output.
+ *
+ * @param data
+ * @see normalise()
+ *
+ */
+ @Override
+ public void init(ExperimentalData data, DiscretePulse pulse) {
+ //generate baseline-subtracted numeric data from ExperimentalData
+ baselineSubtractedFrom(data);
+
+ //notify host pulse object of a new pulse width
+ var problem = ((Calculation) ((SearchTask) data.getParent())
+ .getResponse()).getProblem();
+ setPulseWidthOf(problem);
+
+ //convert to dimensionless time and interpolate
+ double timeFactor = problem.getProperties().characteristicTime();
+ doInterpolation(timeFactor);
+ }
+
+ /**
+ * Copies the numeric pulse from metadata and subtracts a horizontal
+ * baseline from the data points assigned to {@code pulseData}.
+ *
+ * @param data the experimental data containing the metadata with numeric
+ * pulse data.
+ */
+ private void baselineSubtractedFrom(ExperimentalData data) {
+ pulseData = new NumericPulseData(data.getMetadata().getPulseData());
+
+ //subtracts a horizontal baseline from the pulse data
+ var baseline = new FlatBaseline();
+ baseline.fitTo(pulseData);
+
+ for (int i = 0, size = pulseData.getTimeSequence().size(); i < size; i++) {
+ pulseData.setSignalAt(i,
+ pulseData.signalAt(i) - baseline.valueAt(pulseData.timeAt(i)));
+ }
+ }
+
+ private void setPulseWidthOf(Problem problem) {
+ var timeSequence = pulseData.getTimeSequence();
+ double pulseWidth = timeSequence.get(timeSequence.size() - 1);
+
+ var pulseObject = problem.getPulse();
+ pulseObject.setPulseWidth(derive(PULSE_WIDTH, pulseWidth));
+
+ }
+
+ private void doInterpolation(double timeFactor) {
+ var interpolator = new AkimaSplineInterpolator();
+
+ var timeList = pulseData.getTimeSequence().stream().mapToDouble(d -> d / timeFactor).toArray();
+ var powerList = pulseData.getSignalData();
+
+ this.setPulseWidth(timeList[timeList.length - 1]);
+
+ interpolation = interpolator.interpolate(timeList,
+ powerList.stream().mapToDouble(d -> d).toArray());
+ }
+
+ /**
+ * If the argument is less than the pulse width, uses a spline to
+ * interpolate the pulse data at {@code time}. Otherwise returns zero.
+ */
+ @Override
+ public double evaluateAt(double time) {
+ return time > getPulseWidth() ? 0.0 : interpolation.value(time);
+ }
+
+ @Override
+ public PulseTemporalShape copy() {
+ return new NumericPulse();
+ }
+
+ /**
+ * Does not define any property.
+ */
+ @Override
+ public void set(NumericPropertyKeyword type, NumericProperty property) {
+ // TODO Auto-generated method stub
+ }
+
+ public NumericPulseData getData() {
+ return pulseData;
+ }
+
+ public void setData(NumericPulseData pulseData) {
+ this.pulseData = pulseData;
+
+ }
+
+ public UnivariateFunction getInterpolation() {
+ return interpolation;
+ }
+
+ @Override
+ public int getRequiredDiscretisation() {
+ return MIN_POINTS;
+ }
+
+ /*
+ Serialization
+ */
+ private void writeObject(ObjectOutputStream oos)
+ throws IOException {
+ // default serialization
+ oos.defaultWriteObject();
+ // write the object
+ FunctionSerializer.writeSplineFunction((PolynomialSplineFunction) interpolation, oos);
+ }
+
+ private void readObject(ObjectInputStream ois)
+ throws ClassNotFoundException, IOException {
+ // default deserialization
+ ois.defaultReadObject();
+ this.interpolation = FunctionSerializer.readSplineFunction(ois);
+ }
+
+}
diff --git a/src/main/java/pulse/problem/laser/NumericPulseData.java b/src/main/java/pulse/problem/laser/NumericPulseData.java
new file mode 100644
index 00000000..7bbfeeed
--- /dev/null
+++ b/src/main/java/pulse/problem/laser/NumericPulseData.java
@@ -0,0 +1,78 @@
+package pulse.problem.laser;
+
+import java.util.List;
+import pulse.AbstractData;
+import pulse.DiscreteInput;
+import pulse.input.IndexRange;
+import pulse.input.Range;
+
+/**
+ * An instance of the {@code AbstractData} class, which also declares an
+ * {@code externalID}. Use to store numeric data of the pulse for each
+ * measurement imported from an external source.
+ *
+ */
+public class NumericPulseData extends AbstractData implements DiscreteInput {
+
+ private static final long serialVersionUID = 8142129124831241206L;
+ private final int externalID;
+
+ /**
+ * Stores {@code id} and calls super-constructor
+ *
+ * @param id an external ID defined in the imported file
+ */
+ public NumericPulseData(int id) {
+ super();
+ this.externalID = id;
+ }
+
+ /**
+ * Copies everything, including the id number.
+ *
+ * @param data another object
+ */
+ public NumericPulseData(NumericPulseData data) {
+ super(data);
+ this.externalID = data.externalID;
+ }
+
+ /**
+ * Adds a data point to the internal storage and increments counter.
+ */
+ @Override
+ public void addPoint(double time, double power) {
+ super.addPoint(time, power);
+ super.incrementCount();
+ }
+
+ /**
+ * Gets the external ID usually specified in the experimental files. Note
+ * this is not a {@code NumericProperty}
+ *
+ * @return an integer, representing the external ID
+ */
+ public int getExternalID() {
+ return externalID;
+ }
+
+ public double pulseWidth() {
+ return super.timeLimit();
+ }
+
+ @Override
+ public List getX() {
+ return getTimeSequence();
+ }
+
+ @Override
+ public List getY() {
+ return getSignalData();
+ }
+
+ @Override
+ public IndexRange getIndexRange() {
+ return new IndexRange(this.getTimeSequence(), Range.UNLIMITED);
+ }
+
+}
diff --git a/src/main/java/pulse/problem/laser/PulseTemporalShape.java b/src/main/java/pulse/problem/laser/PulseTemporalShape.java
index ed24ec38..9a74c09c 100644
--- a/src/main/java/pulse/problem/laser/PulseTemporalShape.java
+++ b/src/main/java/pulse/problem/laser/PulseTemporalShape.java
@@ -1,50 +1,69 @@
package pulse.problem.laser;
+import pulse.input.ExperimentalData;
import pulse.util.PropertyHolder;
import pulse.util.Reflexive;
/**
* An abstract time-dependent pulse shape. Declares the abstract method to
* calculate the laser power function at a given moment of time. This generally
- * utilises a discrete pulse width.
+ * utilises a discrete pulse width. By default, uses a midpoint-rule numeric
+ * integrator to calculate the pulse integral.
*
*/
-
public abstract class PulseTemporalShape extends PropertyHolder implements Reflexive {
- private double width;
+ private double width;
+
+ public PulseTemporalShape() {
+ //intentionlly blank
+ }
+
+ public PulseTemporalShape(PulseTemporalShape another) {
+ this.width = another.width;
+ }
+
+ /**
+ * This evaluates the dimensionless, discretised pulse function on a
+ * {@code grid} needed to evaluate the heat source in the difference scheme.
+ *
+ * @param time the dimensionless time (a multiplier of {@code tau}), at
+ * which calculation should be performed
+ * @return a double value, representing the pulse function at {@code time}
+ */
+ public abstract double evaluateAt(double time);
- /**
- * This evaluates the dimensionless, discretised pulse function on a
- * {@code grid} needed to evaluate the heat source in the difference scheme.
- *
- * @param time the dimensionless time (a multiplier of {@code tau}), at which
- * calculation should be performed
- * @return a double value, representing the pulse function at {@code time}
- */
+ /**
+ * Stores the pulse width from {@code pulse} and initialises area
+ * integration.
+ *
+ * @param data
+ * @param pulse the discrete pulse containing the pulse width
+ */
+ public void init(ExperimentalData data, DiscretePulse pulse) {
+ width = pulse.getDiscreteWidth();
+ }
- public abstract double evaluateAt(double time);
+ public abstract PulseTemporalShape copy();
- public void init(DiscretePulse pulse) {
- width = pulse.getDiscreteWidth();
- }
+ @Override
+ public String getPrefix() {
+ return "Pulse temporal shape";
+ }
- @Override
- public String getPrefix() {
- return "Pulse temporal shape";
- }
+ @Override
+ public String toString() {
+ return getClass().getSimpleName();
+ }
- @Override
- public String toString() {
- return getClass().getSimpleName();
- }
+ public double getPulseWidth() {
+ return width;
+ }
- public double getPulseWidth() {
- return width;
- }
+ public void setPulseWidth(double width) {
+ this.width = width;
+ }
- public void setPulseWidth(double width) {
- this.width = width;
- }
+ public abstract int getRequiredDiscretisation();
-}
\ No newline at end of file
+}
diff --git a/src/main/java/pulse/problem/laser/RectangularPulse.java b/src/main/java/pulse/problem/laser/RectangularPulse.java
index fa336720..7bc5dfbd 100644
--- a/src/main/java/pulse/problem/laser/RectangularPulse.java
+++ b/src/main/java/pulse/problem/laser/RectangularPulse.java
@@ -9,25 +9,37 @@
* The simplest pulse shape defined as , where
* is the signum function, pulse is the pulse width.
- *
+ *
* @see java.lang.Math.signum(double)
*/
-
public class RectangularPulse extends PulseTemporalShape {
- /**
- * @param time the time measured from the start of the laser pulse.
- */
-
- @Override
- public double evaluateAt(double time) {
- var width = getPulseWidth();
- return 0.5 / width * (1 + signum(width - time));
- }
-
- @Override
- public void set(NumericPropertyKeyword type, NumericProperty property) {
- // intentionally blak
- }
-
-}
\ No newline at end of file
+ private static final long serialVersionUID = 8207478409316696745L;
+ private final static int MIN_POINTS = 2;
+
+ /**
+ * @param time the time measured from the start of the laser pulse.
+ * @return
+ */
+ @Override
+ public double evaluateAt(double time) {
+ var width = getPulseWidth();
+ return 0.5 / width * (1 + signum(width - time));
+ }
+
+ @Override
+ public void set(NumericPropertyKeyword type, NumericProperty property) {
+ // intentionally blak
+ }
+
+ @Override
+ public PulseTemporalShape copy() {
+ return new RectangularPulse();
+ }
+
+ @Override
+ public int getRequiredDiscretisation() {
+ return MIN_POINTS;
+ }
+
+}
diff --git a/src/main/java/pulse/problem/laser/TrapezoidalPulse.java b/src/main/java/pulse/problem/laser/TrapezoidalPulse.java
index 570af2f9..908e2cfc 100644
--- a/src/main/java/pulse/problem/laser/TrapezoidalPulse.java
+++ b/src/main/java/pulse/problem/laser/TrapezoidalPulse.java
@@ -6,122 +6,130 @@
import static pulse.properties.NumericPropertyKeyword.TRAPEZOIDAL_FALL_PERCENTAGE;
import static pulse.properties.NumericPropertyKeyword.TRAPEZOIDAL_RISE_PERCENTAGE;
-import java.util.List;
+import java.util.Set;
import pulse.properties.NumericProperty;
import pulse.properties.NumericPropertyKeyword;
-import pulse.properties.Property;
/**
* A trapezoidal pulse shape, which combines a rise segment, a constant-power
* segment, and a fall segment. The rise and fall ratios can be changed.
*
*/
-
public class TrapezoidalPulse extends PulseTemporalShape {
- private double rise;
- private double fall;
- private double h;
-
- /**
- * Constructs a trapezoidal pulse using a default segmentation principle. The
- * reader is referred to the {@code .xml} file containing the default values of
- * {@code TRAPEZOIDAL_RISE_PERCENTAGE} and {@code TRAPEZOIDAL_FALL_PERCENTAGE}.
- * The maximum laser power is adjusted to ensure the area of the shape is equal
- * to unity.
- */
-
- public TrapezoidalPulse() {
- rise = (int) def(TRAPEZOIDAL_RISE_PERCENTAGE).getValue() / 100.0;
- fall = (int) def(TRAPEZOIDAL_FALL_PERCENTAGE).getValue() / 100.0;
- h = height();
- }
-
- @Override
- public void init(DiscretePulse pulse) {
- super.init(pulse);
- h = height();
- }
-
- /**
- * Calculates the height of the trapezium which under current segmentation will
- * yield an area of unity.
- *
- * @return the calculated height of the constant segmment
- */
-
- private double height() {
- return 2.0 / (getPulseWidth() * (2.0 - rise - fall));
- }
-
- /**
- * Calculates power using a piecewise function, which corresponds to either a
- * linearly changing, a constant laser power or zero.
- *
- * @param time the time measured from the start of the laser pulse.
- */
-
- @Override
- public double evaluateAt(double time) {
- final var reducedTime = time / getPulseWidth();
-
- double result = 0;
-
- if (reducedTime < rise) { // triangular
- result = reducedTime * h / rise;
- } else if (reducedTime < 1.0 - fall) { // rectangular
- result = h;
- } else if (reducedTime < 1.0) { // triangular
- final var t2 = (reducedTime - (1.0 - fall));
- result = (fall - t2) * h / fall;
- }
-
- return result;
-
- }
-
- @Override
- public List listedTypes() {
- var list = super.listedTypes();
- list.add(def(TRAPEZOIDAL_RISE_PERCENTAGE));
- list.add(def(TRAPEZOIDAL_FALL_PERCENTAGE));
- return list;
- }
-
- public NumericProperty getTrapezoidalRise() {
- return derive(TRAPEZOIDAL_RISE_PERCENTAGE, (int) (rise * 100));
- }
-
- public NumericProperty getTrapezoidalFall() {
- return derive(TRAPEZOIDAL_FALL_PERCENTAGE, (int) (fall * 100));
- }
-
- public void setTrapezoidalRise(NumericProperty p) {
- requireType(p, TRAPEZOIDAL_RISE_PERCENTAGE);
- this.rise = (int) p.getValue() / 100.0;
- h = height();
- }
-
- public void setTrapezoidalFall(NumericProperty p) {
- requireType(p, TRAPEZOIDAL_FALL_PERCENTAGE);
- this.fall = (int) p.getValue() / 100.0;
- h = height();
- }
-
- @Override
- public void set(NumericPropertyKeyword type, NumericProperty property) {
- switch (type) {
- case TRAPEZOIDAL_RISE_PERCENTAGE:
- setTrapezoidalRise(property);
- break;
- case TRAPEZOIDAL_FALL_PERCENTAGE:
- setTrapezoidalFall(property);
- break;
- default:
- break;
- }
- firePropertyChanged(this, property);
- }
-
-}
\ No newline at end of file
+ private static final long serialVersionUID = 2089809680713225034L;
+ private double rise;
+ private double fall;
+ private double h;
+
+ private final static int MIN_POINTS = 8;
+
+ /**
+ * Constructs a trapezoidal pulse using a default segmentation principle.
+ * The reader is referred to the {@code .xml} file containing the default
+ * values of {@code TRAPEZOIDAL_RISE_PERCENTAGE} and
+ * {@code TRAPEZOIDAL_FALL_PERCENTAGE}. The maximum laser power is adjusted
+ * to ensure the area of the shape is equal to unity.
+ */
+ public TrapezoidalPulse() {
+ rise = (int) def(TRAPEZOIDAL_RISE_PERCENTAGE).getValue() / 100.0;
+ fall = (int) def(TRAPEZOIDAL_FALL_PERCENTAGE).getValue() / 100.0;
+ h = 1.0;
+ }
+
+ public TrapezoidalPulse(TrapezoidalPulse another) {
+ this.rise = another.rise;
+ this.fall = another.fall;
+ this.h = another.h;
+ }
+
+ /**
+ * Calculates the height of the trapezium which under current segmentation
+ * will yield an area of unity.
+ *
+ * @return the calculated height of the constant segmment
+ */
+ private double height() {
+ return 2.0 / (getPulseWidth() * (2.0 - rise - fall));
+ }
+
+ /**
+ * Calculates power using a piecewise function, which corresponds to either
+ * a linearly changing, a constant laser power or zero.
+ *
+ * @param time the time measured from the start of the laser pulse.
+ */
+ @Override
+ public double evaluateAt(double time) {
+ final var reducedTime = time / getPulseWidth();
+
+ double result = 0;
+
+ if (reducedTime < rise) { // triangular
+ result = reducedTime * h / rise;
+ } else if (reducedTime < 1.0 - fall) { // rectangular
+ result = h;
+ } else if (reducedTime < 1.0) { // triangular
+ final var t2 = (reducedTime - (1.0 - fall));
+ result = (fall - t2) * h / fall;
+ }
+
+ return result;
+
+ }
+
+ @Override
+ public Set listedKeywords() {
+ var set = super.listedKeywords();
+ set.add(TRAPEZOIDAL_RISE_PERCENTAGE);
+ set.add(TRAPEZOIDAL_FALL_PERCENTAGE);
+ return set;
+ }
+
+ public NumericProperty getTrapezoidalRise() {
+ return derive(TRAPEZOIDAL_RISE_PERCENTAGE, (int) (rise * 100));
+ }
+
+ public NumericProperty getTrapezoidalFall() {
+ return derive(TRAPEZOIDAL_FALL_PERCENTAGE, (int) (fall * 100));
+ }
+
+ public void setTrapezoidalRise(NumericProperty p) {
+ requireType(p, TRAPEZOIDAL_RISE_PERCENTAGE);
+ this.rise = (int) p.getValue() / 100.0;
+ h = height();
+ }
+
+ public void setTrapezoidalFall(NumericProperty p) {
+ requireType(p, TRAPEZOIDAL_FALL_PERCENTAGE);
+ this.fall = (int) p.getValue() / 100.0;
+ h = height();
+ }
+
+ @Override
+ public void set(NumericPropertyKeyword type, NumericProperty property) {
+ switch (type) {
+ case TRAPEZOIDAL_RISE_PERCENTAGE:
+ setTrapezoidalRise(property);
+ break;
+ case TRAPEZOIDAL_FALL_PERCENTAGE:
+ setTrapezoidalFall(property);
+ break;
+ default:
+ break;
+ }
+ firePropertyChanged(this, property);
+ }
+
+ @Override
+ public PulseTemporalShape copy() {
+ return new TrapezoidalPulse(this);
+ }
+
+ @Override
+ public int getRequiredDiscretisation() {
+ return MIN_POINTS;
+ }
+
+}
diff --git a/src/main/java/pulse/problem/laser/package-info.java b/src/main/java/pulse/problem/laser/package-info.java
index 21573e60..3e4cf687 100644
--- a/src/main/java/pulse/problem/laser/package-info.java
+++ b/src/main/java/pulse/problem/laser/package-info.java
@@ -1,5 +1,4 @@
/**
* This package deals with discrete laser pulse representation and their various temporal shapes.
*/
-
-package pulse.problem.laser;
\ No newline at end of file
+package pulse.problem.laser;
diff --git a/src/main/java/pulse/problem/schemes/ADIScheme.java b/src/main/java/pulse/problem/schemes/ADIScheme.java
index 75e989c1..b20b34db 100644
--- a/src/main/java/pulse/problem/schemes/ADIScheme.java
+++ b/src/main/java/pulse/problem/schemes/ADIScheme.java
@@ -12,54 +12,66 @@
* needed to solve a {@code Problem}.
*
*/
-
public abstract class ADIScheme extends DifferenceScheme {
- /**
- * Creates a new {@code ADIScheme} with default values of grid density and time
- * factor.
- */
-
- public ADIScheme() {
- this(derive(GRID_DENSITY, 30), derive(TAU_FACTOR, 1.0));
- }
-
- /**
- * Creates an {@code ADIScheme} with the specified arguments. This creates an
- * associated {@code Grid2D} object.
- *
- * @param N the grid density
- * @param timeFactor the time factor (τF)
- */
+ /**
+ *
+ */
+ private static final long serialVersionUID = 4772650159522354367L;
- public ADIScheme(NumericProperty N, NumericProperty timeFactor) {
- super();
- setGrid(new Grid2D(N, timeFactor));
- }
+ /**
+ * Creates a new {@code ADIScheme} with default values of grid density and
+ * time factor.
+ */
+ public ADIScheme() {
+ this(derive(GRID_DENSITY, 30), derive(TAU_FACTOR, 0.5));
+ }
- /**
- * Creates an {@code ADIScheme} with the specified arguments. This creates an
- * associated {@code Grid2D} object.
- *
- * @param N the grid density
- * @param timeFactor the time factor (τF)
- * @param timeLimit a custom time limit (tlim)
- */
+ /**
+ * Creates an {@code ADIScheme} with the specified arguments. This creates
+ * an associated {@code Grid2D} object.
+ *
+ * @param N the grid density
+ * @param timeFactor the time factor (τF)
+ */
+ public ADIScheme(NumericProperty N, NumericProperty timeFactor) {
+ super();
+ setGrid(new Grid2D(N, timeFactor));
+ }
- public ADIScheme(NumericProperty N, NumericProperty timeFactor, NumericProperty timeLimit) {
- setTimeLimit(timeLimit);
- setGrid(new Grid2D(N, timeFactor));
- }
+ /**
+ * Creates an {@code ADIScheme} with the specified arguments. This creates
+ * an associated {@code Grid2D} object.
+ *
+ * @param N the grid density
+ * @param timeFactor the time factor (τF)
+ * @param timeLimit a custom time limit (tlim)
+ */
+ public ADIScheme(NumericProperty N, NumericProperty timeFactor, NumericProperty timeLimit) {
+ setTimeLimit(timeLimit);
+ setGrid(new Grid2D(N, timeFactor));
+ }
- /**
- * Prints out the description of this problem type.
- *
- * @return a verbose description of the problem.
- */
+ /**
+ * Prints out the description of this problem type.
+ *
+ * @return a verbose description of the problem.
+ */
+ @Override
+ public String toString() {
+ return getString("ADIScheme.4");
+ }
- @Override
- public String toString() {
- return getString("ADIScheme.4");
- }
+ /**
+ * Contains only an empty statement, as the pulse needs to be calculated not
+ * only for the time step {@code m} but also accounting for the radial
+ * coordinate
+ *
+ * @param m thte time step
+ */
+ @Override
+ public void prepareStep(int m) {
+ //do nothing
+ }
-}
\ No newline at end of file
+}
diff --git a/src/main/java/pulse/problem/schemes/BlockMatrixAlgorithm.java b/src/main/java/pulse/problem/schemes/BlockMatrixAlgorithm.java
index a1530d34..450db9db 100644
--- a/src/main/java/pulse/problem/schemes/BlockMatrixAlgorithm.java
+++ b/src/main/java/pulse/problem/schemes/BlockMatrixAlgorithm.java
@@ -1,79 +1,85 @@
package pulse.problem.schemes;
/**
- * A modification of the algorithm for solving a system of linear equations, where the first and last equation contains references
- * to the last and first elements of the solution, respectively. The corresponding matrix is composed of an inner tridiagonal block
- * and a border formed by an extra row and column. This block system is solved using the Sherman-Morrison-Woodbury identity and the
- * Thomas algorithm for the main block.
+ * A modification of the algorithm for solving a system of linear equations,
+ * where the first and last equation contains references to the last and first
+ * elements of the solution, respectively. The corresponding matrix is composed
+ * of an inner tridiagonal block and a border formed by an extra row and column.
+ * This block system is solved using the Sherman-Morrison-Woodbury identity and
+ * the Thomas algorithm for the main block.
*
*/
-
public class BlockMatrixAlgorithm extends TridiagonalMatrixAlgorithm {
- private double[] gamma;
- private double[] p;
- private double[] q;
-
- public BlockMatrixAlgorithm(Grid grid) {
- super(grid);
- gamma = new double[getAlpha().length];
- p = new double[gamma.length - 1];
- q = new double[gamma.length - 1];
- }
-
- @Override
- public void sweep(double[] V) {
- final int N = V.length - 1;
- for (int j = N - 1; j >= 0; j--)
- V[j] = p[j] + V[N] * q[j];
- }
-
- @Override
- public void evaluateBeta(final double[] U) {
- super.evaluateBeta(U);
- final int N = getGrid().getGridDensityValue();
- var alpha = getAlpha();
- var beta = getBeta();
-
- p[N - 1] = beta[N];
- q[N - 1] = alpha[N] + gamma[N];
-
- for (int i = N - 2; i >= 0; i--) {
- p[i] = alpha[i + 1] * p[i + 1] + beta[i + 1];
- q[i] = alpha[i + 1] * q[i + 1] + gamma[i + 1];
- }
- }
-
- @Override
- public void evaluateBeta(final double[] U, final int start, final int endExclusive) {
- var alpha = getAlpha();
- var grid = getGrid();
- final double HX2_TAU = grid.getXStep() * grid.getXStep() / getGrid().getTimeStep();
-
- final double a = getCoefA();
- final double b = getCoefB();
-
- for (int i = start; i < endExclusive; i++) {
- setBeta(i, beta(U[i - 1] * HX2_TAU, phi(i - 1), i));
- setGamma(i, a * gamma[i - 1] / (b - a * alpha[i - 1]));
- }
-
- }
-
- public double[] getP() {
- return p;
- }
-
- public double[] getQ() {
- return q;
- }
-
- public void setGamma(final int i, final double g) {
- this.gamma[i] = g;
- }
-
- public double[] getGamma() {
- return gamma;
- }
+ private static final long serialVersionUID = -6553638438386098008L;
+ private final double[] gamma;
+ private final double[] p;
+ private final double[] q;
+
+ public BlockMatrixAlgorithm(Grid grid) {
+ super(grid);
+ final int N = this.getGridPoints();
+ gamma = new double[N + 2];
+ p = new double[N];
+ q = new double[N];
+ }
+
+ @Override
+ public void sweep(double[] V) {
+ final int N = V.length - 1;
+ for (int j = N - 1; j >= 0; j--) {
+ V[j] = p[j] + V[N] * q[j];
+ }
+ }
+
+ @Override
+ public void evaluateBeta(final double[] U) {
+ super.evaluateBeta(U);
+ var alpha = getAlpha();
+ var beta = getBeta();
+
+ final int N = getGridPoints();
+
+ p[N - 1] = beta[N];
+ q[N - 1] = alpha[N] + gamma[N];
+
+ for (int i = N - 2; i >= 0; i--) {
+ p[i] = alpha[i + 1] * p[i + 1] + beta[i + 1];
+ q[i] = alpha[i + 1] * q[i + 1] + gamma[i + 1];
+ }
+ }
+
+ @Override
+ public void evaluateBeta(final double[] U, final int start, final int endExclusive) {
+ var alpha = getAlpha();
+
+ final double h = this.getGridStep();
+ final double HX2_TAU = h * h / this.getTimeStep();
+
+ final double a = getCoefA();
+ final double b = getCoefB();
+
+ for (int i = start; i < endExclusive; i++) {
+ setBeta(i, beta(U[i - 1] * HX2_TAU, phi(i - 1), i));
+ setGamma(i, a * gamma[i - 1] / (b - a * alpha[i - 1]));
+ }
+
+ }
+
+ public double[] getP() {
+ return p;
+ }
+
+ public double[] getQ() {
+ return q;
+ }
+
+ public void setGamma(final int i, final double g) {
+ this.gamma[i] = g;
+ }
+
+ public double[] getGamma() {
+ return gamma;
+ }
}
diff --git a/src/main/java/pulse/problem/schemes/CoupledImplicitScheme.java b/src/main/java/pulse/problem/schemes/CoupledImplicitScheme.java
index 06ff68bd..17976a27 100644
--- a/src/main/java/pulse/problem/schemes/CoupledImplicitScheme.java
+++ b/src/main/java/pulse/problem/schemes/CoupledImplicitScheme.java
@@ -1,113 +1,90 @@
package pulse.problem.schemes;
-import static pulse.properties.NumericProperties.def;
-import static pulse.properties.NumericProperties.derive;
import static pulse.properties.NumericPropertyKeyword.NONLINEAR_PRECISION;
-import java.util.List;
+import java.util.Set;
import pulse.problem.schemes.rte.RTECalculationStatus;
+import pulse.problem.schemes.solvers.SolverException;
+import static pulse.problem.schemes.solvers.SolverException.SolverExceptionType.RTE_SOLVER_ERROR;
import pulse.problem.statements.ParticipatingMedium;
import pulse.problem.statements.Problem;
import pulse.properties.NumericProperty;
import pulse.properties.NumericPropertyKeyword;
-import pulse.properties.Property;
-
-public abstract class CoupledImplicitScheme extends ImplicitScheme implements FixedPointIterations {
-
- private RadiativeTransferCoupling coupling;
- private RTECalculationStatus calculationStatus;
- private double nonlinearPrecision;
-
- private double pls;
-
- public CoupledImplicitScheme(NumericProperty N, NumericProperty timeFactor) {
- super();
- setGrid(new Grid(N, timeFactor));
- nonlinearPrecision = (double) def(NONLINEAR_PRECISION).getValue();
- setCoupling(new RadiativeTransferCoupling());
- calculationStatus = RTECalculationStatus.NORMAL;
- }
-
- public CoupledImplicitScheme(NumericProperty N, NumericProperty timeFactor, NumericProperty timeLimit) {
- this(N, timeFactor);
- setTimeLimit(timeLimit);
- }
-
- @Override
- public void timeStep(final int m) {
- pls = pulse(m);
- doIterations(getCurrentSolution(), nonlinearPrecision, m);
- }
-
- @Override
- public void iteration(final int m) {
- super.timeStep(m);
- }
-
- @Override
- public void finaliseIteration(double[] V) {
- setCalculationStatus( coupling.getRadiativeTransferEquation().compute(V) );
- }
-
- public RadiativeTransferCoupling getCoupling() {
- return coupling;
- }
-
- public void setCoupling(RadiativeTransferCoupling coupling) {
- this.coupling = coupling;
- this.coupling.setParent(this);
- }
-
- @Override
- public void finaliseStep() {
- super.finaliseStep();
- coupling.getRadiativeTransferEquation().getFluxes().store();
- }
-
- @Override
- public List listedTypes() {
- List list = super.listedTypes();
- list.add(def(NONLINEAR_PRECISION));
- return list;
- }
-
- public NumericProperty getNonlinearPrecision() {
- return derive(NONLINEAR_PRECISION, nonlinearPrecision);
- }
-
- public void setNonlinearPrecision(NumericProperty nonlinearPrecision) {
- this.nonlinearPrecision = (double) nonlinearPrecision.getValue();
- }
-
- @Override
- public Class domain() {
- return ParticipatingMedium.class;
- }
-
- @Override
- public boolean normalOperation() {
- return super.normalOperation() && (getCalculationStatus() == RTECalculationStatus.NORMAL);
- }
-
- @Override
- public void set(NumericPropertyKeyword type, NumericProperty property) {
- if (type == NONLINEAR_PRECISION) {
- setNonlinearPrecision(property);
- } else
- super.set(type, property);
- }
-
- public RTECalculationStatus getCalculationStatus() {
- return calculationStatus;
- }
-
- public void setCalculationStatus(RTECalculationStatus calculationStatus) {
- this.calculationStatus = calculationStatus;
- }
-
- public double getCurrentPulseValue() {
- return pls;
- }
-
-}
\ No newline at end of file
+
+public abstract class CoupledImplicitScheme extends ImplicitScheme {
+
+ private static final long serialVersionUID = 1974655675470727643L;
+ private RadiativeTransferCoupling coupling;
+ private RTECalculationStatus calculationStatus;
+ private boolean autoUpdateFluxes = true; //should be false for nonlinear solvers
+
+ public CoupledImplicitScheme(NumericProperty N, NumericProperty timeFactor) {
+ super();
+ setGrid(new Grid(N, timeFactor));
+ setCoupling(new RadiativeTransferCoupling());
+ calculationStatus = RTECalculationStatus.NORMAL;
+ }
+
+ public CoupledImplicitScheme(NumericProperty N, NumericProperty timeFactor, NumericProperty timeLimit) {
+ this(N, timeFactor);
+ setTimeLimit(timeLimit);
+ }
+
+ @Override
+ public void finaliseStep() throws SolverException {
+ super.finaliseStep();
+ if (autoUpdateFluxes) {
+ var rte = this.getCoupling().getRadiativeTransferEquation();
+ setCalculationStatus(rte.compute(getCurrentSolution()));
+ }
+ coupling.getRadiativeTransferEquation().getFluxes().store();
+ }
+
+ @Override
+ public Set listedKeywords() {
+ var set = super.listedKeywords();
+ set.add(NONLINEAR_PRECISION);
+ return set;
+ }
+
+ @Override
+ public boolean normalOperation() {
+ return super.normalOperation() && (getCalculationStatus() == RTECalculationStatus.NORMAL);
+ }
+
+ public final RTECalculationStatus getCalculationStatus() {
+ return calculationStatus;
+ }
+
+ public final void setCalculationStatus(RTECalculationStatus calculationStatus) throws SolverException {
+ this.calculationStatus = calculationStatus;
+ if (calculationStatus != RTECalculationStatus.NORMAL) {
+ throw new SolverException(calculationStatus.toString(),
+ RTE_SOLVER_ERROR);
+ }
+ }
+
+ public final RadiativeTransferCoupling getCoupling() {
+ return coupling;
+ }
+
+ public final void setCoupling(RadiativeTransferCoupling coupling) {
+ this.coupling = coupling;
+ this.coupling.setParent(this);
+ }
+
+ public final boolean isAutoUpdateFluxes() {
+ return this.autoUpdateFluxes;
+ }
+
+ public final void setAutoUpdateFluxes(boolean auto) {
+ this.autoUpdateFluxes = auto;
+ }
+
+ @Override
+ public Class[] domain() {
+ return new Class[]{ParticipatingMedium.class};
+ }
+
+}
diff --git a/src/main/java/pulse/problem/schemes/DifferenceScheme.java b/src/main/java/pulse/problem/schemes/DifferenceScheme.java
index 5bb773df..64de3f8c 100644
--- a/src/main/java/pulse/problem/schemes/DifferenceScheme.java
+++ b/src/main/java/pulse/problem/schemes/DifferenceScheme.java
@@ -1,18 +1,18 @@
package pulse.problem.schemes;
+import java.util.Objects;
import static pulse.properties.NumericProperties.def;
import static pulse.properties.NumericProperties.derive;
import static pulse.properties.NumericProperty.requireType;
import static pulse.properties.NumericPropertyKeyword.TIME_LIMIT;
-import java.util.ArrayList;
-import java.util.List;
+import java.util.Set;
import pulse.problem.laser.DiscretePulse;
+import pulse.problem.schemes.solvers.SolverException;
import pulse.problem.statements.Problem;
import pulse.properties.NumericProperty;
import pulse.properties.NumericPropertyKeyword;
-import pulse.properties.Property;
import pulse.util.PropertyHolder;
import pulse.util.Reflexive;
@@ -23,288 +23,304 @@
* partitioning, adjusted to ensure a stable or conditionally-stable behaviour
* of the solution. The {@code Grid} is also used to define a
* {@code DiscretePulse} function.
- *
+ *
* @see pulse.problem.schemes.Grid
* @see pulse.problem.laser.DiscretePulse
*/
-
public abstract class DifferenceScheme extends PropertyHolder implements Reflexive {
- private DiscretePulse discretePulse;
- private Grid grid;
-
- private double timeLimit;
- private int timeInterval;
-
- private static boolean hideDetailedAdjustment = true;
-
- private final static double EPS = 1e-7; // a small value ensuring numeric stability
-
- /**
- * A constructor which merely sets the time limit to its default value.
- */
-
- protected DifferenceScheme() {
- setTimeLimit(def(TIME_LIMIT));
- }
-
- /**
- * A constructor for setting the time limit to a pre-set value.
- *
- * @param timeLimit the calculation time limit
- */
-
- protected DifferenceScheme(NumericProperty timeLimit) {
- setTimeLimit(timeLimit);
- }
-
- /**
- * Used to get a class of problems on which this difference scheme is
- * applicable.
- *
- * @return a subclass of the {@code Problem} class which can be used as input
- * for this difference scheme.
- */
-
- public abstract Class domain();
-
- /**
- * Creates a {@code DifferenceScheme}, which is an exact copy of this object.
- *
- * @return an exact copy of this {@code DifferenceScheme}.
- */
-
- public abstract DifferenceScheme copy();
-
- /**
- * Copies the {@code Grid} and {@code timeLimit} from {@code df}.
- *
- * @param df the DifferenceScheme to copy from
- */
-
- public void copyFrom(DifferenceScheme df) {
- this.grid = df.getGrid().copy();
- discretePulse = null;
- timeLimit = df.timeLimit;
- }
-
- /**
- *
- * Contains preparatory steps to ensure smooth running of the solver. This
- * includes creating a {@code DiscretePulse} object and calculating the
- * {@code timeInterval}. The latter determines the real-time calculation of a
- * {@code HeatingCurve} based on the numerical solution of {@code problem}; it
- * thus takes into account the difference between the scheme timestep and the
- * {@code HeatingCurve} point spacing. All subclasses of
- * {@code DifferenceScheme} should override and explicitly call this superclass
- * method where appropriate.
- *
- *
- * @param problem the heat problem to be solved
- */
-
- protected void prepare(Problem problem) {
- discretePulse = problem.discretePulseOn(grid);
- grid.adjustTo(discretePulse);
-
- var hc = problem.getHeatingCurve();
- hc.clear();
- }
-
- public void runTimeSequence(Problem problem) {
- runTimeSequence(problem, 0, timeLimit);
- var curve = problem.getHeatingCurve();
- final double maxTemp = (double) problem.getProperties().getMaximumTemperature().getValue();
- curve.scale(maxTemp / curve.apparentMaximum() );
- }
-
- public void runTimeSequence(Problem problem, final double offset, final double endTime) {
- final var grid = getGrid();
-
- var curve = problem.getHeatingCurve();
-
- int adjustedNumPoints = (int)curve.getNumPoints().getValue();
-
- final double startTime = (double)curve.getTimeShift().getValue();
- final double timeSegment = (endTime - startTime - offset) / problem.getProperties().timeFactor();
- final double tau = grid.getTimeStep();
-
- for (double dt = 0, factor = 1.0; dt < tau; adjustedNumPoints *= factor) {
- dt = timeSegment / (adjustedNumPoints - 1);
- factor = dt / tau;
- timeInterval = (int) factor;
- }
-
- final double wFactor = timeInterval * tau * problem.getProperties().timeFactor();
-
- // First point (index = 0) is always (0.0, 0.0)
-
- /*
- * The outer cycle iterates over the number of points of the HeatingCurve
- */
-
- double nextTime = offset + wFactor;
- curve.addPoint(0.0, 0.0);
-
- for (int w = 1; nextTime < 1.01*endTime; nextTime = offset + (++w)*wFactor) {
-
- /*
- * Two adjacent points of the heating curves are separated by timeInterval on
- * the time grid. Thus, to calculate the next point on the heating curve,
- * timeInterval/tau time steps have to be made first.
- */
-
- timeSegment((w - 1) * timeInterval + 1, w * timeInterval + 1);
- curve.addPoint(nextTime, signal());
-
- }
-
- }
-
- private void timeSegment(final int m1, final int m2) {
- for (int m = m1; m < m2 && normalOperation(); m++) {
- timeStep(m);
- finaliseStep();
- }
- }
-
- public double pulse(final int m) {
- return getDiscretePulse().laserPowerAt((m - EPS) * getGrid().getTimeStep());
- }
-
- public abstract double signal();
-
- public abstract void timeStep(final int m);
-
- public abstract void finaliseStep();
-
- public boolean normalOperation() {
- return true;
- }
-
- /**
- * The superclass only lists the {@code TIME_LIMIT} property.
- */
-
- @Override
- public List listedTypes() {
- List list = new ArrayList<>();
- list.add(def(TIME_LIMIT));
- return list;
- }
-
- @Override
- public String toString() {
- return this.getClass().getSimpleName();
- }
-
- /**
- * Gets the discrete representation of {@code Pulse} on the {@code Grid}.
- *
- * @return the discrete pulse
- * @see pulse.problem.statements.Pulse
- */
-
- public DiscretePulse getDiscretePulse() {
- return discretePulse;
- }
-
- /**
- * Gets the {@code Grid} object defining partioning used in this
- * {@code DifferenceScheme}
- *
- * @return the grid
- */
-
- public Grid getGrid() {
- return grid;
- }
-
- /**
- * Sets the grid and adopts it as its child.
- *
- * @param grid the grid
- */
-
- public void setGrid(Grid grid) {
- this.grid = grid;
- this.grid.setParent(this);
- }
-
- /**
- * The time interval is the number of discrete timesteps that will be discarded
- * when storing the resulting solution into a {@code HeatingCurve} object, thus
- * ensuring that only a limited set of points is stored.
- *
- * @return the time interval
- */
-
- public int getTimeInterval() {
- return timeInterval;
- }
-
- /**
- * Sets the time interval to the argument of this method.
- *
- * @param timeInterval a positive integer.
- */
-
- public void setTimeInterval(int timeInterval) {
- this.timeInterval = timeInterval;
- }
-
- /**
- * If true, Lets the UI know that the user only wants to have the most important
- * properties displayed. Otherwise this will signal all properties need to be
- * displayed.
- */
-
- @Override
- public boolean areDetailsHidden() {
- return hideDetailedAdjustment;
- }
-
- /**
- * Changes the policy of displaying a detailed information about this scheme.
- *
- * @param b a boolean.
- */
-
- public static void setDetailsHidden(boolean b) {
- hideDetailedAdjustment = b;
- }
-
- /**
- * The time limit (in whatever units this {@code DifferenceScheme} uses to
- * process the solution), which serves as the ultimate breakpoint for the
- * calculations.
- *
- * @return the {@code NumericProperty} with the type {@code TIME_LIMIT}
- * @see pulse.properties.NumericPropertyKeyword
- */
-
- public NumericProperty getTimeLimit() {
- return derive(TIME_LIMIT, timeLimit);
- }
-
- /**
- * Sets the time limit (in units defined by the corresponding
- * {@code NumericProperty}), which serves as the breakpoint for the
- * calculations.
- *
- * @param timeLimit the {@code NumericProperty} with the type {@code TIME_LIMIT}
- * @see pulse.properties.NumericPropertyKeyword
- */
-
- public void setTimeLimit(NumericProperty timeLimit) {
- requireType(timeLimit, TIME_LIMIT);
- this.timeLimit = (double) timeLimit.getValue();
- }
-
- @Override
- public void set(NumericPropertyKeyword type, NumericProperty property) {
- if (type == TIME_LIMIT)
- setTimeLimit(property);
- }
-
-}
\ No newline at end of file
+ private transient DiscretePulse discretePulse;
+ private Grid grid;
+
+ private double timeLimit;
+ private double pls;
+ private int timeInterval;
+
+ private static boolean hideDetailedAdjustment = true;
+
+ private final static double EPS = 1e-7; // a small value ensuring numeric stability
+
+ /**
+ * A constructor which merely sets the time limit to its default value.
+ */
+ protected DifferenceScheme() {
+ setTimeLimit(def(TIME_LIMIT));
+ }
+
+ /**
+ * A constructor for setting the time limit to a pre-set value.
+ *
+ * @param timeLimit the calculation time limit
+ */
+ protected DifferenceScheme(NumericProperty timeLimit) {
+ setTimeLimit(timeLimit);
+ }
+
+ public void initFrom(DifferenceScheme another) {
+ this.grid = grid.copy();
+ this.timeLimit = another.timeLimit;
+ this.timeInterval = another.timeInterval;
+ }
+
+ /**
+ * Copies the {@code Grid} and {@code timeLimit} from {@code df}.
+ *
+ * @param df the DifferenceScheme to copy from
+ */
+ public void copyFrom(DifferenceScheme df) {
+ this.grid = df.getGrid().copy();
+ discretePulse = null;
+ timeLimit = df.timeLimit;
+ }
+
+ /**
+ *
+ * Contains preparatory steps to ensure smooth running of the solver.This
+ * includes creating a {@code DiscretePulse}object and adjusting the grid of
+ * this scheme to match the {@code DiscretePulse}created for this
+ * {@code problem} Finally, a heating curve is cleared from the previously
+ * calculated values.
+ *
+ * All subclasses of {@code DifferenceScheme} should override and explicitly
+ * call this superclass method where appropriate.
+ *
+ *
+ * @param problem the heat problem to be solved
+ * @throws pulse.problem.schemes.solvers.SolverException
+ * @see pulse.problem.schemes.Grid.adjustTo()
+ */
+ protected void prepare(Problem problem) throws SolverException {
+ if (discretePulse == null) {
+ discretePulse = problem.discretePulseOn(grid);
+ }
+ discretePulse.init();
+ clearArrays();
+ }
+
+ public void runTimeSequence(Problem problem) throws SolverException {
+ runTimeSequence(problem, 0, timeLimit);
+ }
+
+ public void scaleSolution(Problem problem) {
+ var curve = problem.getHeatingCurve();
+ final double maxTemp = (double) problem.getProperties().getMaximumTemperature().getValue();
+ //curve.scale(maxTemp / curve.apparentMaximum());
+ curve.scale(maxTemp);
+ }
+
+ public void runTimeSequence(Problem problem, final double offset, final double endTime) throws SolverException {
+ var curve = problem.getHeatingCurve();
+ curve.clear();
+
+ int numPoints = (int) curve.getNumPoints().getValue();
+
+ final double startTime = (double) curve.getTimeShift().getValue();
+ final double timeSegment = (endTime - startTime - offset) / problem.getProperties().characteristicTime();
+
+ double tau = grid.getTimeStep();
+ final double dt = timeSegment / (numPoints - 1);
+ timeInterval = Math.max((int) (dt / tau), 1);
+
+ double wFactor = timeInterval * tau * problem.getProperties().characteristicTime();
+
+ // First point (index = 0) is always (0.0, 0.0)
+ curve.addPoint(0.0, 0.0);
+
+ double nextTime;
+ int previous;
+
+ /*
+ * The outer cycle iterates over the number of points of the HeatingCurve
+ */
+ for (previous = 1, nextTime = offset; nextTime < endTime || !curve.isFull();
+ previous += timeInterval) {
+
+ /*
+ * Two adjacent points of the heating curves are separated by timeInterval on
+ * the time grid. Thus, to calculate the next point on the heating curve,
+ * timeInterval/tau time steps have to be made first.
+ */
+ timeSegment(previous, previous + timeInterval);
+ nextTime += wFactor;
+ curve.addPoint(nextTime, signal());
+ }
+
+ curve.copyToLastCalculation();
+ scaleSolution(problem);
+ }
+
+ private void timeSegment(final int m1, final int m2) throws SolverException {
+ for (int m = m1; m < m2 && normalOperation(); m++) {
+ prepareStep(m); //prepare
+ timeStep(m); //calculate
+ finaliseStep(); //finalise
+ }
+ }
+
+ public double pulse(final int m) {
+ return getDiscretePulse().laserPowerAt((m - EPS) * getGrid().getTimeStep());
+ }
+
+ /**
+ * Do preparatory calculations that depend only on the time variable, e.g.,
+ * calculate the pulse power.
+ *
+ * @param m the time step number
+ */
+ public void prepareStep(int m) {
+ pls = pulse(m);
+ }
+
+ public boolean normalOperation() {
+ return true;
+ }
+
+ /**
+ * The superclass only lists the {@code TIME_LIMIT} property.
+ */
+ @Override
+ public Set listedKeywords() {
+ var set = super.listedKeywords();
+ set.add(TIME_LIMIT);
+ return set;
+ }
+
+ @Override
+ public String toString() {
+ return this.getClass().getSimpleName();
+ }
+
+ /**
+ * Gets the discrete representation of {@code Pulse} on the {@code Grid}.
+ *
+ * @return the discrete pulse
+ * @see pulse.problem.statements.Pulse
+ */
+ public final DiscretePulse getDiscretePulse() {
+ return discretePulse;
+ }
+
+ /**
+ * Gets the {@code Grid} object defining partioning used in this
+ * {@code DifferenceScheme}
+ *
+ * @return the grid
+ */
+ public final Grid getGrid() {
+ return grid;
+ }
+
+ /**
+ * Sets the grid and adopts it as its child.
+ *
+ * @param grid the grid
+ */
+ public final void setGrid(Grid grid) {
+ this.grid = grid;
+ this.grid.setParent(this);
+ }
+
+ /**
+ * The time interval is the number of discrete timesteps that will be
+ * discarded when storing the resulting solution into a {@code HeatingCurve}
+ * object, thus ensuring that only a limited set of points is stored.
+ *
+ * @return the time interval
+ */
+ public final int getTimeInterval() {
+ return timeInterval;
+ }
+
+ /**
+ * Sets the time interval to the argument of this method.
+ *
+ * @param timeInterval a positive integer.
+ */
+ public final void setTimeInterval(int timeInterval) {
+ this.timeInterval = timeInterval;
+ }
+
+ /**
+ * If true, Lets the UI know that the user only wants to have the most
+ * important properties displayed. Otherwise this will signal all properties
+ * need to be displayed.
+ */
+ @Override
+ public final boolean areDetailsHidden() {
+ return hideDetailedAdjustment;
+ }
+
+ /**
+ * Changes the policy of displaying a detailed information about this
+ * scheme.
+ *
+ * @param b a boolean.
+ */
+ public final static void setDetailsHidden(boolean b) {
+ hideDetailedAdjustment = b;
+ }
+
+ /**
+ * The time limit (in whatever units this {@code DifferenceScheme} uses to
+ * process the solution), which serves as the ultimate breakpoint for the
+ * calculations.
+ *
+ * @return the {@code NumericProperty} with the type {@code TIME_LIMIT}
+ * @see pulse.properties.NumericPropertyKeyword
+ */
+ public final NumericProperty getTimeLimit() {
+ return derive(TIME_LIMIT, timeLimit);
+ }
+
+ public double getCurrentPulseValue() {
+ return pls;
+ }
+
+ /**
+ * Sets the time limit (in units defined by the corresponding
+ * {@code NumericProperty}), which serves as the breakpoint for the
+ * calculations.
+ *
+ * @param timeLimit the {@code NumericProperty} with the type
+ * {@code TIME_LIMIT}
+ * @see pulse.properties.NumericPropertyKeyword
+ */
+ public final void setTimeLimit(NumericProperty timeLimit) {
+ requireType(timeLimit, TIME_LIMIT);
+ this.timeLimit = (double) timeLimit.getValue();
+ firePropertyChanged(this, timeLimit);
+ }
+
+ @Override
+ public void set(NumericPropertyKeyword type, NumericProperty property) {
+ if (type == TIME_LIMIT) {
+ setTimeLimit(property);
+ }
+ }
+
+ public abstract double signal();
+
+ public abstract void clearArrays();
+
+ public abstract void timeStep(int m) throws SolverException;
+
+ public abstract void finaliseStep() throws SolverException;
+
+ /**
+ * Retrieves all problem statements that can be solved with this
+ * implementation of the difference scheme.
+ *
+ * @return an array containing subclasses of the {@code Problem} class which
+ * can be used as input for this difference scheme.
+ */
+ public abstract Class[] domain();
+
+ /**
+ * Creates a {@code DifferenceScheme}, which is an exact copy of this
+ * object.
+ *
+ * @return an exact copy of this {@code DifferenceScheme}.
+ */
+ public abstract DifferenceScheme copy();
+
+}
diff --git a/src/main/java/pulse/problem/schemes/DistributedDetection.java b/src/main/java/pulse/problem/schemes/DistributedDetection.java
index a32f06d5..3d179465 100644
--- a/src/main/java/pulse/problem/schemes/DistributedDetection.java
+++ b/src/main/java/pulse/problem/schemes/DistributedDetection.java
@@ -1,38 +1,41 @@
package pulse.problem.schemes;
+import java.io.Serializable;
import java.util.stream.IntStream;
-import pulse.problem.statements.penetration.AbsorptionModel;
-import pulse.problem.statements.penetration.SpectralRange;
+import pulse.problem.statements.model.AbsorptionModel;
+import pulse.problem.statements.model.SpectralRange;
/**
- * An interface providing the ability to calculate the integral signal
- * out from a finite-depth material layer. The depth is governed by
- * the current {@code AbsorptionModel}.
+ * An interface providing the ability to calculate the integral signal out from
+ * a finite-depth material layer. The depth is governed by the current
+ * {@code AbsorptionModel}.
*
*/
-
-public class DistributedDetection {
-
- /**
- * Calculates the effective signal registered by the detector, which takes into account
- * a distributed emission pattern. The emissivity is assumed equal to the average absorptivity
- * in the thermal region of the spectrum, as per the Kirchhoff's law.
- * @param absorption the absorption model
- * @param V the current time-temperature profile
- * @return the effective detector signal (arbitrary units)
- */
-
- public static double evaluateSignal(final AbsorptionModel absorption, final Grid grid, final double[] V) {
- final double hx = grid.getXStep();
- final int N = grid.getGridDensityValue();
-
- double signal = IntStream.range(0, N)
- .mapToDouble(i -> V[N - i] * absorption.absorption(SpectralRange.THERMAL, i * hx)
- + V[N - 1 - i] * absorption.absorption(SpectralRange.THERMAL, (i + 1) * hx))
- .reduce((a, b) -> a + b).getAsDouble();
-
- return signal * 0.5 * hx;
- }
-
-}
\ No newline at end of file
+public class DistributedDetection implements Serializable {
+
+ private static final long serialVersionUID = 3587781877001360511L;
+
+ /**
+ * Calculates the effective signal registered by the detector, which takes
+ * into account a distributed emission pattern. The emissivity is assumed
+ * equal to the average absorptivity in the thermal region of the spectrum,
+ * as per the Kirchhoff's law.
+ *
+ * @param absorption the absorption model
+ * @param V the current time-temperature profile
+ * @return the effective detector signal (arbitrary units)
+ */
+ public static double evaluateSignal(final AbsorptionModel absorption, final Grid grid, final double[] V) {
+ final double hx = grid.getXStep();
+ final int N = grid.getGridDensityValue();
+
+ double signal = IntStream.range(0, N)
+ .mapToDouble(i -> V[N - i] * absorption.absorption(SpectralRange.THERMAL, i * hx)
+ + V[N - 1 - i] * absorption.absorption(SpectralRange.THERMAL, (i + 1) * hx))
+ .reduce((a, b) -> a + b).getAsDouble();
+
+ return signal * 0.5 * hx;
+ }
+
+}
diff --git a/src/main/java/pulse/problem/schemes/ExplicitScheme.java b/src/main/java/pulse/problem/schemes/ExplicitScheme.java
index e5fb12a4..e4b858cd 100644
--- a/src/main/java/pulse/problem/schemes/ExplicitScheme.java
+++ b/src/main/java/pulse/problem/schemes/ExplicitScheme.java
@@ -12,85 +12,85 @@
* This class provides the necessary framework to enable a simple explicit
* finite-difference scheme (also called the forward-time centred space scheme)
* for solving the one-dimensional heat conduction problem.
- *
+ *
* @see pulse.problem.statements.ClassicalProblem
* @see pulse.problem.statements.NonlinearProblem
*
*/
-
public abstract class ExplicitScheme extends OneDimensionalScheme {
- /**
- * Constructs a default explicit scheme using the default values of
- * {@code GRID_DENSITY} and {@code TAU_FACTOR}.
- */
-
- public ExplicitScheme() {
- this(derive(GRID_DENSITY, 80), derive(TAU_FACTOR, 0.5));
- }
+ /**
+ *
+ */
+ private static final long serialVersionUID = -3025913683505686334L;
- /**
- * Constructs an explicit scheme on a one-dimensional grid that is specified by
- * the values {@code N} and {@code timeFactor}.
- *
- * @see pulse.problem.schemes.DifferenceScheme
- * @param N the {@code NumericProperty} with the type
- * {@code GRID_DENSITY}
- * @param timeFactor the {@code NumericProperty} with the type
- * {@code TAU_FACTOR}
- */
+ /**
+ * Constructs a default explicit scheme using the default values of
+ * {@code GRID_DENSITY} and {@code TAU_FACTOR}.
+ */
+ public ExplicitScheme() {
+ this(derive(GRID_DENSITY, 80), derive(TAU_FACTOR, 0.5));
+ }
- public ExplicitScheme(NumericProperty N, NumericProperty timeFactor) {
- super();
- setGrid(new Grid(N, timeFactor));
- }
+ /**
+ * Constructs an explicit scheme on a one-dimensional grid that is specified
+ * by the values {@code N} and {@code timeFactor}.
+ *
+ * @see pulse.problem.schemes.DifferenceScheme
+ * @param N the {@code NumericProperty} with the type {@code GRID_DENSITY}
+ * @param timeFactor the {@code NumericProperty} with the type
+ * {@code TAU_FACTOR}
+ */
+ public ExplicitScheme(NumericProperty N, NumericProperty timeFactor) {
+ super();
+ setGrid(new Grid(N, timeFactor));
+ }
- /**
- *
- * Constructs an explicit scheme on a one-dimensional grid that is specified by
- * the values {@code N} and {@code timeFactor}. Sets the time limit of this
- * scheme to {@code timeLimit}
- *
- * @param N the {@code NumericProperty} with the type
- * {@code GRID_DENSITY}
- * @param timeFactor the {@code NumericProperty} with the type
- * {@code TAU_FACTOR}
- * @param timeLimit the {@code NumericProperty} with the type
- * {@code TIME_LIMIT}
- * @see pulse.problem.schemes.DifferenceScheme
- */
+ /**
+ *
+ * Constructs an explicit scheme on a one-dimensional grid that is specified
+ * by the values {@code N} and {@code timeFactor}. Sets the time limit of
+ * this scheme to {@code timeLimit}
+ *
+ * @param N the {@code NumericProperty} with the type {@code GRID_DENSITY}
+ * @param timeFactor the {@code NumericProperty} with the type
+ * {@code TAU_FACTOR}
+ * @param timeLimit the {@code NumericProperty} with the type
+ * {@code TIME_LIMIT}
+ * @see pulse.problem.schemes.DifferenceScheme
+ */
+ public ExplicitScheme(NumericProperty N, NumericProperty timeFactor, NumericProperty timeLimit) {
+ super(timeLimit);
+ setGrid(new Grid(N, timeFactor));
+ }
- public ExplicitScheme(NumericProperty N, NumericProperty timeFactor, NumericProperty timeLimit) {
- super(timeLimit);
- setGrid(new Grid(N, timeFactor));
- }
-
- /**
- * Uses the explicit finite-difference representation of the heat equation to calculate the grid-function everywhere
- * except for the boundaries. This will update the current solution using the solution from previous time step.
- */
-
- public void explicitSolution() {
- var grid = getGrid();
- var U = getPreviousSolution();
- final double TAU_HH = grid.getTimeStep()/(fastPowLoop(grid.getXStep(), 2));
- for (int i = 1, N = grid.getGridDensityValue(); i < N; i++)
- setSolutionAt(i, U[i] + TAU_HH * (U[i + 1] - 2. * U[i] + U[i - 1]) + phi(i) );
- }
-
- public double phi(final int i) {
- return 0;
- }
+ /**
+ * Uses the explicit finite-difference representation of the heat equation
+ * to calculate the grid-function everywhere except for the boundaries. This
+ * will update the current solution using the solution from previous time
+ * step.
+ */
+ public void explicitSolution() {
+ var grid = getGrid();
+ var U = getPreviousSolution();
+ final double TAU_HH = grid.getTimeStep() / (fastPowLoop(grid.getXStep(), 2));
+ for (int i = 1, N = grid.getGridDensityValue(); i < N; i++) {
+ setSolutionAt(i, U[i] + TAU_HH * (U[i + 1] - 2. * U[i] + U[i - 1]) + phi(i));
+ }
+ }
- /**
- * Prints out the description of this problem type.
- *
- * @return a verbose description of the problem.
- */
+ public double phi(final int i) {
+ return 0;
+ }
- @Override
- public String toString() {
- return getString("ExplicitScheme.4");
- }
+ /**
+ * Prints out the description of this problem type.
+ *
+ * @return a verbose description of the problem.
+ */
+ @Override
+ public String toString() {
+ return getString("ExplicitScheme.4");
+ }
-}
\ No newline at end of file
+}
diff --git a/src/main/java/pulse/problem/schemes/FixedPointIterations.java b/src/main/java/pulse/problem/schemes/FixedPointIterations.java
index a093e67f..dbf5992f 100644
--- a/src/main/java/pulse/problem/schemes/FixedPointIterations.java
+++ b/src/main/java/pulse/problem/schemes/FixedPointIterations.java
@@ -2,26 +2,70 @@
import static java.lang.Math.abs;
-public interface FixedPointIterations {
+import java.io.Serializable;
+import java.util.Arrays;
+import pulse.problem.schemes.solvers.SolverException;
+import static pulse.problem.schemes.solvers.SolverException.SolverExceptionType.FINITE_DIFFERENCE_ERROR;
- public default void doIterations(double[] V, final double error, final int m) {
+/**
+ * @see Wiki
+ * page
+ *
+ */
+public interface FixedPointIterations extends Serializable {
- final int N = V.length - 1;
+ /**
+ * Performs iterations until the convergence criterion is satisfied.The
+ * latter consists in having a difference two consequent iterations of V
+ * less than the specified error. At the end of each iteration, calls
+ * {@code finaliseIteration()}.
+ *
+ * @param V the calculation array
+ * @param error used in the convergence criterion
+ * @param m time step
+ * @throws pulse.problem.schemes.solvers.SolverException if the calculation
+ * failed
+ * @see finaliseIteration()
+ * @see iteration()
+ */
+ public default void doIterations(double[] V, final double error, final int m) throws SolverException {
- for (double V_0 = error + 1, V_N = error + 1; abs(V[0] - V_0) > error
- || abs(V[N] - V_N) > error; finaliseIteration(V)) {
+ final int N = V.length - 1;
- V_N = V[N];
- V_0 = V[0];
- iteration(m);
+ for (double V_0 = error + 1, V_N = error + 1;
+ abs(V[0] - V_0) / abs(V[0] + V_0 + 1e-16) > error
+ || abs(V[N] - V_N) / abs(V[N] + V_N + 1e-16) > error; finaliseIteration(V)) {
- }
- }
+ V_N = V[N];
+ V_0 = V[0];
+ iteration(m);
- public void iteration(final int m);
+ }
+ }
- public default void finaliseIteration(double[] V) {
- // do nothing
- }
+ /**
+ * Performs an iteration at time {@code m}
+ *
+ * @param m time step
+ * @throws pulse.problem.schemes.solvers.SolverException if the calculation
+ * failed
+ */
+ public void iteration(final int m) throws SolverException;
-}
\ No newline at end of file
+ /**
+ * Finalises the current iteration.By default, does nothing.
+ *
+ * @param V the current iteration
+ * @throws pulse.problem.schemes.solvers.SolverException if the calculation
+ * failed
+ */
+ public default void finaliseIteration(double[] V) throws SolverException {
+ final double threshold = 1E6;
+ double sum = Arrays.stream(V).sum();
+ if (sum > threshold || !Double.isFinite(sum)) {
+ throw new SolverException("Invalid solution values in V array",
+ FINITE_DIFFERENCE_ERROR);
+ }
+ }
+
+}
diff --git a/src/main/java/pulse/problem/schemes/Grid.java b/src/main/java/pulse/problem/schemes/Grid.java
index 20537fb8..d8a0b331 100644
--- a/src/main/java/pulse/problem/schemes/Grid.java
+++ b/src/main/java/pulse/problem/schemes/Grid.java
@@ -2,20 +2,15 @@
import static java.lang.Math.pow;
import static java.lang.Math.rint;
-import static java.lang.String.format;
-import static pulse.properties.NumericProperties.def;
import static pulse.properties.NumericProperties.derive;
import static pulse.properties.NumericProperty.requireType;
import static pulse.properties.NumericPropertyKeyword.GRID_DENSITY;
import static pulse.properties.NumericPropertyKeyword.TAU_FACTOR;
-import java.util.ArrayList;
-import java.util.List;
+import java.util.Set;
-import pulse.problem.laser.DiscretePulse;
import pulse.properties.NumericProperty;
import pulse.properties.NumericPropertyKeyword;
-import pulse.properties.Property;
import pulse.util.PropertyHolder;
/**
@@ -28,218 +23,196 @@
*
*
*/
-
public class Grid extends PropertyHolder {
- private double hx;
- private double tau;
- private double tauFactor;
- private int N;
-
- /**
- * Creates a {@code Grid} object with the specified {@code gridDensity} and
- * {@code timeFactor}.
- *
- * @param gridDensity a {@code NumericProperty} of the type {@code GRID_DENSITY}
- * @param timeFactor a {@code NumericProperty} of the type {@code TIME_FACTOR}
- * @see pulse.properties.NumericPropertyKeyword
- */
-
- public Grid(NumericProperty gridDensity, NumericProperty timeFactor) {
- setGridDensity(gridDensity);
- setTimeFactor(timeFactor);
- }
-
- protected Grid() {
- // intentionally blank
- }
-
- /**
- * Creates a new {@code Grid} object with exactly the same parameters as this
- * one.
- *
- * @return a new {@code Grid} object replicating this {@code Grid}
- */
-
- public Grid copy() {
- return new Grid(getGridDensity(), getTimeFactor());
- }
-
- /**
- * Optimises the {@code Grid} parameters.
- *
- * This can change the {@code tauFactor} and {@code tau} variables in the
- * {@code Grid} object if {@code discretePulseWidth < grid.tau}.
- *
- *
- * @param pulse the discrete pulse representation
- */
-
- public void adjustTo(DiscretePulse pulse) {
- final double ADJUSTMENT_FACTOR = 0.75;
- for (final double factor = 0.95; factor * tau > pulse.getDiscreteWidth(); pulse.recalculate()) {
- tauFactor *= ADJUSTMENT_FACTOR;
- tau = tauFactor * pow(hx, 2);
- }
- }
-
- /**
- * The listed properties include {@code GRID_DENSITY} and {@code TAU_FACTOR}.
- */
-
- @Override
- public List listedTypes() {
- List list = new ArrayList<>(2);
- list.add(def(GRID_DENSITY));
- list.add(def(TAU_FACTOR));
- return list;
- }
-
- @Override
- public void set(NumericPropertyKeyword type, NumericProperty property) {
- switch (type) {
- case TAU_FACTOR:
- setTimeFactor(property);
- break;
- case GRID_DENSITY:
- setGridDensity(property);
- break;
- default:
- break;
- }
- }
-
- /**
- * Retrieves the value of the coordinate step
- * used in finite-difference calculation.
- *
- * @return a double, representing the {@code hx} value.
- */
-
- public double getXStep() {
- return hx;
- }
-
- /**
- * Sets the value of the coordinate step.
- *
- * @param hx a double, representing the new {@code hx} value.
- */
-
- public void setXStep(double hx) {
- this.hx = hx;
- }
-
- /**
- * Retrieves the value of the τ time step used in finite-difference
- * calculation.
- *
- * @return a double, representing the {@code tau} value.
- */
-
- public double getTimeStep() {
- return tau;
- }
-
- protected void setTimeStep(double tau) {
- this.tau = tau;
- }
-
- /**
- * Retrieves the value of the τ-factor, or the time factor, used in
- * finite-difference calculation. This factor determines the proportionally
- * coefficient between τ and .
- *
- * @return a NumericProperty of the {@code TAU_FACTOR} type, representing the
- * {@code tauFactor} value.
- */
-
- public NumericProperty getTimeFactor() {
- return derive(TAU_FACTOR, tauFactor);
- }
-
- /**
- * Retrieves the value of the {@code gridDensity} used to calculate the
- * {@code hx} and {@code tau}.
- *
- * @return a NumericProperty of the {@code GRID_DENSITY} type, representing the
- * {@code gridDensity} value.
- */
-
- public NumericProperty getGridDensity() {
- return derive(GRID_DENSITY, N);
- }
-
- protected int getGridDensityValue() {
- return N;
- }
-
- protected void setGridDensityValue(int N) {
- this.N = N;
- }
-
- /**
- * Sets the value of the {@code gridDensity}. Automatically recalculates the
- * {@code hx} value.
- *
- * @param gridDensity a NumericProperty of the {@code GRID_DENSITY} type
- */
-
- public void setGridDensity(NumericProperty gridDensity) {
- requireType(gridDensity, GRID_DENSITY);
- this.N = (int) gridDensity.getValue();
- hx = 1. / N;
- setTimeFactor(derive(TAU_FACTOR, 1.0));
- }
-
- /**
- * Sets the value of the {@code tauFactor}. Automatically recalculates the
- * {@code tau} value.
- *
- * @param timeFactor a NumericProperty of the {@code TAU_FACTOR} type
- */
-
- public void setTimeFactor(NumericProperty timeFactor) {
- requireType(timeFactor, TAU_FACTOR);
- this.tauFactor = (double) timeFactor.getValue();
- setTimeStep(tauFactor * pow(hx, 2));
- }
-
- /**
- * The dimensionless time on this {@code Grid}, which is the
- * {@code time/dimensionFactor} rounded up to a factor of the time step
- * {@code tau}.
- *
- * @param time the time
- * @param dimensionFactor a conversion factor with the dimension of time
- * @return a double representing the time on the finite grid
- */
-
- public double gridTime(double time, double dimensionFactor) {
- return rint((time / dimensionFactor) / tau) * tau;
- }
-
- /**
- * The dimensionless axial distance on this {@code Grid}, which is the
- * {@code distance/lengthFactor} rounded up to a factor of the coordinate step
- * {@code hx}.
- *
- * @param distance the distance along the axial direction
- * @param lengthFactor a conversion factor with the dimension of length
- * @return a double representing the axial distance on the finite grid
- */
-
- public double gridAxialDistance(double distance, double lengthFactor) {
- return rint((distance / lengthFactor) / hx) * hx;
- }
-
- @Override
- public String toString() {
- var sb = new StringBuilder();
- sb.append("");
- sb.append(getClass().getSimpleName() + ":
*/
-
public class Grid2D extends Grid {
- private double hy;
-
- protected Grid2D() {
- super();
- }
-
- /**
- * Creates a {@code Grid2D} where the radial and axial spatial steps are equal
- * to the inverse {@code gridDensity}. Otherwise, calls the superclass
- * constructor.
- *
- * @param gridDensity the grid density
- * @param timeFactor the {@code τ