Author: luc
Date: Wed Feb 6 20:08:33 2013
New Revision: 1443178
URL: http://svn.apache.org/viewvc?rev=1443178&view=rev
Log:
Added conversion of gradients and Hessians from spherical to Cartesian
coordinates in 3D.
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java
(with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java
(with props)
Modified:
commons/proper/math/trunk/src/changes/changes.xml
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1443178&r1=1443177&r2=1443178&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Wed Feb 6 20:08:33 2013
@@ -55,6 +55,10 @@ This is a minor release: It combines bug
Changes to existing features were made in a backwards-compatible
way such as to allow drop-in replacement of the v3.1[.1] JAR file.
">
+ <action dev="luc" type="add" >
+ Added conversion of gradients and Hessians from spherical to Cartesian
+ coordinates in 3D.
+ </action>
<action dev="erans" type="update" issue="MATH-931" due-to="Sean Owen">
Greater efficiency in "UnitSphereRandomVectorGenerator".
</action>
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java?rev=1443178&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java
Wed Feb 6 20:08:33 2013
@@ -0,0 +1,396 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.geometry.euclidean.threed;
+
+
+import java.io.Serializable;
+
+import org.apache.commons.math3.util.FastMath;
+
+/** This class provides conversions related to <a
+ * href="http://mathworld.wolfram.com/SphericalCoordinates.html";>spherical
coordinates</a>.
+ * <p>
+ * The conventions used here are the mathematical ones, i.e. spherical
coordinates are
+ * related to Cartesian coordinates as follows:
+ * </p>
+ * <ul>
+ * <li>x = r cos(θ) sin(Φ)</li>
+ * <li>y = r sin(θ) sin(Φ)</li>
+ * <li>z = r cos(Φ)</li>
+ * </ul>
+ * <ul>
+ * <li>r = √(x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>)</li>
+ * <li>θ = atan2(y, x)</li>
+ * <li>Φ = acos(z/r)</li>
+ * </ul>
+ * <p>
+ * r is the radius, θ is the azimuthal angle in the x-y plane and Φ
is the polar
+ * (co-latitude) angle. These conventions are <em>different</em> from the
conventions used
+ * in physics (and in particular in spherical harmonics) where the meanings of
θ and
+ * Φ are reversed.
+ * </p>
+ * <p>
+ * This class provides conversion of coordinates and also of gradient and
Hessian
+ * between spherical and Cartesian coordinates.
+ * </p>
+ * @since 3.2
+ * @version $Id$
+ */
+public class SphericalCoordinates implements Serializable {
+
+ /** Serializable UID. */
+ private static final long serialVersionUID = 20130206L;
+
+ /** Cartesian coordinates. */
+ private final Vector3D v;
+
+ /** Radius. */
+ private final double r;
+
+ /** Azimuthal angle in the x-y plane θ. */
+ private final double theta;
+
+ /** Polar angle (co-latitude) φ. */
+ private final double phi;
+
+ /** Jacobian of (r, θ &phi). */
+ private double[][] jacobian;
+
+ /** Hessian of radius. */
+ private double[][] rHessian;
+
+ /** Hessian of azimuthal angle in the x-y plane θ. */
+ private double[][] thetaHessian;
+
+ /** Hessian of polar (co-latitude) angle Φ. */
+ private double[][] phiHessian;
+
+ /** Build a spherical coordinates transformer from Cartesian coordinates.
+ * @param v Cartesian coordinates
+ */
+ public SphericalCoordinates(final Vector3D v) {
+
+ // Cartesian coordinates
+ this.v = v;
+
+ // remaining spherical coordinates
+ this.r = v.getNorm();
+ this.theta = v.getAlpha();
+ this.phi = FastMath.acos(v.getZ() / r);
+
+ }
+
+ /** Build a spherical coordinates transformer from spherical coordinates.
+ * @param r radius
+ * @param theta azimuthal angle in x-y place
+ * @param phi polar (co-latitude) angle
+ */
+ public SphericalCoordinates(final double r, final double theta, final
double phi) {
+
+ final double cosTheta = FastMath.cos(theta);
+ final double sinTheta = FastMath.sin(theta);
+ final double cosPhi = FastMath.cos(phi);
+ final double sinPhi = FastMath.sin(phi);
+
+ // spherical coordinates
+ this.r = r;
+ this.theta = theta;
+ this.phi = phi;
+
+ // Cartesian coordinates
+ this.v = new Vector3D(r * cosTheta * sinPhi,
+ r * sinTheta * sinPhi,
+ r * cosPhi);
+
+ }
+
+ /** Get the Cartesian coordinates.
+ * @return Cartesian coordinates
+ */
+ public Vector3D getCartesian() {
+ return v;
+ }
+
+ /** Get the radius.
+ * @return radius r
+ * @see #getTheta()
+ * @see #getPhi()
+ */
+ public double getR() {
+ return r;
+ }
+
+ /** Get the azimuthal angle in x-y plane.
+ * @return azimuthal angle in x-y plane θ
+ * @see #getR()
+ * @see #getPhi()
+ */
+ public double getTheta() {
+ return theta;
+ }
+
+ /** Get the polar (co-latitude) angle.
+ * @return polar (co-latitude) angle Φ
+ * @see #getR()
+ * @see #getTheta()
+ */
+ public double getPhi() {
+ return phi;
+ }
+
+ /** Convert a gradient with respect to spherical coordinates into a
gradient
+ * with respect to Cartesian coordinates.
+ * @param sGradient gradient with respect to spherical coordinates
+ * {df/dr, df/dθ, df/dΦ}
+ * @return gradient with respect to Cartesian coordinates
+ * {df/dx, df/dy, df/dz}
+ */
+ public double[] toCartesianGradient(final double[] sGradient) {
+
+ // lazy evaluation of Jacobian
+ computeJacobian();
+
+ // compose derivatives as gradient^T . J
+ // the expressions have been simplified since we know jacobian[1][2] =
dTheta/dZ = 0
+ return new double[] {
+ sGradient[0] * jacobian[0][0] + sGradient[1] * jacobian[1][0] +
sGradient[2] * jacobian[2][0],
+ sGradient[0] * jacobian[0][1] + sGradient[1] * jacobian[1][1] +
sGradient[2] * jacobian[2][1],
+ sGradient[0] * jacobian[0][2] +
sGradient[2] * jacobian[2][2]
+ };
+
+ }
+
+ /** Convert a Hessian with respect to spherical coordinates into a Hessian
+ * with respect to Cartesian coordinates.
+ * <p>
+ * As Hessian are always symmetric, we use only the lower left part of the
provided
+ * spherical Hessian, so the upper part may not be initialized. However,
we still
+ * do fill up the complete array we create, with guaranteed symmetry.
+ * </p>
+ * @param sHessian Hessian with respect to spherical coordinates
+ * {{d<sup>2</sup>f/dr<sup>2</sup>, d<sup>2</sup>f/drdθ,
d<sup>2</sup>f/drdΦ},
+ * {d<sup>2</sup>f/drdθ, d<sup>2</sup>f/dθ<sup>2</sup>,
d<sup>2</sup>f/dθdΦ},
+ * {d<sup>2</sup>f/drdΦ, d<sup>2</sup>f/dθdΦ,
d<sup>2</sup>f/dΦ<sup>2</sup>}
+ * @param sGradient gradient with respect to spherical coordinates
+ * {df/dr, df/dθ, df/dΦ}
+ * @return Hessian with respect to Cartesian coordinates
+ * {{d<sup>2</sup>f/dx<sup>2</sup>, d<sup>2</sup>f/rGradient.getY(),
d<sup>2</sup>f/dxdz},
+ * {d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dy<sup>2</sup>,
d<sup>2</sup>f/dydz},
+ * {d<sup>2</sup>f/dxdz, d<sup>2</sup>f/dydz,
d<sup>2</sup>f/dz<sup>2</sup>}}
+ */
+ public double[][] toCartesianHessian(final double[][] sHessian, final
double[] sGradient) {
+
+ computeJacobian();
+ computeHessians();
+
+ // compose derivative as J^T . H_f . J + df/dr H_r + df/dtheta H_theta
+ df/dphi H_phi
+ // the expressions have been simplified since we know jacobian[1][2] =
dTheta/dZ = 0
+ // and H_theta is only a 2x2 matrix as it does not depend on z
+ final double[][] hj = new double[3][3];
+ final double[][] cHessian = new double[3][3];
+
+ // compute H_f . J
+ // beware we use ONLY the lower-left part of sHessian
+ hj[0][0] = sHessian[0][0] * jacobian[0][0] + sHessian[1][0] *
jacobian[1][0] + sHessian[2][0] * jacobian[2][0];
+ hj[0][1] = sHessian[0][0] * jacobian[0][1] + sHessian[1][0] *
jacobian[1][1] + sHessian[2][0] * jacobian[2][1];
+ hj[0][2] = sHessian[0][0] * jacobian[0][2]
+ sHessian[2][0] * jacobian[2][2];
+ hj[1][0] = sHessian[1][0] * jacobian[0][0] + sHessian[1][1] *
jacobian[1][0] + sHessian[2][1] * jacobian[2][0];
+ hj[1][1] = sHessian[1][0] * jacobian[0][1] + sHessian[1][1] *
jacobian[1][1] + sHessian[2][1] * jacobian[2][1];
+ hj[1][2] = sHessian[1][0] * jacobian[0][2]
+ sHessian[2][1] * jacobian[2][2];
+ hj[2][0] = sHessian[2][0] * jacobian[0][0] + sHessian[2][1] *
jacobian[1][0] + sHessian[2][2] * jacobian[2][0];
+ hj[2][1] = sHessian[2][0] * jacobian[0][1] + sHessian[2][1] *
jacobian[1][1] + sHessian[2][2] * jacobian[2][1];
+ hj[2][2] = sHessian[2][0] * jacobian[0][2]
+ sHessian[2][2] * jacobian[2][2];
+
+ // compute lower-left part of J^T . H_f . J
+ cHessian[0][0] = jacobian[0][0] * hj[0][0] + jacobian[1][0] * hj[1][0]
+ jacobian[2][0] * hj[2][0];
+ cHessian[1][0] = jacobian[0][1] * hj[0][0] + jacobian[1][1] * hj[1][0]
+ jacobian[2][1] * hj[2][0];
+ cHessian[2][0] = jacobian[0][2] * hj[0][0]
+ jacobian[2][2] * hj[2][0];
+ cHessian[1][1] = jacobian[0][1] * hj[0][1] + jacobian[1][1] * hj[1][1]
+ jacobian[2][1] * hj[2][1];
+ cHessian[2][1] = jacobian[0][2] * hj[0][1]
+ jacobian[2][2] * hj[2][1];
+ cHessian[2][2] = jacobian[0][2] * hj[0][2]
+ jacobian[2][2] * hj[2][2];
+
+ // add gradient contribution
+ cHessian[0][0] += sGradient[0] * rHessian[0][0] + sGradient[1] *
thetaHessian[0][0] + sGradient[2] * phiHessian[0][0];
+ cHessian[1][0] += sGradient[0] * rHessian[1][0] + sGradient[1] *
thetaHessian[1][0] + sGradient[2] * phiHessian[1][0];
+ cHessian[2][0] += sGradient[0] * rHessian[2][0]
+ sGradient[2] * phiHessian[2][0];
+ cHessian[1][1] += sGradient[0] * rHessian[1][1] + sGradient[1] *
thetaHessian[1][1] + sGradient[2] * phiHessian[1][1];
+ cHessian[2][1] += sGradient[0] * rHessian[2][1]
+ sGradient[2] * phiHessian[2][1];
+ cHessian[2][2] += sGradient[0] * rHessian[2][2]
+ sGradient[2] * phiHessian[2][2];
+
+ // ensure symmetry
+ cHessian[0][1] = cHessian[1][0];
+ cHessian[0][2] = cHessian[2][0];
+ cHessian[1][2] = cHessian[2][1];
+
+ return cHessian;
+
+ }
+
+ /** Lazy evaluation of (r, θ, φ) Jacobian.
+ */
+ private void computeJacobian() {
+ if (jacobian == null) {
+
+ // intermediate variables
+ final double x = v.getX();
+ final double y = v.getY();
+ final double z = v.getZ();
+ final double rho2 = x * x + y * y;
+ final double rho = FastMath.sqrt(rho2);
+ final double r2 = rho2 + z * z;
+
+ jacobian = new double[3][3];
+
+ // row representing the gradient of r
+ jacobian[0][0] = x / r;
+ jacobian[0][1] = y / r;
+ jacobian[0][2] = z / r;
+
+ // row representing the gradient of theta
+ jacobian[1][0] = -y / rho2;
+ jacobian[1][1] = x / rho2;
+ // jacobian[1][2] is already set to 0 at allocation time
+
+ // row representing the gradient of phi
+ jacobian[2][0] = x * z / (rho * r2);
+ jacobian[2][1] = y * z / (rho * r2);
+ jacobian[2][2] = -rho / r2;
+
+ }
+ }
+
+ /** Lazy evaluation of Hessians.
+ */
+ private void computeHessians() {
+
+ if (rHessian == null) {
+
+ // intermediate variables
+ final double x = v.getX();
+ final double y = v.getY();
+ final double z = v.getZ();
+ final double x2 = x * x;
+ final double y2 = y * y;
+ final double z2 = z * z;
+ final double rho2 = x2 + y2;
+ final double rho = FastMath.sqrt(rho2);
+ final double r2 = rho2 + z2;
+ final double xOr = x / r;
+ final double yOr = y / r;
+ final double zOr = z / r;
+ final double xOrho2 = x / rho2;
+ final double yOrho2 = y / rho2;
+ final double xOr3 = xOr / r2;
+ final double yOr3 = yOr / r2;
+ final double zOr3 = zOr / r2;
+
+ // lower-left part of Hessian of r
+ rHessian = new double[3][3];
+ rHessian[0][0] = y * yOr3 + z * zOr3;
+ rHessian[1][0] = -x * yOr3;
+ rHessian[2][0] = -z * xOr3;
+ rHessian[1][1] = x * xOr3 + z * zOr3;
+ rHessian[2][1] = -y * zOr3;
+ rHessian[2][2] = x * xOr3 + y * yOr3;
+
+ // upper-right part is symmetric
+ rHessian[0][1] = rHessian[1][0];
+ rHessian[0][2] = rHessian[2][0];
+ rHessian[1][2] = rHessian[2][1];
+
+ // lower-left part of Hessian of azimuthal angle theta
+ thetaHessian = new double[2][2];
+ thetaHessian[0][0] = 2 * xOrho2 * yOrho2;
+ thetaHessian[1][0] = yOrho2 * yOrho2 - xOrho2 * xOrho2;
+ thetaHessian[1][1] = -2 * xOrho2 * yOrho2;
+
+ // upper-right part is symmetric
+ thetaHessian[0][1] = thetaHessian[1][0];
+
+ // lower-left part of Hessian of polar (co-latitude) angle phi
+ final double rhor2 = rho * r2;
+ final double rho2r2 = rho * rhor2;
+ final double rhor4 = rhor2 * r2;
+ final double rho3r4 = rhor4 * rho2;
+ final double r2P2rho2 = 3 * rho2 + z2;
+ phiHessian = new double[3][3];
+ phiHessian[0][0] = z * (rho2r2 - x2 * r2P2rho2) / rho3r4;
+ phiHessian[1][0] = -x * y * z * r2P2rho2 / rho3r4;
+ phiHessian[2][0] = x * (rho2 - z2) / rhor4;
+ phiHessian[1][1] = z * (rho2r2 - y2 * r2P2rho2) / rho3r4;
+ phiHessian[2][1] = y * (rho2 - z2) / rhor4;
+ phiHessian[2][2] = 2 * rho * zOr3 / r;
+
+ // upper-right part is symmetric
+ phiHessian[0][1] = phiHessian[1][0];
+ phiHessian[0][2] = phiHessian[2][0];
+ phiHessian[1][2] = phiHessian[2][1];
+
+ }
+
+ }
+
+ /**
+ * Replace the instance with a data transfer object for serialization.
+ * @return data transfer object that will be serialized
+ */
+ private Object writeReplace() {
+ return new DataTransferObject(v.getX(), v.getY(), v.getZ());
+ }
+
+ /** Internal class used only for serialization. */
+ private static class DataTransferObject implements Serializable {
+
+ /** Serializable UID. */
+ private static final long serialVersionUID = 20130206L;
+
+ /** Abscissa.
+ * @serial
+ */
+ private final double x;
+
+ /** Ordinate.
+ * @serial
+ */
+ private final double y;
+
+ /** Height.
+ * @serial
+ */
+ private final double z;
+
+ /** Simple constructor.
+ * @param x abscissa
+ * @param y ordinate
+ * @param z height
+ */
+ public DataTransferObject(final double x, final double y, final double
z) {
+ this.x = x;
+ this.y = y;
+ this.z = z;
+ }
+
+ /** Replace the deserialized data transfer object with a {@link
SphericalCoordinates}.
+ * @return replacement {@link SphericalCoordinates}
+ */
+ private Object readResolve() {
+ return new SphericalCoordinates(new Vector3D(x, y, z));
+ }
+
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinates.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java?rev=1443178&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java
Wed Feb 6 20:08:33 2013
@@ -0,0 +1,186 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.geometry.euclidean.threed;
+
+import org.apache.commons.math3.TestUtils;
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class SphericalCoordinatesTest {
+
+ @Test
+ public void testCoordinatesStoC() throws DimensionMismatchException {
+ double piO2 = 0.5 * FastMath.PI;
+ SphericalCoordinates sc1 = new SphericalCoordinates(2.0, 0, piO2);
+ Assert.assertEquals(0, sc1.getCartesian().distance(new Vector3D(2, 0,
0)), 1.0e-10);
+ SphericalCoordinates sc2 = new SphericalCoordinates(2.0, piO2, piO2);
+ Assert.assertEquals(0, sc2.getCartesian().distance(new Vector3D(0, 2,
0)), 1.0e-10);
+ SphericalCoordinates sc3 = new SphericalCoordinates(2.0, FastMath.PI,
piO2);
+ Assert.assertEquals(0, sc3.getCartesian().distance(new Vector3D(-2, 0,
0)), 1.0e-10);
+ SphericalCoordinates sc4 = new SphericalCoordinates(2.0, -piO2, piO2);
+ Assert.assertEquals(0, sc4.getCartesian().distance(new Vector3D(0, -2,
0)), 1.0e-10);
+ SphericalCoordinates sc5 = new SphericalCoordinates(2.0, 1.23456, 0);
+ Assert.assertEquals(0, sc5.getCartesian().distance(new Vector3D(0, 0,
2)), 1.0e-10);
+ SphericalCoordinates sc6 = new SphericalCoordinates(2.0, 6.54321,
FastMath.PI);
+ Assert.assertEquals(0, sc6.getCartesian().distance(new Vector3D(0, 0,
-2)), 1.0e-10);
+ }
+
+ @Test
+ public void testCoordinatesCtoS() throws DimensionMismatchException {
+ double piO2 = 0.5 * FastMath.PI;
+ SphericalCoordinates sc1 = new SphericalCoordinates(new Vector3D(2, 0,
0));
+ Assert.assertEquals(2, sc1.getR(), 1.0e-10);
+ Assert.assertEquals(0, sc1.getTheta(), 1.0e-10);
+ Assert.assertEquals(piO2, sc1.getPhi(), 1.0e-10);
+ SphericalCoordinates sc2 = new SphericalCoordinates(new Vector3D(0, 2,
0));
+ Assert.assertEquals(2, sc2.getR(), 1.0e-10);
+ Assert.assertEquals(piO2, sc2.getTheta(), 1.0e-10);
+ Assert.assertEquals(piO2, sc2.getPhi(), 1.0e-10);
+ SphericalCoordinates sc3 = new SphericalCoordinates(new Vector3D(-2,
0, 0));
+ Assert.assertEquals(2, sc3.getR(), 1.0e-10);
+ Assert.assertEquals(FastMath.PI, sc3.getTheta(), 1.0e-10);
+ Assert.assertEquals(piO2, sc3.getPhi(), 1.0e-10);
+ SphericalCoordinates sc4 = new SphericalCoordinates(new Vector3D(0,
-2, 0));
+ Assert.assertEquals(2, sc4.getR(), 1.0e-10);
+ Assert.assertEquals(-piO2, sc4.getTheta(), 1.0e-10);
+ Assert.assertEquals(piO2, sc4.getPhi(), 1.0e-10);
+ SphericalCoordinates sc5 = new SphericalCoordinates(new Vector3D(0, 0,
2));
+ Assert.assertEquals(2, sc5.getR(), 1.0e-10);
+ // don't check theta on poles, as it is singular
+ Assert.assertEquals(0, sc5.getPhi(), 1.0e-10);
+ SphericalCoordinates sc6 = new SphericalCoordinates(new Vector3D(0, 0,
-2));
+ Assert.assertEquals(2, sc6.getR(), 1.0e-10);
+ // don't check theta on poles, as it is singular
+ Assert.assertEquals(FastMath.PI, sc6.getPhi(), 1.0e-10);
+ }
+
+ @Test
+ public void testGradient() {
+ for (double r = 0.2; r < 10; r += 0.5) {
+ for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
+ for (double phi = 0.1; phi < FastMath.PI; phi += 0.1) {
+ SphericalCoordinates sc = new SphericalCoordinates(r,
theta, phi);
+
+ DerivativeStructure svalue = valueSpherical(new
DerivativeStructure(3, 1, 0, r),
+ new
DerivativeStructure(3, 1, 1, theta),
+ new
DerivativeStructure(3, 1, 2, phi));
+ double[] sGradient = new double[] {
+ svalue.getPartialDerivative(1, 0, 0),
+ svalue.getPartialDerivative(0, 1, 0),
+ svalue.getPartialDerivative(0, 0, 1),
+ };
+
+ DerivativeStructure cvalue = valueCartesian(new
DerivativeStructure(3, 1, 0, sc.getCartesian().getX()),
+ new
DerivativeStructure(3, 1, 1, sc.getCartesian().getY()),
+ new
DerivativeStructure(3, 1, 2, sc.getCartesian().getZ()));
+ Vector3D refCGradient = new
Vector3D(cvalue.getPartialDerivative(1, 0, 0),
+
cvalue.getPartialDerivative(0, 1, 0),
+
cvalue.getPartialDerivative(0, 0, 1));
+
+ Vector3D testCGradient = new
Vector3D(sc.toCartesianGradient(sGradient));
+
+ Assert.assertEquals(0,
testCGradient.distance(refCGradient) / refCGradient.getNorm(), 5.0e-14);
+
+ }
+ }
+ }
+ }
+
+ @Test
+ public void testHessian() {
+ for (double r = 0.2; r < 10; r += 0.5) {
+ for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.2) {
+ for (double phi = 0.1; phi < FastMath.PI; phi += 0.2) {
+ SphericalCoordinates sc = new SphericalCoordinates(r,
theta, phi);
+
+ DerivativeStructure svalue = valueSpherical(new
DerivativeStructure(3, 2, 0, r),
+ new
DerivativeStructure(3, 2, 1, theta),
+ new
DerivativeStructure(3, 2, 2, phi));
+ double[] sGradient = new double[] {
+ svalue.getPartialDerivative(1, 0, 0),
+ svalue.getPartialDerivative(0, 1, 0),
+ svalue.getPartialDerivative(0, 0, 1),
+ };
+ double[][] sHessian = new double[3][3];
+ sHessian[0][0] = svalue.getPartialDerivative(2, 0, 0); //
d2F/dR2
+ sHessian[1][0] = svalue.getPartialDerivative(1, 1, 0); //
d2F/dRdTheta
+ sHessian[2][0] = svalue.getPartialDerivative(1, 0, 1); //
d2F/dRdPhi
+ sHessian[0][1] = Double.NaN; // just to check upper-right
part is not used
+ sHessian[1][1] = svalue.getPartialDerivative(0, 2, 0); //
d2F/dTheta2
+ sHessian[2][1] = svalue.getPartialDerivative(0, 1, 1); //
d2F/dThetadPhi
+ sHessian[0][2] = Double.NaN; // just to check upper-right
part is not used
+ sHessian[1][2] = Double.NaN; // just to check upper-right
part is not used
+ sHessian[2][2] = svalue.getPartialDerivative(0, 0, 2); //
d2F/dPhi2
+
+ DerivativeStructure cvalue = valueCartesian(new
DerivativeStructure(3, 2, 0, sc.getCartesian().getX()),
+ new
DerivativeStructure(3, 2, 1, sc.getCartesian().getY()),
+ new
DerivativeStructure(3, 2, 2, sc.getCartesian().getZ()));
+ double[][] refCHessian = new double[3][3];
+ refCHessian[0][0] = cvalue.getPartialDerivative(2, 0, 0);
// d2F/dX2
+ refCHessian[1][0] = cvalue.getPartialDerivative(1, 1, 0);
// d2F/dXdY
+ refCHessian[2][0] = cvalue.getPartialDerivative(1, 0, 1);
// d2F/dXdZ
+ refCHessian[0][1] = refCHessian[1][0];
+ refCHessian[1][1] = cvalue.getPartialDerivative(0, 2, 0);
// d2F/dY2
+ refCHessian[2][1] = cvalue.getPartialDerivative(0, 1, 1);
// d2F/dYdZ
+ refCHessian[0][2] = refCHessian[2][0];
+ refCHessian[1][2] = refCHessian[2][1];
+ refCHessian[2][2] = cvalue.getPartialDerivative(0, 0, 2);
// d2F/dZ2
+ double norm = 0;
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ norm = FastMath.max(norm,
FastMath.abs(refCHessian[i][j]));
+ }
+ }
+
+ double[][] testCHessian = sc.toCartesianHessian(sHessian,
sGradient);
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ Assert.assertEquals("" +
FastMath.abs((refCHessian[i][j] - testCHessian[i][j]) / norm),
+ refCHessian[i][j],
testCHessian[i][j], 1.0e-14 * norm);
+ }
+ }
+
+ }
+ }
+ }
+ }
+
+ public DerivativeStructure valueCartesian(DerivativeStructure x,
DerivativeStructure y, DerivativeStructure z) {
+ return x.divide(y.multiply(5).add(10)).multiply(z.pow(3));
+ }
+
+ public DerivativeStructure valueSpherical(DerivativeStructure r,
DerivativeStructure theta, DerivativeStructure phi) {
+ return valueCartesian(r.multiply(theta.cos()).multiply(phi.sin()),
+ r.multiply(theta.sin()).multiply(phi.sin()),
+ r.multiply(phi.cos()));
+ }
+
+ @Test
+ public void testSerialization() {
+ SphericalCoordinates a = new SphericalCoordinates(3, 2, 1);
+ SphericalCoordinates b = (SphericalCoordinates)
TestUtils.serializeAndRecover(a);
+ Assert.assertEquals(0, a.getCartesian().distance(b.getCartesian()),
1.0e-10);
+ Assert.assertEquals(a.getR(), b.getR(), 1.0e-10);
+ Assert.assertEquals(a.getTheta(), b.getTheta(), 1.0e-10);
+ Assert.assertEquals(a.getPhi(), b.getPhi(), 1.0e-10);
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/geometry/euclidean/threed/SphericalCoordinatesTest.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
