commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From pste...@apache.org
Subject svn commit: r478458 - in /jakarta/commons/proper/math/trunk/src/mantissa: src/org/spaceroots/mantissa/ src/org/spaceroots/mantissa/geometry/ src/org/spaceroots/mantissa/ode/ tests-src/org/spaceroots/mantissa/geometry/ tests-src/org/spaceroots/mantissa/...
Date Thu, 23 Nov 2006 04:17:14 GMT
Author: psteitz
Date: Wed Nov 22 20:17:13 2006
New Revision: 478458

URL: http://svn.apache.org/viewvc?view=rev&rev=478458
Log:
Collection of patches to initial Mantissa sources:
- Fixed a problem when switching functions triggered derivatives
  discontinuities
- Removed methods and classes that were deprecated in Mantissa
  and do not need to be preserved in commons-math as backward compatibility
  is not a problem for this newly integrated code
- Changed Vector3D and Rotation to immutable classes for ease of use
- Improved some javadoc in class Rotation
JIRA: MATH-161
Submitted (with patches) by Luc Maisonobe

Modified:
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java
    jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java
    jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java Wed Nov 22 20:17:13 2006
@@ -91,7 +91,7 @@
     { "dimensions mismatch: ODE problem has dimension {0},"
     + " state vector has dimension {1}",
       "incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0}),"
-    + " et le vecteur d'\u00e9tat ({1})" },
+    + " et le vecteur d''\u00e9tat ({1})" },
     { "too small integration interval: length = {0}",
       "intervalle d''int\u00e9gration trop petit : {0}" },
 

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java Wed Nov 22 20:17:13 2006
@@ -1,204 +0,0 @@
-// 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.spaceroots.mantissa.geometry;
-
-/**
- * This class implements immutable vectors in a three-dimensional space.
-
- * @version $Id: ImmutableVector3D.java 1705 2006-09-17 19:57:39Z luc $
- * @author L. Maisonobe
-
- */
-
-public class ImmutableVector3D
-  extends Vector3D {
-
-  /** Simple constructor.
-   * Build a vector from its coordinates
-   * @param x abscissa
-   * @param y ordinate
-   * @param z height
-   */
-  public ImmutableVector3D(double x, double y, double z) {
-    super(x, y, z);
-    computeNorm();
-  }
-
-  /** Simple constructor.
-   * Build a vector from its azimuthal coordinates
-   * @param alpha azimuth around Z
-   *              (0 is +X, PI/2 is +Y, PI is -X and 3PI/2 is -Y)
-   * @param delta elevation above (XY) plane, from -PI to +PI
-   */
-  public ImmutableVector3D(double alpha, double delta) {
-    super(alpha, delta);
-    computeNorm();
-  }
-
-  /** Copy constructor.
-   * Build a copy of a vector
-   * @param v vector to copy
-   */
-  public ImmutableVector3D(Vector3D v) {
-    super(v);
-    computeNorm();
-  }
-
-  /** Multiplicative constructor
-   * Build a vector from another one and a scale factor. 
-   * The vector built will be a * u
-   * @param a scale factor
-   * @param u base (unscaled) vector
-   */
-  public ImmutableVector3D(double a, Vector3D u) {
-    super(a, u);
-    computeNorm();
-  }
-
-  /** Linear constructor
-   * Build a vector from two other ones and corresponding scale factors.
-   * The vector built will be a * u +  b * v
-   * @param a first scale factor
-   * @param u first base (unscaled) vector
-   * @param b second scale factor
-   * @param v second base (unscaled) vector
-   */
-  public ImmutableVector3D(double a, Vector3D u, double b, Vector3D v) {
-    super(a, u, b, v);
-    computeNorm();
-  }
-
-  /** Set the abscissa of the vector.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param x new abscissa for the vector
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void setX(double x) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Set the ordinate of the vector.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param y new ordinate for the vector
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void setY(double y) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Set the height of the vector.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param z new height for the vector
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void setZ(double z) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Set all coordinates of the vector.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param x new abscissa for the vector
-   * @param y new ordinate for the vector
-   * @param z new height for the vector
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void setCoordinates(double x, double y, double z) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Compute the norm once and for all. */
-  private void computeNorm() {
-    norm = Math.sqrt(x * x + y * y + z * z);
-  }
-
-  /** Get the norm for the vector.
-   * @return euclidian norm for the vector
-   */
-  public double getNorm() {
-    return norm;
-  }
-
-  /** Add a vector to the instance.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param v vector to add
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void addToSelf(Vector3D v) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Subtract a vector from the instance.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param v vector to subtract
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void subtractFromSelf(Vector3D v) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Normalize the instance.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void normalizeSelf() {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Revert the instance.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void negateSelf() {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Multiply the instance by a scalar
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param a scalar by which the instance should be multiplied
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void multiplySelf(double a) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Reinitialize internal state from the specified array slice data.
-   * This method should not be called for immutable vectors, it always
-   * throws an <code>UnsupportedOperationException</code> exception
-   * @param start start index in the array
-   * @param array array holding the data to extract
-   * @exception UnsupportedOperationException thrown in every case
-   */
-  public void mapStateFromArray(int start, double[] array) {
-    throw new UnsupportedOperationException("vector is immutable");
-  }
-
-  /** Norm of the vector. */
-  private double norm;
-
-  private static final long serialVersionUID = 5377895850033895270L;
-
-}

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java Wed Nov 22 20:17:13 2006
@@ -17,7 +17,6 @@
 
 package org.spaceroots.mantissa.geometry;
 
-import org.spaceroots.mantissa.utilities.ArraySliceMappable;
 import java.io.Serializable;
 
 /**
@@ -25,59 +24,61 @@
 
  * <p>Rotations can be represented by several different mathematical
  * entities (matrices, axe and angle, Cardan or Euler angles,
- * quaternions). This class is an higher level abstraction, more
- * user-oriented and hiding this implementation detail. Well, for the
+ * quaternions). This class presents an higher level abstraction, more
+ * user-oriented and hiding this implementation details. Well, for the
  * curious, we use quaternions for the internal representation. The
  * user can build a rotation from any of these representations, and
  * any of these representations can be retrieved from a
  * <code>Rotation</code> instance (see the various constructors and
  * getters). In addition, a rotation can also be built implicitely
- * from a set of vectors before and after it has been applied. This
- * means that this class can be used to compute transformations from
- * one representation to another one. For example, extracting a set of
- * Cardan angles from a rotation matrix can be done using one single
- * line of code:</p>
+ * from a set of vectors and their image.</p>
+ * <p>This implies that this class can be used to transform one
+ * representation into another one. For example, converting a rotation
+ * matrix into a set of Cardan angles from can be done using the
+ * followong single line of code:</p>
  * <pre>
  * double[] angles = new Rotation(matrix, 1.0e-10).getAngles(RotationOrder.XYZ);
  * </pre>
- * <p>Focus is more oriented on what a rotation <em>do</em>. Once it
- * has been built, and regardless of its representation, a rotation is
- * an <em>operator</em> which basically transforms three dimensional
- * {@link Vector3D vectors} into other three dimensional {@link
- * Vector3D vectors}. Depending on the application, the meaning of
- * these vectors can vary. For example in an attitude simulation tool,
- * you will often consider the vector is fixed and you transform its
- * coordinates in one frame into its coordinates in another frame. In
- * this case, the rotation implicitely defines the relation between
- * the two frames. Another example could be a telescope control application,
- * where the rotation would transform the sighting direction at rest
- * into the desired observing direction. In this case the frame is the
- * same (probably a topocentric one) and the raw and transformed
- * vectors are different. In many case, both approaches will be
- * combined, in our telescope example, we will probably also need to
- * transform the observing direction in the topocentric frame into the
- * observing direction in inertial frame taking into account the
- * observatory location and the earth rotation.</p>
-
- * <p>These examples show that a rotation is what the user wants it to
- * be, so this class does not push the user towards one specific
- * definition and hence does not provide methods like
- * <code>projectVectorIntoDestinationFrame</code> or
- * <code>computeTransformedDirection</code>. It provides simpler and
- * more generic methods: {@link #applyTo(Vector3D) applyTo(Vector3D)}
- * and {@link #applyInverseTo(Vector3D) applyInverseTo(Vector3D)}.</p>
-
- * <p>Since a rotation is basically a vectorial operator, several
- * rotations can be composed together and the composite operation
- * <code>r = r1 o r2</code> (which means that for each vector
- * <code>u</code>, <code>r(u) = r1(r2(u))</code>) is also a
- * rotation. Hence we can consider that in addition to vectors, a
- * rotation can be applied to other rotations (or to itself). With our
- * previous notations, we would say we can apply <code>r1</code> to
- * <code>r2</code> and the result we get is <code>r = r1 o
- * r2</code>. For this purpose, the class provides the methods: {@link
- * #applyTo(Rotation) applyTo(Rotation)} and {@link
- * #applyInverseTo(Rotation) applyInverseTo(Rotation)}.</p>
+ * <p>Focus is oriented on what a rotation <em>do</em> rather than on its
+ * underlying representation. Once it has been built, and regardless of its
+ * internal representation, a rotation is an <em>operator</em> which basically
+ * transforms three dimensional {@link Vector3D vectors} into other three
+ * dimensional {@link Vector3D vectors}. Depending on the application, the
+ * meaning of these vectors may vary and the semantics of the rotation also.</p>
+ * <p>For example in an spacecraft attitude simulation tool, users will often
+ * consider the vectors are fixed (say the Earth direction for example) and the
+ * rotation transforms the coordinates coordinates of this vector in inertial
+ * frame into the coordinates of the same vector in satellite frame. In this
+ * case, the rotation implicitely defines the relation between the two frames.
+ * Another example could be a telescope control application, where the rotation
+ * would transform the sighting direction at rest into the desired observing
+ * direction when the telescope is pointed towards an object of interest. In this
+ * case the rotation transforms the directionf at rest in a topocentric frame
+ * into the sighting direction in the same topocentric frame. In many case, both
+ * approaches will be combined, in our telescope example, we will probably also
+ * need to transform the observing direction in the topocentric frame into the
+ * observing direction in inertial frame taking into account the observatory
+ * location and the Earth rotation.</p>
+
+ * <p>These examples show that a rotation is what the user wants it to be, so this
+ * class does not push the user towards one specific definition and hence does not
+ * provide methods like <code>projectVectorIntoDestinationFrame</code> or
+ * <code>computeTransformedDirection</code>. It provides simpler and more generic
+ * methods: {@link #applyTo(Vector3D) applyTo(Vector3D)} and {@link
+ * #applyInverseTo(Vector3D) applyInverseTo(Vector3D)}.</p>
+
+ * <p>Since a rotation is basically a vectorial operator, several rotations can be
+ * composed together and the composite operation <code>r = r<sub>1</sub> o
+ * r<sub>2</sub></code> (which means that for each vector <code>u</code>,
+ * <code>r(u) = r<sub>1</sub>(r<sub>2</sub>(u))</code>) is also a rotation. Hence
+ * we can consider that in addition to vectors, a rotation can be applied to other
+ * rotations as well (or to itself). With our previous notations, we would say we
+ * can apply <code>r<sub>1</sub></code> to <code>r<sub>2</sub></code> and the result
+ * we get is <code>r = r<sub>1</sub> o r<sub>2</sub></code>. For this purpose, the
+ * class provides the methods: {@link #applyTo(Rotation) applyTo(Rotation)} and
+ * {@link #applyInverseTo(Rotation) applyInverseTo(Rotation)}.</p>
+
+ * <p>Rotations are guaranteed to be immutable objects.</p>
 
  * @version $Id: Rotation.java 1705 2006-09-17 19:57:39Z luc $
  * @author L. Maisonobe
@@ -86,8 +87,7 @@
 
  */
 
-public class Rotation
-  implements ArraySliceMappable, Serializable {
+public class Rotation implements Serializable {
 
   /** Build the identity rotation.
    */
@@ -99,30 +99,11 @@
   }
 
   /** Build a rotation from the quaternion coordinates.
-   * @param q0 scalar part of the quaternion
-   * @param q1 first coordinate of the vectorial part of the quaternion
-   * @param q2 second coordinate of the vectorial part of the quaternion
-   * @param q3 third coordinate of the vectorial part of the quaternion
-   * @deprecated since Mantissa 6.3, this method as been deprecated as it
-   * does not properly handles non-normalized quaternions, it should be
-   * replaced by {@link #Rotation(double, double, double, double, boolean)}
-   */
-  public Rotation(double q0, double q1, double q2, double q3) {
-    this.q0 = q0;
-    this.q1 = q1;
-    this.q2 = q2;
-    this.q3 = q3;
-  }
-
-  /** Build a rotation from the quaternion coordinates.
    * <p>A rotation can be built from a <em>normalized</em> quaternion,
    * i.e. a quaternion for which q<sub>0</sub><sup>2</sup> +
    * q<sub>1</sub><sup>2</sup> + q<sub>2</sub><sup>2</sup> +
    * q<sub>3</sub><sup>2</sup> = 1. If the quaternion is not normalized,
    * the constructor can normalize it in a preprocessing step.</p>
-   * <p>This method replaces the {@link #Rotation(double, double,
-   * double, double) constructor using only 4 doubles} which was deprecated
-   * as of version 6.3 of Mantissa.</p>
    * @param q0 scalar part of the quaternion
    * @param q1 first coordinate of the vectorial part of the quaternion
    * @param q2 second coordinate of the vectorial part of the quaternion
@@ -181,7 +162,7 @@
   /** Build a rotation from a 3X3 matrix.
 
    * <p>Rotation matrices are orthogonal matrices, i.e. unit matrices
-   * (which are matrices for which m.mT = I) with real
+   * (which are matrices for which m.m<sup>T</sup> = I) with real
    * coefficients. The module of the determinant of unit matrices is
    * 1, among the orthogonal 3X3 matrices, only the ones having a
    * positive determinant (+1) are rotation matrices.</p>
@@ -290,14 +271,15 @@
 
   /** Build the rotation that transforms a pair of vector into another pair.
 
-   * <p>Except for possible scale factors, if the instance were
-   * applied to the pair (u1, u2) it will produce the pair (v1,
-   * v2).</p>
-
-   * <p>If the angular separation between u1 and u2 is not the same as
-   * the angular separation between v1 and v2, then a corrected v2'
-   * will be used rather than v2, the corrected vector will be in the
-   * (v1, v2) plane.</p>
+   * <p>Except for possible scale factors, if the instance were applied to
+   * the pair (u<sub>1</sub>, u<sub>2</sub>) it will produce the pair
+   * (v<sub>1</sub>, v<sub>2</sub>).</p>
+
+   * <p>If the angular separation between u<sub>1</sub> and u<sub>2</sub> is
+   * not the same as the angular separation between v<sub>1</sub> and
+   * v<sub>2</sub>, then a corrected v'<sub>2</sub> will be used rather than
+   * v<sub>2</sub>, the corrected vector will be in the (v<sub>1</sub>,
+   * v<sub>2</sub>) plane.</p>
 
    * @param u1 first vector of the origin pair
    * @param u2 second vector of the origin pair
@@ -477,11 +459,11 @@
    * rotations are three successive rotations around the canonical
    * axes X, Y and Z, the first and last rotations beeing around the
    * same axis. There are 6 such sets of rotations (XYX, XZX, YXY,
-   * YZY, ZXZ and ZYZ), the most popular one being ZXZ. Beware that
-   * many people routinely use the term Euler angles even for what
-   * really are Cardan angles (this confusion is especially widespread
-   * in the aerospace business where Roll, Pitch and Yaw angles are
-   * often tagged as Euler angles).</p>
+   * YZY, ZXZ and ZYZ), the most popular one being ZXZ.</p>
+   * <p>Beware that many people routinely use the term Euler angles even
+   * for what really are Cardan angles (this confusion is especially
+   * widespread in the aerospace business where Roll, Pitch and Yaw angles
+   * are often tagged as Euler angles).</p>
 
    * @param order order of rotations to use
    * @param alpha1 angle of the first elementary rotation
@@ -490,105 +472,89 @@
    */
   public Rotation(RotationOrder order,
                   double alpha1, double alpha2, double alpha3) {
+    this(computeRotation(order, alpha1, alpha2, alpha3));
+  }
 
-    if (order == RotationOrder.XYZ) {
-
-      compose(new Rotation(Vector3D.plusI, alpha1),
-              new Rotation(Vector3D.plusJ, alpha2),
-              new Rotation(Vector3D.plusK, alpha3));
+  /** Copy constructor.
+   * <p>This constructor is used only for the sake of Cardan/Euler
+   * angles handling.</p>
+   * @param r rotation to copy
+   */
+  private Rotation(Rotation r) {
+    q0 = r.q0;
+    q1 = r.q1;
+    q2 = r.q2;
+    q3 = r.q3;
+  }
 
+  /** Build a rotation from three Cardan or Euler elementary rotations.
+   * @param order order of rotations to use
+   * @param alpha1 angle of the first elementary rotation
+   * @param alpha2 angle of the second elementary rotation
+   * @param alpha3 angle of the third elementary rotation
+   */
+  public static Rotation computeRotation(RotationOrder order,
+		                                 double alpha1,
+		                                 double alpha2,
+		                                 double alpha3) {
+    if (order == RotationOrder.XYZ) {
+      return compose(new Rotation(Vector3D.plusI, alpha1),
+                     new Rotation(Vector3D.plusJ, alpha2),
+                     new Rotation(Vector3D.plusK, alpha3));
     } else if (order == RotationOrder.XZY) {
-
-      compose(new Rotation(Vector3D.plusI, alpha1),
-              new Rotation(Vector3D.plusK, alpha2),
-              new Rotation(Vector3D.plusJ, alpha3));
-
+      return compose(new Rotation(Vector3D.plusI, alpha1),
+                     new Rotation(Vector3D.plusK, alpha2),
+                     new Rotation(Vector3D.plusJ, alpha3));
     } else if (order == RotationOrder.YXZ) {
-
-      compose(new Rotation(Vector3D.plusJ, alpha1),
-              new Rotation(Vector3D.plusI, alpha2),
-              new Rotation(Vector3D.plusK, alpha3));
-
+      return compose(new Rotation(Vector3D.plusJ, alpha1),
+                     new Rotation(Vector3D.plusI, alpha2),
+                     new Rotation(Vector3D.plusK, alpha3));
     } else if (order == RotationOrder.YZX) {
-
-      compose(new Rotation(Vector3D.plusJ, alpha1),
-              new Rotation(Vector3D.plusK, alpha2),
-              new Rotation(Vector3D.plusI, alpha3));
-
+      return compose(new Rotation(Vector3D.plusJ, alpha1),
+                     new Rotation(Vector3D.plusK, alpha2),
+                     new Rotation(Vector3D.plusI, alpha3));
     } else if (order == RotationOrder.ZXY) {
-
-      compose(new Rotation(Vector3D.plusK, alpha1),
-              new Rotation(Vector3D.plusI, alpha2),
-              new Rotation(Vector3D.plusJ, alpha3));
-
+      return compose(new Rotation(Vector3D.plusK, alpha1),
+                     new Rotation(Vector3D.plusI, alpha2),
+                     new Rotation(Vector3D.plusJ, alpha3));
     } else if (order == RotationOrder.ZYX) {
-
-      compose(new Rotation(Vector3D.plusK, alpha1),
-              new Rotation(Vector3D.plusJ, alpha2),
-              new Rotation(Vector3D.plusI, alpha3));
-
+     return compose(new Rotation(Vector3D.plusK, alpha1),
+                    new Rotation(Vector3D.plusJ, alpha2),
+                    new Rotation(Vector3D.plusI, alpha3));
     } else if (order == RotationOrder.XYX) {
-
-      compose(new Rotation(Vector3D.plusI, alpha1),
-              new Rotation(Vector3D.plusJ, alpha2),
-              new Rotation(Vector3D.plusI, alpha3));
-
+     return compose(new Rotation(Vector3D.plusI, alpha1),
+                    new Rotation(Vector3D.plusJ, alpha2),
+                    new Rotation(Vector3D.plusI, alpha3));
     } else if (order == RotationOrder.XZX) {
-
-      compose(new Rotation(Vector3D.plusI, alpha1),
-              new Rotation(Vector3D.plusK, alpha2),
-              new Rotation(Vector3D.plusI, alpha3));
-
+     return compose(new Rotation(Vector3D.plusI, alpha1),
+                    new Rotation(Vector3D.plusK, alpha2),
+                    new Rotation(Vector3D.plusI, alpha3));
     } else if (order == RotationOrder.YXY) {
-
-      compose(new Rotation(Vector3D.plusJ, alpha1),
-              new Rotation(Vector3D.plusI, alpha2),
-              new Rotation(Vector3D.plusJ, alpha3));
-
+     return compose(new Rotation(Vector3D.plusJ, alpha1),
+                    new Rotation(Vector3D.plusI, alpha2),
+                    new Rotation(Vector3D.plusJ, alpha3));
     } else if (order == RotationOrder.YZY) {
-
-      compose(new Rotation(Vector3D.plusJ, alpha1),
-              new Rotation(Vector3D.plusK, alpha2),
-              new Rotation(Vector3D.plusJ, alpha3));
-
+     return compose(new Rotation(Vector3D.plusJ, alpha1),
+                    new Rotation(Vector3D.plusK, alpha2),
+                    new Rotation(Vector3D.plusJ, alpha3));
     } else if (order == RotationOrder.ZXZ) {
-
-      compose(new Rotation(Vector3D.plusK, alpha1),
-              new Rotation(Vector3D.plusI, alpha2),
-              new Rotation(Vector3D.plusK, alpha3));
-
+     return compose(new Rotation(Vector3D.plusK, alpha1),
+                    new Rotation(Vector3D.plusI, alpha2),
+                    new Rotation(Vector3D.plusK, alpha3));
     } else { // last possibility is ZYZ
-
-      compose(new Rotation(Vector3D.plusK, alpha1),
-              new Rotation(Vector3D.plusJ, alpha2),
-              new Rotation(Vector3D.plusK, alpha3));
-
+     return compose(new Rotation(Vector3D.plusK, alpha1),
+                    new Rotation(Vector3D.plusJ, alpha2),
+                    new Rotation(Vector3D.plusK, alpha3));
     }
-
   }
 
   /** Override the instance by the composition of three rotations.
-   * @param r1 last (outermost) rotation to compose
+   * @param r3 last (outermost) rotation to compose
    * @param r2 intermediate rotation to compose
-   * @param r3 first (innermost) rotation to compose
-   */
-  private void compose(Rotation r1, Rotation r2, Rotation r3) {
-    Rotation composed = r1.applyTo(r2.applyTo(r3));
-    q0 = composed.q0;
-    q1 = composed.q1;
-    q2 = composed.q2;
-    q3 = composed.q3;
-  }
-
-  /** Copy constructor.
-   * Build a copy of a rotation
-   * @param r rotation to copy
+   * @param r1 first (innermost) rotation to compose
    */
-  public Rotation(Rotation r) {
-    q0 = r.q0;
-    q1 = r.q1;
-    q2 = r.q2;
-    q3 = r.q3;
+  private static Rotation compose(Rotation r3, Rotation r2, Rotation r1) {
+    return r3.applyTo(r2.applyTo(r1));
   }
 
   /** Revert a rotation.
@@ -599,7 +565,7 @@
    * of the instance
    */
   public Rotation revert() {
-    return new Rotation(-q0, q1, q2, q3);
+    return new Rotation(-q0, q1, q2, q3, false);
   }
 
   /** Get the scalar coordinate of the quaternion.
@@ -647,7 +613,7 @@
   }
 
   /** Get the angle of the rotation.
-   * @return angle of the rotation (between 0 and PI)
+   * @return angle of the rotation (between 0 and &pi;)
    */
   public double getAngle() {
     if ((q0 < -0.1) || (q0 > 0.1)) {
@@ -663,14 +629,15 @@
 
    * <p>The equations show that each rotation can be defined by two
    * different values of the Cardan or Euler angles set. For example
-   * if Cardan angles are used, the rotation defined by the angles a1,
-   * a2 and a3 is the same as the rotation defined by the angles PI +
-   * a1, PI - a2 and PI + a3. This method implements the following
-   * arbitrary choices. For Cardan angles, the chosen set is the one
-   * for which the second angle is between -PI/2 and PI/2 (i.e its
-   * cosine is positive). For Euler angles, the chosen set is the one
-   * for which the second angle is between 0 and PI (i.e its sine is
-   * positive).</p>
+   * if Cardan angles are used, the rotation defined by the angles
+   * a<sub>1</sub>, a<sub>2</sub> and a<sub>3</sub> is the same as
+   * the rotation defined by the angles &pi; + a<sub>1</sub>, &pi;
+   * - a<sub>2</sub> and &pi; + a<sub>3</sub>. This method implements
+   * the following arbitrary choices. For Cardan angles, the chosen
+   * set is the one for which the second angle is between -&pi;/2 and
+   * &pi;/2 (i.e its cosine is positive). For Euler angles, the chosen
+   * set is the one for which the second angle is between 0 and &pi;
+   * (i.e its sine is positive).</p>
 
    * <p>Cardan and Euler angle have a very disappointing drawback: all
    * of them have singularities. This means that if the instance is
@@ -678,12 +645,12 @@
    * rotation order, it will be impossible to retrieve the angles. For
    * Cardan angles, this is often called gimbal lock. There is
    * <em>nothing</em> to do to prevent this, it is an intrisic problem
-   * of Cardan and Euler representation (but not a problem with the
+   * with Cardan and Euler representation (but not a problem with the
    * rotation itself, which is perfectly well defined). For Cardan
    * angles, singularities occur when the second angle is close to
-   * -PI/2 or +PI/2, for Euler angle singularities occur when the
-   * second angle is close to 0 or PI, this means that the identity
-   * rotation is always singular for Euler angles !
+   * -&pi;/2 or +&pi;/2, for Euler angle singularities occur when the
+   * second angle is close to 0 or &pi;, this implies that the identity
+   * rotation is always singular for Euler angles!
 
    * @param order rotation order to use
    * @return an array of three angles, in the order specified by the set
@@ -988,7 +955,8 @@
     return new Rotation(r.q0 * q0 - (r.q1 * q1 + r.q2 * q2 + r.q3 * q3),
                         r.q1 * q0 + r.q0 * q1 + (r.q2 * q3 - r.q3 * q2),
                         r.q2 * q0 + r.q0 * q2 + (r.q3 * q1 - r.q1 * q3),
-                        r.q3 * q0 + r.q0 * q3 + (r.q1 * q2 - r.q2 * q1));
+                        r.q3 * q0 + r.q0 * q3 + (r.q1 * q2 - r.q2 * q1),
+                        false);
   }
 
   /** Apply the inverse of the instance to another rotation.
@@ -1006,7 +974,8 @@
     return new Rotation(-r.q0 * q0 - (r.q1 * q1 + r.q2 * q2 + r.q3 * q3),
                         -r.q1 * q0 + r.q0 * q1 + (r.q2 * q3 - r.q3 * q2),
                         -r.q2 * q0 + r.q0 * q2 + (r.q3 * q1 - r.q1 * q3),
-                        -r.q3 * q0 + r.q0 * q3 + (r.q1 * q2 - r.q2 * q1));
+                        -r.q3 * q0 + r.q0 * q3 + (r.q1 * q2 - r.q2 * q1),
+                        false);
   }
 
   /** Perfect orthogonality on a 3X3 matrix.
@@ -1106,35 +1075,17 @@
                                           });
   }
 
-  public int getStateDimension() {
-    return 4;
-  }
-    
-  public void mapStateFromArray(int start, double[] array) {
-    q0 = array[start];
-    q1 = array[start + 1];
-    q2 = array[start + 2];
-    q3 = array[start + 3];
-  }
-
-  public void mapStateToArray(int start, double[] array) {
-    array[start]     = q0;
-    array[start + 1] = q1;
-    array[start + 2] = q2;
-    array[start + 3] = q3;
-  }
-
   /** Scalar coordinate of the quaternion. */
-  private double q0;
+  private final double q0;
 
   /** First coordinate of the vectorial part of the quaternion. */
-  private double q1;
+  private final double q1;
 
   /** Second coordinate of the vectorial part of the quaternion. */
-  private double q2;
+  private final double q2;
 
   /** Third coordinate of the vectorial part of the quaternion. */
-  private double q3;
+  private final double q3;
 
   private static final long serialVersionUID = 7264384082212242475L;
 

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java Wed Nov 22 20:17:13 2006
@@ -18,52 +18,32 @@
 package org.spaceroots.mantissa.geometry;
 
 import java.io.Serializable;
-import org.spaceroots.mantissa.utilities.ArraySliceMappable;
 
 
 /** This class implements vectors in a three-dimensional space.
+ * <p>Vector3D are guaranteed to be immutable objects.</p>
  * @version $Id: Vector3D.java 1705 2006-09-17 19:57:39Z luc $
  * @author L. Maisonobe
  */
+public class Vector3D implements Serializable {
 
-public class Vector3D
-  implements ArraySliceMappable, Serializable {
+  /** First canonical vector (coordinates : 1, 0, 0). */
+  public static final Vector3D plusI = new Vector3D(1, 0, 0);
 
-  /** First canonical vector (coordinates : 1, 0, 0).
-   * This is really an {@link ImmutableVector3D ImmutableVector3D},
-   * hence it can't be changed in any way.
-   */
-  public static final Vector3D plusI = new ImmutableVector3D(1, 0, 0);
-
-  /** Opposite of the first canonical vector (coordinates : -1, 0, 0).
-   * This is really an {@link ImmutableVector3D ImmutableVector3D},
-   * hence it can't be changed in any way.
-   */
-  public static final Vector3D minusI = new ImmutableVector3D(-1, 0, 0);
+  /** Opposite of the first canonical vector (coordinates : -1, 0, 0). */
+  public static final Vector3D minusI = new Vector3D(-1, 0, 0);
 
-  /** Second canonical vector (coordinates : 0, 1, 0).
-   * This is really an {@link ImmutableVector3D ImmutableVector3D},
-   * hence it can't be changed in any way.
-   */
-  public static final Vector3D plusJ = new ImmutableVector3D(0, 1, 0);
+  /** Second canonical vector (coordinates : 0, 1, 0). */
+  public static final Vector3D plusJ = new Vector3D(0, 1, 0);
 
-  /** Opposite of the second canonical vector (coordinates : 0, -1, 0).
-   * This is really an {@link ImmutableVector3D ImmutableVector3D},
-   * hence it can't be changed in any way.
-   */
-  public static final Vector3D minusJ = new ImmutableVector3D(0, -1, 0);
+  /** Opposite of the second canonical vector (coordinates : 0, -1, 0). */
+  public static final Vector3D minusJ = new Vector3D(0, -1, 0);
 
-  /** Third canonical vector (coordinates : 0, 0, 1).
-   * This is really an {@link ImmutableVector3D ImmutableVector3D},
-   * hence it can't be changed in any way.
-   */
-  public static final Vector3D plusK = new ImmutableVector3D(0, 0, 1);
+  /** Third canonical vector (coordinates : 0, 0, 1). */
+  public static final Vector3D plusK = new Vector3D(0, 0, 1);
 
-  /** Opposite of the third canonical vector (coordinates : 0, 0, -1).
-   * This is really an {@link ImmutableVector3D ImmutableVector3D},
-   * hence it can't be changed in any way.
-   */
-  public static final Vector3D minusK = new ImmutableVector3D(0, 0, -1);
+  /** Opposite of the third canonical vector (coordinates : 0, 0, -1). */
+  public static final Vector3D minusK = new Vector3D(0, 0, -1);
 
   /** Simple constructor.
    * Build a null vector.
@@ -140,31 +120,30 @@
    * @param c third scale factor
    * @param w third base (unscaled) vector
    */
-  public Vector3D(double a, Vector3D u,
-                  double b, Vector3D v,
+  public Vector3D(double a, Vector3D u, double b, Vector3D v,
                   double c, Vector3D w) {
     this.x = a * u.x + b * v.x + c * w.x;
     this.y = a * u.y + b * v.y + c * w.y;
     this.z = a * u.z + b * v.z + c * w.z;
   }
 
-  /** Copy constructor.
-   * Build a copy of a vector
-   * @param v vector to copy
-   */
-  public Vector3D(Vector3D v) {
-    x = v.x;
-    y = v.y;
-    z = v.z;
-  }
-
-  /** Reset the instance.
-   * @param v vector to copy data from
-   */
-  public void reset(Vector3D v) {
-    x = v.x;
-    y = v.y;
-    z = v.z;
+  /** Linear constructor
+   * Build a vector from four other ones and corresponding scale factors.
+   * The vector built will be a * t +  b * u + c * v + d * w
+   * @param a first scale factor
+   * @param t first base (unscaled) vector
+   * @param b second scale factor
+   * @param u second base (unscaled) vector
+   * @param c third scale factor
+   * @param v third base (unscaled) vector
+   * @param d third scale factor
+   * @param w third base (unscaled) vector
+   */
+  public Vector3D(double a, Vector3D t, double b, Vector3D u,
+                  double c, Vector3D v, double d, Vector3D w) {
+    this.x = a * t.x + b * u.x + c * v.x + d * w.x;
+    this.y = a * t.y + b * u.y + c * v.y + d * w.y;
+    this.z = a * t.z + b * u.z + c * v.z + d * w.z;
   }
 
   /** Get the abscissa of the vector.
@@ -176,15 +155,6 @@
     return x;
   }
 
-  /** Set the abscissa of the vector.
-   * @param x new abscissa for the vector
-   * @see #getX()
-   * @see #setCoordinates(double, double, double)
-   */
-  public void setX(double x) {
-    this.x = x;
-  }
-
   /** Get the ordinate of the vector.
    * @return ordinate of the vector
    * @see #Vector3D(double, double, double)
@@ -194,15 +164,6 @@
     return y;
   }
 
-  /** Set the ordinate of the vector.
-   * @param y new ordinate for the vector
-   * @see #getY()
-   * @see #setCoordinates(double, double, double)
-   */
-  public void setY(double y) {
-    this.y = y;
-  }
-
   /** Get the height of the vector.
    * @return height of the vector
    * @see #Vector3D(double, double, double)
@@ -212,27 +173,6 @@
     return z;
   }
 
-  /** Set the height of the vector.
-   * @param z new height for the vector
-   * @see #getZ()
-   * @see #setCoordinates(double, double, double)
-   */
-  public void setZ(double z) {
-    this.z = z;
-  }
-
-  /** Set all coordinates of the vector.
-   * @param x new abscissa for the vector
-   * @param y new ordinate for the vector
-   * @param z new height for the vector
-   * @see #Vector3D(double, double, double)
-   */
-  public void setCoordinates(double x, double y, double z) {
-    this.x = x;
-    this.y = y;
-    this.z = z;
-  }
-
   /** Get the norm for the vector.
    * @return euclidian norm for the vector
    */
@@ -256,25 +196,18 @@
     return Math.asin(z / getNorm());
   }
 
-  /** Add a vector to the instance.
-   * Add a vector to the instance. The instance is changed.
-   * @param v vector to add
-   */
-  public void addToSelf(Vector3D v) {
-    x += v.x;
-    y += v.y;
-    z += v.z;
-  }
-
-  /** Add a scaled vector to the instance.
-   * Add a scaled vector to the instance. The instance is changed.
-   * @param factor scale factor to apply to v before adding it
-   * @param v vector to add
-   */
-  public void addToSelf(double factor, Vector3D v) {
-    x += factor * v.x;
-    y += factor * v.y;
-    z += factor * v.z;
+  /** Normalize a vector.
+    * @param v vector to normalize
+   * @return a new vector equal to v / ||v||
+   * @exception ArithmeticException if the norm of the instance is null
+   */
+  public static Vector3D normalize(Vector3D v) {
+	double norm = v.getNorm();
+	if (norm == 0) {
+      throw new ArithmeticException("null norm");
+	}
+	double inv = 1.0 / norm;
+	return new Vector3D(inv * v.x, inv * v.y, inv * v.z);
   }
 
   /** Add two vectors.
@@ -287,16 +220,6 @@
     return new Vector3D(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
   }
 
-  /** Subtract a vector from the instance.
-   * Subtract a vector from the instance. The instance is changed.
-   * @param v vector to subtract
-   */
-  public void subtractFromSelf(Vector3D v) {
-    x -= v.x;
-    y -= v.y;
-    z -= v.z;
-  }
-
   /** Subtract two vectors.
    * Subtract two vectors and return the difference as a new vector
    * @param v1 first vector
@@ -307,22 +230,6 @@
     return new Vector3D(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
   }
 
-  /** Normalize the instance.
-   * Divide the instance by its norm in order to have a unit
-   * vector. The instance is changed.
-   * @exception ArithmeticException if the norm is null
-   */
-  public void normalizeSelf() {
-    double s = getNorm();
-    if (s == 0) {
-      throw new ArithmeticException("null norm");
-    }
-    double invNorm = 1 / s;
-    x *= invNorm;
-    y *= invNorm;
-    z *= invNorm;
-  }
-
   /** Get a vector orthogonal to the instance.
    * <p>There are an infinite number of normalized vectors orthogonal
    * to the instance. This method picks up one of them almost
@@ -331,8 +238,7 @@
    * following example shows how to build a frame having the k axis
    * aligned with the known vector u :
    * <pre><code>
-   *   Vector3D k = u;
-   *   k.normalizeSelf();
+   *   Vector3D k = Vector3D.normalize(u);
    *   Vector3D i = k.orthogonal();
    *   Vector3D j = Vector3D.crossProduct(k, i);
    * </code></pre></p>
@@ -392,15 +298,6 @@
 
   }
 
-  /** Revert the instance.
-   * Replace the instance u by -u
-   */
-  public void negateSelf() {
-    x = -x;
-    y = -y;
-    z = -z;
-  }
-
   /** Get the opposite of a vector.
    * @param u vector to revert
    * @return a new vector which is -u
@@ -409,16 +306,6 @@
     return new Vector3D(-u.x, -u.y, -u.z);
   }
 
-  /** Multiply the instance by a scalar
-   * Multiply the instance by a scalar. The instance is changed.
-   * @param a scalar by which the instance should be multiplied
-   */
-  public void multiplySelf(double a) {
-    x *= a;
-    y *= a;
-    z *= a;
-  }
-
   /** Multiply a vector by a scalar
    * Multiply a vectors by a scalar and return the product as a new vector
    * @param a scalar
@@ -438,18 +325,6 @@
     return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
   }
 
-  /** Set the instance to the result of the cross-product of two vectors.
-   * @param v1 first vector (can be the instance)
-   * @param v2 second vector (can be the instance)
-   */
-  public void setToCrossProduct(Vector3D v1, Vector3D v2) {
-    double newX = v1.y * v2.z - v1.z * v2.y;
-    double newY = v1.z * v2.x - v1.x * v2.z;
-    z = v1.x * v2.y - v1.y * v2.x;
-    x = newX;
-    y = newY;
-  }
-
   /** Compute the cross-product of two vectors.
    * @param v1 first vector
    * @param v2 second vector
@@ -461,31 +336,15 @@
                         v1.x * v2.y - v1.y * v2.x);
   }
 
-  public int getStateDimension() {
-    return 3;
-  }
-    
-  public void mapStateFromArray(int start, double[] array) {
-    x = array[start];
-    y = array[start + 1];
-    z = array[start + 2];
-  }
-
-  public void mapStateToArray(int start, double[] array) {
-    array[start]     = x;
-    array[start + 1] = y;
-    array[start + 2] = z;
-  }
-
   /** Abscissa. */
-  protected double x;
+  private final double x;
 
   /** Ordinate. */
-  protected double y;
+  private final double y;
 
   /** Height. */
-  protected double z;
+  private final double z;
 
-   private static final long serialVersionUID = 4115635019045864211L;
+  private static final long serialVersionUID = 484345009325358136L;
 
 }

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java Wed Nov 22 20:17:13 2006
@@ -880,7 +880,11 @@
         interpolator.storeTime(currentT);
         handler.handleStep(interpolator, lastStep);
 
-        switchesHandler.reset(currentT, y);
+        if (switchesHandler.reset(currentT, y) && ! lastStep) {
+          // some switching function has triggered changes that
+          // invalidate the derivatives, we need to recompute them
+          firstStepAlreadyComputed = false;
+        }
 
         int optimalIter;
         if (k == 1) {

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java Wed Nov 22 20:17:13 2006
@@ -311,7 +311,11 @@
         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
       }
 
-      switchesHandler.reset(currentT, y);
+      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+        // some switching function has triggered changes that
+        // invalidate the derivatives, we need to recompute them
+        equations.computeDerivatives(currentT, y, yDotK[0]);
+      }
 
       if (! lastStep) {
         // stepsize control for next step

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java Wed Nov 22 20:17:13 2006
@@ -234,7 +234,11 @@
         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
       }
 
-      switchesHandler.reset(currentT, y);
+      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+        // some switching function has triggered changes that
+        // invalidate the derivatives, we need to recompute them
+        equations.computeDerivatives(currentT, y, yDotK[0]);
+      }
 
       if (needUpdate) {
         // a switching function has changed the step

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java Wed Nov 22 20:17:13 2006
@@ -239,15 +239,23 @@
    * beginning of the next step
    * @param y array were to put the desired state vector at the beginning
    * of the next step
+   * @return true if the integrator should reset the derivatives too
    */
-  public void reset(double t, double[] y) {
-    if (pendingEvent) {
-      if (nextAction == SwitchingFunction.RESET) {
-        function.resetState(t, y);
-      }
-      pendingEvent      = false;
-      pendingEventTime  = Double.NaN;
+  public boolean reset(double t, double[] y) {
+
+    if (! pendingEvent) {
+      return false;
     }
+
+    if (nextAction == SwitchingFunction.RESET_STATE) {
+      function.resetState(t, y);
+    }
+    pendingEvent      = false;
+    pendingEventTime  = Double.NaN;
+
+    return (nextAction == SwitchingFunction.RESET_STATE)
+        || (nextAction == SwitchingFunction.RESET_DERIVATIVES);
+
   }
 
   /** Get the value of the g function at the specified time.

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java Wed Nov 22 20:17:13 2006
@@ -22,18 +22,20 @@
  * <p>A switching function allows to handle discrete events in
  * integration problems. These events occur for example when the
  * integration process should be stopped as some value is reached
- * (G-stop facility), or when the derivatives has
- * discontinuities. These events are traditionally defined as
- * occurring when a <code>g</code> function sign changes.</p>
+ * (G-stop facility), or when the derivatives have
+ * discontinuities, or simply when the user wants to monitor some
+ * states boundaries crossings. These events are traditionally defined
+ * as occurring when a <code>g</code> function sign changes, hence
+ * the name <em>switching functions</em>.</p>
  *
  * <p>Since events are only problem-dependent and are triggered by the
  * independant <i>time</i> variable and the state vector, they can
- * occur at virtually any time. The integrators will take care to
- * avoid sign changes inside the steps, they will reduce the step size
- * when such an event is detected in order to put this event exactly
- * at the end of the current step. This guarantees that step
- * interpolation (which always has a one step scope) is relevant even
- * in presence of discontinuities. This is independent from the
+ * occur at virtually any time, unknown in advance. The integrators will
+ * take care to avoid sign changes inside the steps, they will reduce
+ * the step size when such an event is detected in order to put this
+ * event exactly at the end of the current step. This guarantees that
+ * step interpolation (which always has a one step scope) is relevant
+ * even in presence of discontinuities. This is independent from the
  * stepsize control provided by integrators that monitor the local
  * error (this feature is available on all integrators, including
  * fixed step ones).</p>
@@ -52,29 +54,38 @@
    */
   public static final int STOP = 0;
 
-  /** Reset indicator.
+  /** Reset state indicator.
    * <p>This value should be used as the return value of the {@link
    * #eventOccurred eventOccurred} method when the integration should
    * go on after the event ending the current step, with a new state
-   * vector (which will be retrieved through the {@link #resetState
+   * vector (which will be retrieved thanks to the {@link #resetState
    * resetState} method).</p>
    */
-  public static final int RESET = 1;
+  public static final int RESET_STATE = 1;
+
+  /** Reset derivatives indicator.
+   * <p>This value should be used as the return value of the {@link
+   * #eventOccurred eventOccurred} method when the integration should
+   * go on after the event ending the current step, with a new derivatives
+   * vector (which will be retrieved thanks to the {@link
+   * FirstOrderDifferentialEquations#computeDerivatives} method).</p>
+   */
+  public static final int RESET_DERIVATIVES = 2;
 
   /** Continue indicator.
    * <p>This value should be used as the return value of the {@link
    * #eventOccurred eventOccurred} method when the integration should go
    * on after the event ending the current step.</p>
    */
-  public static final int CONTINUE = 2;
+  public static final int CONTINUE = 3;
 
   /** Compute the value of the switching function.
 
    * <p>Discrete events are generated when the sign of this function
    * changes, the integrator will take care to change the stepsize in
    * such a way these events occur exactly at step boundaries. This
-   * function must be continuous, as the integrator will need to find
-   * its roots to locate the events.</p>
+   * function must be continuous (at least in its roots neighborhood),
+   * as the integrator will need to find its roots to locate the events.</p>
 
    * @param t current value of the independant <i>time</i> variable
    * @param y array containing the current value of the state vector
@@ -88,22 +99,33 @@
    * ending exactly on a sign change of the function, just before the
    * step handler itself is called. It allows the user to update his
    * internal data to acknowledge the fact the event has been handled
-   * (for example setting a flag to switch the derivatives computation
-   * in case of discontinuity), and it allows to direct the integrator
-   * to either stop or continue integration, possibly with a reset
-   * state.</p>
-
-   * <p>If {@link #STOP} is returned, the step handler will be called
-   * with the <code>isLast</code> flag of the {@link
-   * StepHandler#handleStep handleStep} method set to true. If {@link
-   * #RESET} is returned, the {@link #resetState resetState} method
-   * will be called once the step handler has finished its task.</p>
+   * (for example setting a flag in the {@link
+   * FirstOrderDifferentialEquations differential equations} to switch
+   * the derivatives computation in case of discontinuity), or to
+   * direct the integrator to either stop or continue integration,
+   * possibly with a reset state or derivatives.</p>
+
+   * <ul>
+   *   <li>if {@link #STOP} is returned, the step handler will be called
+   *   with the <code>isLast</code> flag of the {@link
+   *   StepHandler#handleStep handleStep} method set to true and the
+   *   integration will be stopped,</li>
+   *   <li>if {@link #RESET_STATE} is returned, the {@link #resetState
+   *   resetState} method will be called once the step handler has
+   *   finished its task, and the integrator will also recompute the
+   *   derivatives,</li>
+   *   <li>if {@link #RESET_DERIVATIVES} is returned, the integrator
+   *   will recompute the derivatives,
+   *   <li>if {@link #CONTINUE} is returned, no specific action will
+   *   be taken (apart from having called this method) and integration
+   *   will continue.</li>
+   * </ul>
 
    * @param t current value of the independant <i>time</i> variable
    * @param y array containing the current value of the state vector
    * @return indication of what the integrator should do next, this
-   * value must be one of {@link #STOP}, {@link #RESET} or {@link
-   * #CONTINUE}
+   * value must be one of {@link #STOP}, {@link #RESET_STATE},
+   * {@link #RESET_DERIVATIVES} or {@link #CONTINUE}
    */
   public int eventOccurred(double t, double[] y);
   
@@ -111,11 +133,11 @@
 
    * <p>This method is called after the step handler has returned and
    * before the next step is started, but only when {@link
-   * #eventOccurred} has itself returned the {@link #RESET}
+   * #eventOccurred} has itself returned the {@link #RESET_STATE}
    * indicator. It allows the user to reset the state vector for the
    * next step, without perturbing the step handler of the finishing
    * step. If the {@link #eventOccurred} never returns the {@link
-   * #RESET} indicator, this function will never be called, and it is
+   * #RESET_STATE} indicator, this function will never be called, and it is
    * safe to leave its body empty.</p>
 
    * @param t current value of the independant <i>time</i> variable

Modified: jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java Wed Nov 22 20:17:13 2006
@@ -166,11 +166,16 @@
    * beginning of the next step
    * @param y array were to put the desired state vector at the beginning
    * of the next step
+   * @return true if the integrator should reset the derivatives too
    */
-  public void reset(double t, double[] y) {
+  public boolean reset(double t, double[] y) {
+    boolean resetDerivatives = false;
     for (Iterator iter = functions.iterator(); iter.hasNext();) {
-      ((SwitchState) iter.next()).reset(t, y);
+      if (((SwitchState) iter.next()).reset(t, y)) {
+        resetDerivatives = true;
+      }
     }
+    return resetDerivatives;
   }
 
   /** Switching functions. */

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java Wed Nov 22 20:17:13 2006
@@ -26,7 +26,6 @@
     TestSuite suite = new TestSuite("org.spaceroots.mantissa.geometry"); 
 
     suite.addTest(Vector3DTest.suite());
-    suite.addTest(ImmutableVector3DTest.suite());
     suite.addTest(RotationTest.suite());
 
     return suite; 

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java Wed Nov 22 20:17:13 2006
@@ -1,43 +0,0 @@
-// 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.spaceroots.mantissa.geometry;
-
-import junit.framework.*;
-
-public class ImmutableVector3DTest
-  extends TestCase {
-
-  public ImmutableVector3DTest(String name) {
-    super(name);
-  }
-
-  public void testCanonical() {
-    try {
-      Vector3D.plusK.normalizeSelf();
-      fail("an exception should have been thrown");
-    } catch (UnsupportedOperationException uoe) {
-    } catch (Exception e) {
-      fail ("wrong exception caught");
-    }
-  }
-  
-  public static Test suite() {
-    return new TestSuite(ImmutableVector3DTest.class);
-  }
-
-}

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java Wed Nov 22 20:17:13 2006
@@ -43,7 +43,7 @@
 
     Vector3D v1 = new Vector3D(1, 2, 3);
     Vector3D v2 = new Vector3D(-3, -2, -1);
-    v1.subtractFromSelf(v2);
+    v1 = Vector3D.subtract(v1, v2);
     checkVector(v1, new Vector3D(4, 4, 4));
 
     checkVector(Vector3D.subtract(v2, v1), new Vector3D(-7, -6, -5));
@@ -53,7 +53,7 @@
   public void testAdd() {
     Vector3D v1 = new Vector3D(1, 2, 3);
     Vector3D v2 = new Vector3D(-3, -2, -1);
-    v1.addToSelf(v2);
+    v1 = Vector3D.add(v1, v2);
     checkVector(v1, new Vector3D(-2, 0, 2));
 
     checkVector(Vector3D.add(v2, v1), new Vector3D(-5, -2, 1));
@@ -62,7 +62,7 @@
 
   public void testScalarProduct() {
     Vector3D v = new Vector3D(1, 2, 3);
-    v.multiplySelf(3);
+    v = Vector3D.multiply(3, v);
     checkVector(v, new Vector3D(3, 6, 9));
 
     checkVector(Vector3D.multiply(0.5, v), new Vector3D(1.5, 3, 4.5));
@@ -101,12 +101,10 @@
   public void testAngularSeparation() {
     Vector3D v1 = new Vector3D(2, -1, 4);
 
-    Vector3D  k = v1;
-    k.normalizeSelf();
+    Vector3D  k = Vector3D.normalize(v1);
     Vector3D  i = k.orthogonal();
 
-    Vector3D v2 = Vector3D.multiply(Math.cos(1.2), k);
-    v2.addToSelf(Vector3D.multiply(Math.sin(1.2), i));
+    Vector3D v2 = new Vector3D(Math.cos(1.2), k, Math.sin(1.2), i);
 
     assertTrue(Math.abs(Vector3D.angle(v1, v2) - 1.2) < 1.0e-12);
 

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java Wed Nov 22 20:17:13 2006
@@ -305,6 +305,17 @@
 
   }
 
+  public void testUnstableDerivative()
+  throws DerivativeException, IntegratorException {
+    final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+    FirstOrderIntegrator integ =
+      new DormandPrince853Integrator(0.1, 10, 1.0e-12, 0.0);
+    integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+    double[] y = { Double.NaN };
+    integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+    assertEquals(8.0, y[0], 1.0e-12);
+  }
+
   public static Test suite() {
     return new TestSuite(DormandPrince853IntegratorTest.class);
   }

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java Wed Nov 22 20:17:13 2006
@@ -196,6 +196,16 @@
                     pb.getFinalTime(), new double[pb.getDimension()]);
   }
 
+  public void testUnstableDerivative()
+  throws DerivativeException, IntegratorException {
+    final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+    FirstOrderIntegrator integ = new GillIntegrator(0.3);
+    integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+    double[] y = { Double.NaN };
+    integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+    assertEquals(8.0, y[0], 1.0e-12);
+  }
+
   public static Test suite() {
     return new TestSuite(GillIntegratorTest.class);
   }

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java Wed Nov 22 20:17:13 2006
@@ -256,6 +256,17 @@
                     pb.getFinalTime(), new double[pb.getDimension()]);
   }
 
+  public void testUnstableDerivative()
+    throws DerivativeException, IntegratorException {
+    final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+    FirstOrderIntegrator integ =
+      new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
+    integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+    double[] y = { Double.NaN };
+    integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+    assertEquals(8.0, y[0], 1.0e-12);
+  }
+
   public static Test suite() {
     return new TestSuite(GraggBulirschStoerIntegratorTest.class);
   }

Modified: jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java?view=diff&rev=478458&r1=478457&r2=478458
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java (original)
+++ jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java Wed Nov 22 20:17:13 2006
@@ -104,7 +104,7 @@
     public int eventOccurred(double t, double[] y) {
       // this sign change is needed because the state will be reset soon
       sign = -sign;
-      return SwitchingFunction.RESET;
+      return SwitchingFunction.RESET_STATE;
     }
   
     public void resetState(double t, double[] y) {



---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message