SphericalCoordinates.java

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

  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. package org.hipparchus.geometry.euclidean.threed;


  25. import java.io.Serializable;

  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.SinCos;

  30. /** This class provides conversions related to <a
  31.  * href="http://mathworld.wolfram.com/SphericalCoordinates.html">spherical coordinates</a>.
  32.  * <p>
  33.  * The conventions used here are the mathematical ones, i.e. spherical coordinates are
  34.  * related to Cartesian coordinates as follows:
  35.  * </p>
  36.  * <ul>
  37.  *   <li>x = r cos(&theta;) sin(&Phi;)</li>
  38.  *   <li>y = r sin(&theta;) sin(&Phi;)</li>
  39.  *   <li>z = r cos(&Phi;)</li>
  40.  * </ul>
  41.  * <ul>
  42.  *   <li>r       = &radic;(x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>)</li>
  43.  *   <li>&theta; = atan2(y, x)</li>
  44.  *   <li>&Phi;   = acos(z/r)</li>
  45.  * </ul>
  46.  * <p>
  47.  * r is the radius, &theta; is the azimuthal angle in the x-y plane and &Phi; is the polar
  48.  * (co-latitude) angle. These conventions are <em>different</em> from the conventions used
  49.  * in physics (and in particular in spherical harmonics) where the meanings of &theta; and
  50.  * &Phi; are reversed.
  51.  * </p>
  52.  * <p>
  53.  * This class provides conversion of coordinates and also of gradient and Hessian
  54.  * between spherical and Cartesian coordinates.
  55.  * </p>
  56.  */
  57. public class SphericalCoordinates implements Serializable {

  59.     /** Serializable UID. */
  60.     private static final long serialVersionUID = 20130206L;

  62.     /** Cartesian coordinates. */
  63.     private final Vector3D v;

  65.     /** Radius. */
  66.     private final double r;

  68.     /** Azimuthal angle in the x-y plane &theta;. */
  69.     private final double theta;

  71.     /** Polar angle (co-latitude) &Phi;. */
  72.     private final double phi;

  74.     /** Jacobian of (r, &theta; &Phi;). */
  75.     private double[][] jacobian;

  77.     /** Hessian of radius. */
  78.     private double[][] rHessian;

  80.     /** Hessian of azimuthal angle in the x-y plane &theta;. */
  81.     private double[][] thetaHessian;

  83.     /** Hessian of polar (co-latitude) angle &Phi;. */
  84.     private double[][] phiHessian;

  86.     /** Build a spherical coordinates transformer from Cartesian coordinates.
  87.      * @param v Cartesian coordinates
  88.      */
  89.     public SphericalCoordinates(final Vector3D v) {

  91.         // Cartesian coordinates
  92.         this.v = v;

  94.         // remaining spherical coordinates
  95.         this.r     = v.getNorm();
  96.         this.theta = v.getAlpha();
  97.         this.phi   = FastMath.acos(v.getZ() / r);

  99.     }

  101.     /** Build a spherical coordinates transformer from spherical coordinates.
  102.      * @param r radius
  103.      * @param theta azimuthal angle in x-y plane
  104.      * @param phi polar (co-latitude) angle
  105.      */
  106.     public SphericalCoordinates(final double r, final double theta, final double phi) {

  108.         final SinCos sinCosTheta = FastMath.sinCos(theta);
  109.         final SinCos sinCosPhi   = FastMath.sinCos(phi);

  111.         // spherical coordinates
  112.         this.r     = r;
  113.         this.theta = theta;
  114.         this.phi   = phi;

  116.         // Cartesian coordinates
  117.         this.v  = new Vector3D(r * sinCosTheta.cos() * sinCosPhi.sin(),
  118.                                r * sinCosTheta.sin() * sinCosPhi.sin(),
  119.                                r * sinCosPhi.cos());

  121.     }

  123.     /** Get the Cartesian coordinates.
  124.      * @return Cartesian coordinates
  125.      */
  126.     public Vector3D getCartesian() {
  127.         return v;
  128.     }

  130.     /** Get the radius.
  131.      * @return radius r
  132.      * @see #getTheta()
  133.      * @see #getPhi()
  134.      */
  135.     public double getR() {
  136.         return r;
  137.     }

  139.     /** Get the azimuthal angle in x-y plane.
  140.      * @return azimuthal angle in x-y plane &theta;
  141.      * @see #getR()
  142.      * @see #getPhi()
  143.      */
  144.     public double getTheta() {
  145.         return theta;
  146.     }

  148.     /** Get the polar (co-latitude) angle.
  149.      * @return polar (co-latitude) angle &Phi;
  150.      * @see #getR()
  151.      * @see #getTheta()
  152.      */
  153.     public double getPhi() {
  154.         return phi;
  155.     }

  157.     /** Convert a gradient with respect to spherical coordinates into a gradient
  158.      * with respect to Cartesian coordinates.
  159.      * @param sGradient gradient with respect to spherical coordinates
  160.      * {df/dr, df/d&theta;, df/d&Phi;}
  161.      * @return gradient with respect to Cartesian coordinates
  162.      * {df/dx, df/dy, df/dz}
  163.      */
  164.     public double[] toCartesianGradient(final double[] sGradient) {

  166.         // lazy evaluation of Jacobian
  167.         computeJacobian();

  169.         // compose derivatives as gradient^T . J
  170.         // the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0
  171.         return new double[] {
  172.             sGradient[0] * jacobian[0][0] + sGradient[1] * jacobian[1][0] + sGradient[2] * jacobian[2][0],
  173.             sGradient[0] * jacobian[0][1] + sGradient[1] * jacobian[1][1] + sGradient[2] * jacobian[2][1],
  174.             sGradient[0] * jacobian[0][2]                                 + sGradient[2] * jacobian[2][2]
  175.         };

  177.     }

  179.     /** Convert a Hessian with respect to spherical coordinates into a Hessian
  180.      * with respect to Cartesian coordinates.
  181.      * <p>
  182.      * As Hessian are always symmetric, we use only the lower left part of the provided
  183.      * spherical Hessian, so the upper part may not be initialized. However, we still
  184.      * do fill up the complete array we create, with guaranteed symmetry.
  185.      * </p>
  186.      * @param sHessian Hessian with respect to spherical coordinates
  187.      * {{d<sup>2</sup>f/dr<sup>2</sup>, d<sup>2</sup>f/drd&theta;, d<sup>2</sup>f/drd&Phi;},
  188.      *  {d<sup>2</sup>f/drd&theta;, d<sup>2</sup>f/d&theta;<sup>2</sup>, d<sup>2</sup>f/d&theta;d&Phi;},
  189.      *  {d<sup>2</sup>f/drd&Phi;, d<sup>2</sup>f/d&theta;d&Phi;, d<sup>2</sup>f/d&Phi;<sup>2</sup>}
  190.      * @param sGradient gradient with respect to spherical coordinates
  191.      * {df/dr, df/d&theta;, df/d&Phi;}
  192.      * @return Hessian with respect to Cartesian coordinates
  193.      * {{d<sup>2</sup>f/dx<sup>2</sup>, d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dxdz},
  194.      *  {d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dy<sup>2</sup>, d<sup>2</sup>f/dydz},
  195.      *  {d<sup>2</sup>f/dxdz, d<sup>2</sup>f/dydz, d<sup>2</sup>f/dz<sup>2</sup>}}
  196.      */
  197.     public double[][] toCartesianHessian(final double[][] sHessian, final double[] sGradient) {

  199.         computeJacobian();
  200.         computeHessians();

  202.         // compose derivative as J^T . H_f . J + df/dr H_r + df/dtheta H_theta + df/dphi H_phi
  203.         // the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0
  204.         // and H_theta is only a 2x2 matrix as it does not depend on z
  205.         final double[][] hj = new double[3][3];
  206.         final double[][] cHessian = new double[3][3];

  208.         // compute H_f . J
  209.         // beware we use ONLY the lower-left part of sHessian
  210.         hj[0][0] = sHessian[0][0] * jacobian[0][0] + sHessian[1][0] * jacobian[1][0] + sHessian[2][0] * jacobian[2][0];
  211.         hj[0][1] = sHessian[0][0] * jacobian[0][1] + sHessian[1][0] * jacobian[1][1] + sHessian[2][0] * jacobian[2][1];
  212.         hj[0][2] = sHessian[0][0] * jacobian[0][2]                                   + sHessian[2][0] * jacobian[2][2];
  213.         hj[1][0] = sHessian[1][0] * jacobian[0][0] + sHessian[1][1] * jacobian[1][0] + sHessian[2][1] * jacobian[2][0];
  214.         hj[1][1] = sHessian[1][0] * jacobian[0][1] + sHessian[1][1] * jacobian[1][1] + sHessian[2][1] * jacobian[2][1];
  215.         // don't compute hj[1][2] as it is not used below
  216.         hj[2][0] = sHessian[2][0] * jacobian[0][0] + sHessian[2][1] * jacobian[1][0] + sHessian[2][2] * jacobian[2][0];
  217.         hj[2][1] = sHessian[2][0] * jacobian[0][1] + sHessian[2][1] * jacobian[1][1] + sHessian[2][2] * jacobian[2][1];
  218.         hj[2][2] = sHessian[2][0] * jacobian[0][2]                                   + sHessian[2][2] * jacobian[2][2];

  220.         // compute lower-left part of J^T . H_f . J
  221.         cHessian[0][0] = jacobian[0][0] * hj[0][0] + jacobian[1][0] * hj[1][0] + jacobian[2][0] * hj[2][0];
  222.         cHessian[1][0] = jacobian[0][1] * hj[0][0] + jacobian[1][1] * hj[1][0] + jacobian[2][1] * hj[2][0];
  223.         cHessian[2][0] = jacobian[0][2] * hj[0][0]                             + jacobian[2][2] * hj[2][0];
  224.         cHessian[1][1] = jacobian[0][1] * hj[0][1] + jacobian[1][1] * hj[1][1] + jacobian[2][1] * hj[2][1];
  225.         cHessian[2][1] = jacobian[0][2] * hj[0][1]                             + jacobian[2][2] * hj[2][1];
  226.         cHessian[2][2] = jacobian[0][2] * hj[0][2]                             + jacobian[2][2] * hj[2][2];

  228.         // add gradient contribution
  229.         cHessian[0][0] += sGradient[0] * rHessian[0][0] + sGradient[1] * thetaHessian[0][0] + sGradient[2] * phiHessian[0][0];
  230.         cHessian[1][0] += sGradient[0] * rHessian[1][0] + sGradient[1] * thetaHessian[1][0] + sGradient[2] * phiHessian[1][0];
  231.         cHessian[2][0] += sGradient[0] * rHessian[2][0]                                     + sGradient[2] * phiHessian[2][0];
  232.         cHessian[1][1] += sGradient[0] * rHessian[1][1] + sGradient[1] * thetaHessian[1][1] + sGradient[2] * phiHessian[1][1];
  233.         cHessian[2][1] += sGradient[0] * rHessian[2][1]                                     + sGradient[2] * phiHessian[2][1];
  234.         cHessian[2][2] += sGradient[0] * rHessian[2][2]                                     + sGradient[2] * phiHessian[2][2];

  236.         // ensure symmetry
  237.         cHessian[0][1] = cHessian[1][0];
  238.         cHessian[0][2] = cHessian[2][0];
  239.         cHessian[1][2] = cHessian[2][1];

  241.         return cHessian;

  243.     }

  245.     /** Lazy evaluation of (r, &theta;, &phi;) Jacobian.
  246.      */
  247.     private void computeJacobian() {
  248.         if (jacobian == null) {

  250.             // intermediate variables
  251.             final double x    = v.getX();
  252.             final double y    = v.getY();
  253.             final double z    = v.getZ();
  254.             final double rho2 = x * x + y * y;
  255.             final double rho  = FastMath.sqrt(rho2);
  256.             final double r2   = rho2 + z * z;

  258.             jacobian = new double[3][3];

  260.             // row representing the gradient of r
  261.             jacobian[0][0] = x / r;
  262.             jacobian[0][1] = y / r;
  263.             jacobian[0][2] = z / r;

  265.             // row representing the gradient of theta
  266.             jacobian[1][0] = -y / rho2;
  267.             jacobian[1][1] =  x / rho2;
  268.             // jacobian[1][2] is already set to 0 at allocation time

  270.             // row representing the gradient of phi
  271.             jacobian[2][0] = x * z / (rho * r2);
  272.             jacobian[2][1] = y * z / (rho * r2);
  273.             jacobian[2][2] = -rho / r2;

  275.         }
  276.     }

  278.     /** Lazy evaluation of Hessians.
  279.      */
  280.     private void computeHessians() {

  282.         if (rHessian == null) {

  284.             // intermediate variables
  285.             final double x      = v.getX();
  286.             final double y      = v.getY();
  287.             final double z      = v.getZ();
  288.             final double x2     = x * x;
  289.             final double y2     = y * y;
  290.             final double z2     = z * z;
  291.             final double rho2   = x2 + y2;
  292.             final double rho    = FastMath.sqrt(rho2);
  293.             final double r2     = rho2 + z2;
  294.             final double xOr    = x / r;
  295.             final double yOr    = y / r;
  296.             final double zOr    = z / r;
  297.             final double xOrho2 = x / rho2;
  298.             final double yOrho2 = y / rho2;
  299.             final double xOr3   = xOr / r2;
  300.             final double yOr3   = yOr / r2;
  301.             final double zOr3   = zOr / r2;

  303.             // lower-left part of Hessian of r
  304.             rHessian = new double[3][3];
  305.             rHessian[0][0] = y * yOr3 + z * zOr3;
  306.             rHessian[1][0] = -x * yOr3;
  307.             rHessian[2][0] = -z * xOr3;
  308.             rHessian[1][1] = x * xOr3 + z * zOr3;
  309.             rHessian[2][1] = -y * zOr3;
  310.             rHessian[2][2] = x * xOr3 + y * yOr3;

  312.             // upper-right part is symmetric
  313.             rHessian[0][1] = rHessian[1][0];
  314.             rHessian[0][2] = rHessian[2][0];
  315.             rHessian[1][2] = rHessian[2][1];

  317.             // lower-left part of Hessian of azimuthal angle theta
  318.             thetaHessian = new double[2][2];
  319.             thetaHessian[0][0] = 2 * xOrho2 * yOrho2;
  320.             thetaHessian[1][0] = yOrho2 * yOrho2 - xOrho2 * xOrho2;
  321.             thetaHessian[1][1] = -2 * xOrho2 * yOrho2;

  323.             // upper-right part is symmetric
  324.             thetaHessian[0][1] = thetaHessian[1][0];

  326.             // lower-left part of Hessian of polar (co-latitude) angle phi
  327.             final double rhor2       = rho * r2;
  328.             final double rho2r2      = rho * rhor2;
  329.             final double rhor4       = rhor2 * r2;
  330.             final double rho3r4      = rhor4 * rho2;
  331.             final double r2P2rho2    = 3 * rho2 + z2;
  332.             phiHessian = new double[3][3];
  333.             phiHessian[0][0] = z * (rho2r2 - x2 * r2P2rho2) / rho3r4;
  334.             phiHessian[1][0] = -x * y * z * r2P2rho2 / rho3r4;
  335.             phiHessian[2][0] = x * (rho2 - z2) / rhor4;
  336.             phiHessian[1][1] = z * (rho2r2 - y2 * r2P2rho2) / rho3r4;
  337.             phiHessian[2][1] = y * (rho2 - z2) / rhor4;
  338.             phiHessian[2][2] = 2 * rho * zOr3 / r;

  340.             // upper-right part is symmetric
  341.             phiHessian[0][1] = phiHessian[1][0];
  342.             phiHessian[0][2] = phiHessian[2][0];
  343.             phiHessian[1][2] = phiHessian[2][1];

  345.         }

  347.     }

  349.     /**
  350.      * Replace the instance with a data transfer object for serialization.
  351.      * @return data transfer object that will be serialized
  352.      */
  353.     private Object writeReplace() {
  354.         return new DataTransferObject(v.getX(), v.getY(), v.getZ());
  355.     }

  357.     /** Internal class used only for serialization. */
  358.     private static class DataTransferObject implements Serializable {

  360.         /** Serializable UID. */
  361.         private static final long serialVersionUID = 20130206L;

  363.         /** Abscissa.
  364.          * @serial
  365.          */
  366.         private final double x;

  368.         /** Ordinate.
  369.          * @serial
  370.          */
  371.         private final double y;

  373.         /** Height.
  374.          * @serial
  375.          */
  376.         private final double z;

  378.         /** Simple constructor.
  379.          * @param x abscissa
  380.          * @param y ordinate
  381.          * @param z height
  382.          */
  383.         DataTransferObject(final double x, final double y, final double z) {
  384.             this.x = x;
  385.             this.y = y;
  386.             this.z = z;
  387.         }

  389.         /** Replace the deserialized data transfer object with a {@link SphericalCoordinates}.
  390.          * @return replacement {@link SphericalCoordinates}
  391.          */
  392.         private Object readResolve() {
  393.             return new SphericalCoordinates(new Vector3D(x, y, z));
  394.         }

  396.     }

  398. }
Created with JaCoCo 0.8.12.202403310830