1 files changed, 116 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolator.java
new file mode 100644
index 0000000..89a5dc1
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolator.java
@@ -0,0 +1,116 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+import org.apache.commons.math3.ode.sampling.StepInterpolator;
+
+/**
+ * This class implements a step interpolator for second order
+ * Runge-Kutta integrator.
+ *
+ * <p>This interpolator computes dense output inside the last
+ * step computed. The interpolation equation is consistent with the
+ * integration scheme :
+ * <ul>
+ *   <li>Using reference point at step start:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>) + &theta; h [(1 - &theta;) y'<sub>1</sub> + &theta; y'<sub>2</sub>]
+ *   </li>
+ *   <li>Using reference point at step end:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h) + (1-&theta;) h [&theta; y'<sub>1</sub> - (1+&theta;) y'<sub>2</sub>]
+ *   </li>
+ * </ul>
+ * </p>
+ *
+ * where &theta; belongs to [0 ; 1] and where y'<sub>1</sub> and y'<sub>2</sub> are the two
+ * evaluations of the derivatives already computed during the
+ * step.</p>
+ *
+ * @see MidpointIntegrator
+ * @since 1.2
+ */
+
+class MidpointStepInterpolator
+  extends RungeKuttaStepInterpolator {
+
+  /** Serializable version identifier */
+  private static final long serialVersionUID = 20111120L;
+
+  /** Simple constructor.
+   * This constructor builds an instance that is not usable yet, the
+   * {@link
+   * org.apache.commons.math3.ode.sampling.AbstractStepInterpolator#reinitialize}
+   * method should be called before using the instance in order to
+   * initialize the internal arrays. This constructor is used only
+   * in order to delay the initialization in some cases. The {@link
+   * RungeKuttaIntegrator} class uses the prototyping design pattern
+   * to create the step interpolators by cloning an uninitialized model
+   * and later initializing the copy.
+   */
+  // CHECKSTYLE: stop RedundantModifier
+  // the public modifier here is needed for serialization
+  public MidpointStepInterpolator() {
+  }
+  // CHECKSTYLE: resume RedundantModifier
+
+  /** Copy constructor.
+   * @param interpolator interpolator to copy from. The copy is a deep
+   * copy: its arrays are separated from the original arrays of the
+   * instance
+   */
+  MidpointStepInterpolator(final MidpointStepInterpolator interpolator) {
+    super(interpolator);
+  }
+
+  /** {@inheritDoc} */
+  @Override
+  protected StepInterpolator doCopy() {
+    return new MidpointStepInterpolator(this);
+  }
+
+
+  /** {@inheritDoc} */
+  @Override
+  protected void computeInterpolatedStateAndDerivatives(final double theta,
+                                          final double oneMinusThetaH) {
+
+    final double coeffDot2 = 2 * theta;
+    final double coeffDot1 = 1 - coeffDot2;
+
+    if ((previousState != null) && (theta <= 0.5)) {
+        final double coeff1    = theta * oneMinusThetaH;
+        final double coeff2    = theta * theta * h;
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot1 = yDotK[0][i];
+            final double yDot2 = yDotK[1][i];
+            interpolatedState[i] = previousState[i] + coeff1 * yDot1 + coeff2 * yDot2;
+            interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2;
+        }
+    } else {
+        final double coeff1    = oneMinusThetaH * theta;
+        final double coeff2    = oneMinusThetaH * (1.0 + theta);
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot1 = yDotK[0][i];
+            final double yDot2 = yDotK[1][i];
+            interpolatedState[i] = currentState[i] + coeff1 * yDot1 - coeff2 * yDot2;
+            interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2;
+        }
+    }
+
+  }
+
+}