From 1d1e19921ccc1ec10ba3a20212efc1f5c1d0f155 Mon Sep 17 00:00:00 2001 From: Phil Steitz Date: Thu, 23 Nov 2006 04:17:13 +0000 Subject: [PATCH] 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 git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@478458 13f79535-47bb-0310-9956-ffa450edef68 --- .../mantissa/MessagesResources_fr.java | 2 +- .../mantissa/geometry/ImmutableVector3D.java | 204 ---------- .../mantissa/geometry/Rotation.java | 363 ++++++++---------- .../mantissa/geometry/Vector3D.java | 235 +++--------- .../ode/GraggBulirschStoerIntegrator.java | 6 +- .../ode/RungeKuttaFehlbergIntegrator.java | 6 +- .../mantissa/ode/RungeKuttaIntegrator.java | 6 +- .../spaceroots/mantissa/ode/SwitchState.java | 22 +- .../mantissa/ode/SwitchingFunction.java | 78 ++-- .../ode/SwitchingFunctionsHandler.java | 9 +- .../mantissa/geometry/AllTests.java | 1 - .../geometry/ImmutableVector3DTest.java | 43 --- .../mantissa/geometry/Vector3DTest.java | 12 +- .../ode/DormandPrince853IntegratorTest.java | 11 + .../mantissa/ode/GillIntegratorTest.java | 10 + .../ode/GraggBulirschStoerIntegratorTest.java | 11 + .../spaceroots/mantissa/ode/TestProblem4.java | 2 +- 17 files changed, 330 insertions(+), 691 deletions(-) diff --git a/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java b/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java index cd166354d..47aba4123 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java +++ b/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java @@ -91,7 +91,7 @@ public class MessagesResources_fr { "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}" }, diff --git a/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java b/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java index 54194f268..e69de29bb 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java +++ b/src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java @@ -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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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 UnsupportedOperationException 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; - -} diff --git a/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java b/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java index 5e84a038b..4a4267ea6 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java +++ b/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java @@ -17,7 +17,6 @@ package org.spaceroots.mantissa.geometry; -import org.spaceroots.mantissa.utilities.ArraySliceMappable; import java.io.Serializable; /** @@ -25,59 +24,61 @@ import java.io.Serializable; *

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 * Rotation 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:

+ * from a set of vectors and their image.

+ *

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:

* 
  * double[] angles = new Rotation(matrix, 1.0e-10).getAngles(RotationOrder.XYZ);
  *
- *

Focus is more oriented on what a rotation do. Once it - * has been built, and regardless of its representation, a rotation is - * an operator 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.

+ *

Focus is oriented on what a rotation do rather than on its + * underlying representation. Once it has been built, and regardless of its + * internal representation, a rotation is an operator 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.

+ *

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.

- *

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 - * projectVectorIntoDestinationFrame or - * computeTransformedDirection. It provides simpler and - * more generic methods: {@link #applyTo(Vector3D) applyTo(Vector3D)} - * and {@link #applyInverseTo(Vector3D) applyInverseTo(Vector3D)}.

+ *

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 projectVectorIntoDestinationFrame or + * computeTransformedDirection. It provides simpler and more generic + * methods: {@link #applyTo(Vector3D) applyTo(Vector3D)} and {@link + * #applyInverseTo(Vector3D) applyInverseTo(Vector3D)}.

- *

Since a rotation is basically a vectorial operator, several - * rotations can be composed together and the composite operation - * r = r1 o r2 (which means that for each vector - * u, r(u) = r1(r2(u))) 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 r1 to - * r2 and the result we get is r = r1 o - * r2. For this purpose, the class provides the methods: {@link - * #applyTo(Rotation) applyTo(Rotation)} and {@link - * #applyInverseTo(Rotation) applyInverseTo(Rotation)}.

+ *

Since a rotation is basically a vectorial operator, several rotations can be + * composed together and the composite operation r = r1 o + * r2 (which means that for each vector u, + * r(u) = r1(r2(u))) 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 r1 to r2 and the result + * we get is r = r1 o r2. For this purpose, the + * class provides the methods: {@link #applyTo(Rotation) applyTo(Rotation)} and + * {@link #applyInverseTo(Rotation) applyInverseTo(Rotation)}.

+ + *

Rotations are guaranteed to be immutable objects.

* @version $Id: Rotation.java 1705 2006-09-17 19:57:39Z luc $ * @author L. Maisonobe @@ -86,8 +87,7 @@ import java.io.Serializable; */ -public class Rotation - implements ArraySliceMappable, Serializable { +public class Rotation implements Serializable { /** Build the identity rotation. */ @@ -98,31 +98,12 @@ public class Rotation q3 = 0; } - /** 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. *

A rotation can be built from a normalized quaternion, * i.e. a quaternion for which q02 + * q12 + q22 + * q32 = 1. If the quaternion is not normalized, * the constructor can normalize it in a preprocessing step.

- *

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.

* @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 @@ public class Rotation /** Build a rotation from a 3X3 matrix. *

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.mT = 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.

@@ -290,14 +271,15 @@ public class Rotation /** Build the rotation that transforms a pair of vector into another pair. - *

Except for possible scale factors, if the instance were - * applied to the pair (u1, u2) it will produce the pair (v1, - * v2).

+ *

Except for possible scale factors, if the instance were applied to + * the pair (u1, u2) it will produce the pair + * (v1, v2).

- *

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.

+ *

If the angular separation between u1 and u2 is + * not the same as the angular separation between v1 and + * v2, then a corrected v'2 will be used rather than + * v2, the corrected vector will be in the (v1, + * v2) plane.

* @param u1 first vector of the origin pair * @param u2 second vector of the origin pair @@ -477,11 +459,11 @@ public class Rotation * 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).

+ * 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).

* @param order order of rotations to use * @param alpha1 angle of the first elementary rotation @@ -490,107 +472,91 @@ public class Rotation */ public Rotation(RotationOrder order, double alpha1, double alpha2, double alpha3) { - - if (order == RotationOrder.XYZ) { - - 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)); - - } else if (order == RotationOrder.YXZ) { - - 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)); - - } else if (order == RotationOrder.ZXY) { - - 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)); - - } else if (order == RotationOrder.XYX) { - - 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)); - - } else if (order == RotationOrder.YXY) { - - 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)); - - } else if (order == RotationOrder.ZXZ) { - - 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)); - - } - - } - - /** Override the instance by the composition of three rotations. - * @param r1 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; + this(computeRotation(order, alpha1, alpha2, alpha3)); } /** Copy constructor. - * Build a copy of a rotation + *

This constructor is used only for the sake of Cardan/Euler + * angles handling.

* @param r rotation to copy */ - public Rotation(Rotation r) { + 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) { + return compose(new Rotation(Vector3D.plusI, alpha1), + new Rotation(Vector3D.plusK, alpha2), + new Rotation(Vector3D.plusJ, alpha3)); + } else if (order == RotationOrder.YXZ) { + return compose(new Rotation(Vector3D.plusJ, alpha1), + new Rotation(Vector3D.plusI, alpha2), + new Rotation(Vector3D.plusK, alpha3)); + } else if (order == RotationOrder.YZX) { + return compose(new Rotation(Vector3D.plusJ, alpha1), + new Rotation(Vector3D.plusK, alpha2), + new Rotation(Vector3D.plusI, alpha3)); + } else if (order == RotationOrder.ZXY) { + return compose(new Rotation(Vector3D.plusK, alpha1), + new Rotation(Vector3D.plusI, alpha2), + new Rotation(Vector3D.plusJ, alpha3)); + } else if (order == RotationOrder.ZYX) { + return compose(new Rotation(Vector3D.plusK, alpha1), + new Rotation(Vector3D.plusJ, alpha2), + new Rotation(Vector3D.plusI, alpha3)); + } else if (order == RotationOrder.XYX) { + return compose(new Rotation(Vector3D.plusI, alpha1), + new Rotation(Vector3D.plusJ, alpha2), + new Rotation(Vector3D.plusI, alpha3)); + } else if (order == RotationOrder.XZX) { + return compose(new Rotation(Vector3D.plusI, alpha1), + new Rotation(Vector3D.plusK, alpha2), + new Rotation(Vector3D.plusI, alpha3)); + } else if (order == RotationOrder.YXY) { + return compose(new Rotation(Vector3D.plusJ, alpha1), + new Rotation(Vector3D.plusI, alpha2), + new Rotation(Vector3D.plusJ, alpha3)); + } else if (order == RotationOrder.YZY) { + return compose(new Rotation(Vector3D.plusJ, alpha1), + new Rotation(Vector3D.plusK, alpha2), + new Rotation(Vector3D.plusJ, alpha3)); + } else if (order == RotationOrder.ZXZ) { + return compose(new Rotation(Vector3D.plusK, alpha1), + new Rotation(Vector3D.plusI, alpha2), + new Rotation(Vector3D.plusK, alpha3)); + } else { // last possibility is ZYZ + 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 r3 last (outermost) rotation to compose + * @param r2 intermediate rotation to compose + * @param r1 first (innermost) rotation to compose + */ + private static Rotation compose(Rotation r3, Rotation r2, Rotation r1) { + return r3.applyTo(r2.applyTo(r1)); + } + /** Revert a rotation. * Build a rotation which reverse the effect of another * rotation. This means that is r(u) = v, then r.revert (v) = u. The @@ -599,7 +565,7 @@ public class Rotation * 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 @@ public class Rotation } /** Get the angle of the rotation. - * @return angle of the rotation (between 0 and PI) + * @return angle of the rotation (between 0 and π) */ public double getAngle() { if ((q0 < -0.1) || (q0 > 0.1)) { @@ -663,14 +629,15 @@ public class Rotation *

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).

+ * 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 π + a1, π + * - a2 and π + a3. This method implements + * the following arbitrary choices. For Cardan angles, the chosen + * set is the one for which the second angle is between -π/2 and + * π/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 π + * (i.e its sine is positive).

*

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 @@ public class Rotation * rotation order, it will be impossible to retrieve the angles. For * Cardan angles, this is often called gimbal lock. There is * nothing 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 ! + * -π/2 or +π/2, for Euler angle singularities occur when the + * second angle is close to 0 or π, 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 @@ public class Rotation 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 @@ public class Rotation 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 class Rotation }); } - 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; diff --git a/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java b/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java index 8bd96d953..1a985ad4b 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java +++ b/src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java @@ -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. + *

Vector3D are guaranteed to be immutable objects.

* @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). */ + public static final Vector3D minusI = new Vector3D(-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); + /** Second canonical vector (coordinates : 0, 1, 0). */ + public static final Vector3D plusJ = new Vector3D(0, 1, 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); + /** Opposite of the second canonical vector (coordinates : 0, -1, 0). */ + public static final Vector3D minusJ = 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); + /** Third canonical vector (coordinates : 0, 0, 1). */ + public static final Vector3D plusK = new Vector3D(0, 0, 1); - /** 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); - - /** 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 @@ public class Vector3D * @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 + /** 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(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; + 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 @@ public class Vector3D 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 @@ public class Vector3D 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 @@ public class Vector3D 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 @@ public class Vector3D 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 + /** 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 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; + 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 @@ public class Vector3D 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 @@ public class Vector3D 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. *

There are an infinite number of normalized vectors orthogonal * to the instance. This method picks up one of them almost @@ -331,8 +238,7 @@ public class Vector3D * following example shows how to build a frame having the k axis * aligned with the known vector u : * 


-   *   Vector3D k = u;
-   *   k.normalizeSelf();
+   *   Vector3D k = Vector3D.normalize(u);
    *   Vector3D i = k.orthogonal();
    *   Vector3D j = Vector3D.crossProduct(k, i);
    *

@@ -392,15 +298,6 @@ public class Vector3D } - /** 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 @@ public class Vector3D 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 @@ public class Vector3D 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 @@ public class Vector3D 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; } diff --git a/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java b/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java index 32c182986..bbc18de5f 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java +++ b/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java @@ -880,7 +880,11 @@ public class GraggBulirschStoerIntegrator 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) { diff --git a/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java b/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java index e27bbb824..05eaa1da8 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java +++ b/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java @@ -311,7 +311,11 @@ public abstract class RungeKuttaFehlbergIntegrator 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 diff --git a/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java b/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java index 1705f6b66..c44234a48 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java +++ b/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java @@ -234,7 +234,11 @@ public abstract class RungeKuttaIntegrator 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 diff --git a/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java b/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java index 7c3879285..1fd262124 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java +++ b/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java @@ -239,15 +239,23 @@ class SwitchState * 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. diff --git a/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java b/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java index cf254eb8b..c0e231fd7 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java +++ b/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java @@ -22,18 +22,20 @@ package org.spaceroots.mantissa.ode; *

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 g function sign changes.

+ * (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 g function sign changes, hence + * the name switching functions.

* *

Since events are only problem-dependent and are triggered by the * independant time 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).

@@ -52,29 +54,38 @@ public interface SwitchingFunction { */ public static final int STOP = 0; - /** Reset indicator. + /** Reset state indicator. *

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).

*/ - public static final int RESET = 1; + public static final int RESET_STATE = 1; + + /** Reset derivatives indicator. + *

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).

+ */ + public static final int RESET_DERIVATIVES = 2; /** Continue indicator. *

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.

*/ - public static final int CONTINUE = 2; + public static final int CONTINUE = 3; /** Compute the value of the switching function. *

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.

+ * function must be continuous (at least in its roots neighborhood), + * as the integrator will need to find its roots to locate the events.

* @param t current value of the independant time variable * @param y array containing the current value of the state vector @@ -88,22 +99,33 @@ public interface SwitchingFunction { * 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.

+ * (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.

- *

If {@link #STOP} is returned, the step handler will be called - * with the isLast 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.

+ * * @param t current value of the independant time 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 @@ public interface SwitchingFunction { *

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.

* @param t current value of the independant time variable diff --git a/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java b/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java index 3c7a1756e..0c459a76e 100644 --- a/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java +++ b/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java @@ -166,11 +166,16 @@ public class SwitchingFunctionsHandler { * 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. */ diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java b/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java index e1ce528cc..ab9082b3e 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java @@ -26,7 +26,6 @@ public class AllTests { TestSuite suite = new TestSuite("org.spaceroots.mantissa.geometry"); suite.addTest(Vector3DTest.suite()); - suite.addTest(ImmutableVector3DTest.suite()); suite.addTest(RotationTest.suite()); return suite; diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java b/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java index b378ef0c9..e69de29bb 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java @@ -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); - } - -} diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java b/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java index 576379fea..324c4ca14 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java @@ -43,7 +43,7 @@ public class Vector3DTest 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 class Vector3DTest 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 class Vector3DTest 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 class Vector3DTest 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); diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java index 3cefb9a07..29287b892 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java @@ -305,6 +305,17 @@ public class DormandPrince853IntegratorTest } + 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); } diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java index 5aa0272d7..26d3c75f9 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java @@ -196,6 +196,16 @@ public class GillIntegratorTest 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); } diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java index 113fb7aa5..b261c4346 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java @@ -256,6 +256,17 @@ public class GraggBulirschStoerIntegratorTest 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); } diff --git a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java index 4fe511b22..ec9b3ea47 100644 --- a/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java +++ b/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java @@ -104,7 +104,7 @@ class TestProblem4 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) {