Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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
This commit is contained in:
Phil Steitz 2006-11-23 04:17:13 +00:00
parent 05195b77ca
commit 1d1e19921c
17 changed files with 330 additions and 691 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java
View File

@ -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}" },

204
src/mantissa/src/org/spaceroots/mantissa/geometry/ImmutableVector3D.java
View File

@ -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;
}

363
src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java
View File

@ -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;
* <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>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>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>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 @@ 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.
* <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 @@ public class Rotation
/** 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 @@ public class Rotation
/** 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>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 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>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 @@ 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).</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,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
* <p>This constructor is used only for the sake of Cardan/Euler
* angles handling.</p>
* @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 &pi;)
*/
public double getAngle() {
if ((q0 < -0.1) || (q0 > 0.1)) {
@ -663,14 +629,15 @@ public class Rotation
* <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 @@ public class Rotation
* 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 @@ 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;

235
src/mantissa/src/org/spaceroots/mantissa/geometry/Vector3D.java
View File

@ -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). */
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.
* <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 @@ public class Vector3D
* 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 @@ 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;
}

6
src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java
View File

@ -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) {

6
src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java
View File

@ -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

6
src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java
View File

@ -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

22
src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java
View File

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

78
src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java
View File

@ -22,18 +22,20 @@ package org.spaceroots.mantissa.ode;
* <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 interface SwitchingFunction {
*/
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 @@ 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.</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>
* <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>
* <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 @@ public interface SwitchingFunction {
* <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

9
src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java
View File

@ -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. */

1
src/mantissa/tests-src/org/spaceroots/mantissa/geometry/AllTests.java
View File

@ -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;

43
src/mantissa/tests-src/org/spaceroots/mantissa/geometry/ImmutableVector3DTest.java
View File

@ -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);
}
}

12
src/mantissa/tests-src/org/spaceroots/mantissa/geometry/Vector3DTest.java
View File

@ -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);

11
src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java
View File

@ -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);
}

10
src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java
View File

@ -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);
}

11
src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java
View File

@ -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);
}

2
src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java
View File

@ -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) {