View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.geometry.euclidean.threed;
24  
25  import org.hipparchus.UnitTestUtils;
26  import org.hipparchus.analysis.differentiation.DSFactory;
27  import org.hipparchus.analysis.differentiation.DerivativeStructure;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.util.FastMath;
30  import org.hipparchus.util.MathUtils;
31  import org.junit.jupiter.api.Test;
32  
33  import static org.junit.jupiter.api.Assertions.assertEquals;
34  
35  public class SphericalCoordinatesTest {
36  
37      @Test
38      void testCoordinatesStoC() throws MathIllegalArgumentException {
39          SphericalCoordinates sc1 = new SphericalCoordinates(2.0, 0, MathUtils.SEMI_PI);
40          assertEquals(0, sc1.getCartesian().distance(new Vector3D(2, 0, 0)), 1.0e-10);
41          SphericalCoordinates sc2 = new SphericalCoordinates(2.0, MathUtils.SEMI_PI, MathUtils.SEMI_PI);
42          assertEquals(0, sc2.getCartesian().distance(new Vector3D(0, 2, 0)), 1.0e-10);
43          SphericalCoordinates sc3 = new SphericalCoordinates(2.0, FastMath.PI, MathUtils.SEMI_PI);
44          assertEquals(0, sc3.getCartesian().distance(new Vector3D(-2, 0, 0)), 1.0e-10);
45          SphericalCoordinates sc4 = new SphericalCoordinates(2.0, -MathUtils.SEMI_PI, MathUtils.SEMI_PI);
46          assertEquals(0, sc4.getCartesian().distance(new Vector3D(0, -2, 0)), 1.0e-10);
47          SphericalCoordinates sc5 = new SphericalCoordinates(2.0, 1.23456, 0);
48          assertEquals(0, sc5.getCartesian().distance(new Vector3D(0, 0, 2)), 1.0e-10);
49          SphericalCoordinates sc6 = new SphericalCoordinates(2.0, 6.54321, FastMath.PI);
50          assertEquals(0, sc6.getCartesian().distance(new Vector3D(0, 0, -2)), 1.0e-10);
51      }
52  
53      @Test
54      void testCoordinatesCtoS() throws MathIllegalArgumentException {
55          double piO2 = MathUtils.SEMI_PI;
56          SphericalCoordinates sc1 = new SphericalCoordinates(new Vector3D(2, 0, 0));
57          assertEquals(2,           sc1.getR(),     1.0e-10);
58          assertEquals(0,           sc1.getTheta(), 1.0e-10);
59          assertEquals(piO2,        sc1.getPhi(),   1.0e-10);
60          SphericalCoordinates sc2 = new SphericalCoordinates(new Vector3D(0, 2, 0));
61          assertEquals(2,           sc2.getR(),     1.0e-10);
62          assertEquals(piO2,        sc2.getTheta(), 1.0e-10);
63          assertEquals(piO2,        sc2.getPhi(),   1.0e-10);
64          SphericalCoordinates sc3 = new SphericalCoordinates(new Vector3D(-2, 0, 0));
65          assertEquals(2,           sc3.getR(),     1.0e-10);
66          assertEquals(FastMath.PI, sc3.getTheta(), 1.0e-10);
67          assertEquals(piO2,        sc3.getPhi(),   1.0e-10);
68          SphericalCoordinates sc4 = new SphericalCoordinates(new Vector3D(0, -2, 0));
69          assertEquals(2,           sc4.getR(),     1.0e-10);
70          assertEquals(-piO2,       sc4.getTheta(), 1.0e-10);
71          assertEquals(piO2,        sc4.getPhi(),   1.0e-10);
72          SphericalCoordinates sc5 = new SphericalCoordinates(new Vector3D(0, 0, 2));
73          assertEquals(2,           sc5.getR(),     1.0e-10);
74          //  don't check theta on poles, as it is singular
75          assertEquals(0,           sc5.getPhi(),   1.0e-10);
76          SphericalCoordinates sc6 = new SphericalCoordinates(new Vector3D(0, 0, -2));
77          assertEquals(2,           sc6.getR(),     1.0e-10);
78          //  don't check theta on poles, as it is singular
79          assertEquals(FastMath.PI, sc6.getPhi(),   1.0e-10);
80      }
81  
82      @Test
83      void testGradient() {
84          DSFactory factory = new DSFactory(3, 1);
85          for (double r = 0.2; r < 10; r += 0.5) {
86              for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
87                  for (double phi = 0.1; phi < FastMath.PI; phi += 0.1) {
88                      SphericalCoordinates sc = new SphericalCoordinates(r, theta, phi);
89  
90                      DerivativeStructure svalue = valueSpherical(factory.variable(0, r),
91                                                                  factory.variable(1, theta),
92                                                                  factory.variable(2, phi));
93                      double[] sGradient = new double[] {
94                          svalue.getPartialDerivative(1, 0, 0),
95                          svalue.getPartialDerivative(0, 1, 0),
96                          svalue.getPartialDerivative(0, 0, 1),
97                      };
98  
99                      DerivativeStructure cvalue = valueCartesian(factory.variable(0, sc.getCartesian().getX()),
100                                                                 factory.variable(1, sc.getCartesian().getY()),
101                                                                 factory.variable(2, sc.getCartesian().getZ()));
102                     Vector3D refCGradient = new Vector3D(cvalue.getPartialDerivative(1, 0, 0),
103                                                          cvalue.getPartialDerivative(0, 1, 0),
104                                                          cvalue.getPartialDerivative(0, 0, 1));
105 
106                     Vector3D testCGradient = new Vector3D(sc.toCartesianGradient(sGradient));
107 
108                     assertEquals(0, testCGradient.distance(refCGradient) / refCGradient.getNorm(), 5.0e-14);
109 
110                 }
111             }
112         }
113     }
114 
115     @Test
116     void testHessian() {
117         DSFactory factory = new DSFactory(3, 2);
118         for (double r = 0.2; r < 10; r += 0.5) {
119             for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.2) {
120                 for (double phi = 0.1; phi < FastMath.PI; phi += 0.2) {
121                     SphericalCoordinates sc = new SphericalCoordinates(r, theta, phi);
122 
123                     DerivativeStructure svalue = valueSpherical(factory.variable(0, r),
124                                                                 factory.variable(1, theta),
125                                                                 factory.variable(2, phi));
126                     double[] sGradient = new double[] {
127                         svalue.getPartialDerivative(1, 0, 0),
128                         svalue.getPartialDerivative(0, 1, 0),
129                         svalue.getPartialDerivative(0, 0, 1),
130                     };
131                     double[][] sHessian = new double[3][3];
132                     sHessian[0][0] = svalue.getPartialDerivative(2, 0, 0); // d2F/dR2
133                     sHessian[1][0] = svalue.getPartialDerivative(1, 1, 0); // d2F/dRdTheta
134                     sHessian[2][0] = svalue.getPartialDerivative(1, 0, 1); // d2F/dRdPhi
135                     sHessian[0][1] = Double.NaN; // just to check upper-right part is not used
136                     sHessian[1][1] = svalue.getPartialDerivative(0, 2, 0); // d2F/dTheta2
137                     sHessian[2][1] = svalue.getPartialDerivative(0, 1, 1); // d2F/dThetadPhi
138                     sHessian[0][2] = Double.NaN; // just to check upper-right part is not used
139                     sHessian[1][2] = Double.NaN; // just to check upper-right part is not used
140                     sHessian[2][2] = svalue.getPartialDerivative(0, 0, 2); // d2F/dPhi2
141 
142                     DerivativeStructure cvalue = valueCartesian(factory.variable(0, sc.getCartesian().getX()),
143                                                                 factory.variable(1, sc.getCartesian().getY()),
144                                                                 factory.variable(2, sc.getCartesian().getZ()));
145                     double[][] refCHessian = new double[3][3];
146                     refCHessian[0][0] = cvalue.getPartialDerivative(2, 0, 0); // d2F/dX2
147                     refCHessian[1][0] = cvalue.getPartialDerivative(1, 1, 0); // d2F/dXdY
148                     refCHessian[2][0] = cvalue.getPartialDerivative(1, 0, 1); // d2F/dXdZ
149                     refCHessian[0][1] = refCHessian[1][0];
150                     refCHessian[1][1] = cvalue.getPartialDerivative(0, 2, 0); // d2F/dY2
151                     refCHessian[2][1] = cvalue.getPartialDerivative(0, 1, 1); // d2F/dYdZ
152                     refCHessian[0][2] = refCHessian[2][0];
153                     refCHessian[1][2] = refCHessian[2][1];
154                     refCHessian[2][2] = cvalue.getPartialDerivative(0, 0, 2); // d2F/dZ2
155                     double norm =  0;
156                     for (int i = 0; i < 3; ++i) {
157                         for (int j = 0; j < 3; ++j) {
158                             norm = FastMath.max(norm, FastMath.abs(refCHessian[i][j]));
159                         }
160                     }
161 
162                     double[][] testCHessian = sc.toCartesianHessian(sHessian, sGradient);
163                     for (int i = 0; i < 3; ++i) {
164                         for (int j = 0; j < 3; ++j) {
165                             assertEquals(refCHessian[i][j], testCHessian[i][j], 1.0e-14 * norm, "" + FastMath.abs((refCHessian[i][j] - testCHessian[i][j]) / norm));
166                         }
167                     }
168 
169                 }
170             }
171         }
172     }
173 
174     public DerivativeStructure valueCartesian(DerivativeStructure x, DerivativeStructure y, DerivativeStructure z) {
175         return x.divide(y.multiply(5).add(10)).multiply(z.pow(3));
176     }
177 
178     public DerivativeStructure valueSpherical(DerivativeStructure r, DerivativeStructure theta, DerivativeStructure phi) {
179         return valueCartesian(r.multiply(theta.cos()).multiply(phi.sin()),
180                               r.multiply(theta.sin()).multiply(phi.sin()),
181                               r.multiply(phi.cos()));
182     }
183 
184     @Test
185     void testSerialization() {
186         SphericalCoordinates a = new SphericalCoordinates(3, 2, 1);
187         SphericalCoordinates b = (SphericalCoordinates) UnitTestUtils.serializeAndRecover(a);
188         assertEquals(0, a.getCartesian().distance(b.getCartesian()), 1.0e-10);
189         assertEquals(a.getR(),     b.getR(),     1.0e-10);
190         assertEquals(a.getTheta(), b.getTheta(), 1.0e-10);
191         assertEquals(a.getPhi(),   b.getPhi(),   1.0e-10);
192     }
193 
194 }