Ellipsoid.java

  1. /* Copyright 2002-2024 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS 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.  *   http://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. package org.orekit.bodies;

  18. import java.io.Serializable;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.exception.MathRuntimeException;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.geometry.euclidean.twod.FieldVector2D;
  24. import org.hipparchus.geometry.euclidean.twod.Vector2D;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.MathArrays;
  27. import org.hipparchus.util.Precision;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.frames.Frame;

  31. /**
  32.  * Modeling of a general three-axes ellipsoid.
  33.  * @since 7.0
  34.  * @author Luc Maisonobe
  35.  */
  36. public class Ellipsoid implements Serializable {

  37.     /** Serializable UID. */
  38.     private static final long serialVersionUID = 20140924L;

  39.     /** Frame at the ellipsoid center, aligned with principal axes. */
  40.     private final Frame frame;

  41.     /** First semi-axis length. */
  42.     private final double a;

  43.     /** Second semi-axis length. */
  44.     private final double b;

  45.     /** Third semi-axis length. */
  46.     private final double c;

  47.     /** Simple constructor.
  48.      * @param frame at the ellipsoid center, aligned with principal axes
  49.      * @param a first semi-axis length
  50.      * @param b second semi-axis length
  51.      * @param c third semi-axis length
  52.      */
  53.     public Ellipsoid(final Frame frame, final double a, final double b, final double c) {
  54.         this.frame = frame;
  55.         this.a     = a;
  56.         this.b     = b;
  57.         this.c     = c;
  58.     }

  59.     /** Get the length of the first semi-axis.
  60.      * @return length of the first semi-axis (m)
  61.      */
  62.     public double getA() {
  63.         return a;
  64.     }

  65.     /** Get the length of the second semi-axis.
  66.      * @return length of the second semi-axis (m)
  67.      */
  68.     public double getB() {
  69.         return b;
  70.     }

  71.     /** Get the length of the third semi-axis.
  72.      * @return length of the third semi-axis (m)
  73.      */
  74.     public double getC() {
  75.         return c;
  76.     }

  77.     /** Get the ellipsoid central frame.
  78.      * @return ellipsoid central frame
  79.      */
  80.     public Frame getFrame() {
  81.         return frame;
  82.     }

  83.     /** Check if a point is inside the ellipsoid.
  84.      * @param point point to check, in the ellipsoid frame
  85.      * @return true if the point is inside the ellipsoid
  86.      * (or exactly on ellipsoid surface)
  87.      * @since 7.1
  88.      */
  89.     public boolean isInside(final Vector3D point) {
  90.         final double scaledX = point.getX() / a;
  91.         final double scaledY = point.getY() / b;
  92.         final double scaledZ = point.getZ() / c;
  93.         return scaledX * scaledX + scaledY * scaledY + scaledZ * scaledZ <= 1.0;
  94.     }

  95.     /** Check if a point is inside the ellipsoid.
  96.      * @param point point to check, in the ellipsoid frame
  97.      * @return true if the point is inside the ellipsoid
  98.      * (or exactly on ellipsoid surface)
  99.      * @param <T> the type of the field elements
  100.      * @since 12.0
  101.      */
  102.     public <T extends CalculusFieldElement<T>> boolean isInside(final FieldVector3D<T> point) {
  103.         final T scaledX = point.getX().divide(a);
  104.         final T scaledY = point.getY().divide(b);
  105.         final T scaledZ = point.getZ().divide(c);
  106.         final T d2      = scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).add(scaledZ.multiply(scaledZ));
  107.         return d2.subtract(1.0).getReal() <= 0.0;
  108.     }

  109.     /** Compute the 2D ellipse at the intersection of the 3D ellipsoid and a plane.
  110.      * @param planePoint point belonging to the plane, in the ellipsoid frame
  111.      * @param planeNormal normal of the plane, in the ellipsoid frame
  112.      * @return plane section or null if there are no intersections
  113.      * @exception MathRuntimeException if the norm of planeNormal is null
  114.      */
  115.     public Ellipse getPlaneSection(final Vector3D planePoint, final Vector3D planeNormal)
  116.         throws MathRuntimeException {

  117.         // we define the points Q in the plane using two free variables τ and υ as:
  118.         // Q = P + τ u + υ v
  119.         // where u and v are two unit vectors belonging to the plane
  120.         // Q belongs to the 3D ellipsoid so:
  121.         // (xQ / a)² + (yQ / b)² + (zQ / c)² = 1
  122.         // combining both equations, we get:
  123.         //   (xP² + 2 xP (τ xU + υ xV) + (τ xU + υ xV)²) / a²
  124.         // + (yP² + 2 yP (τ yU + υ yV) + (τ yU + υ yV)²) / b²
  125.         // + (zP² + 2 zP (τ zU + υ zV) + (τ zU + υ zV)²) / c²
  126.         // = 1
  127.         // which can be rewritten:
  128.         // α τ² + β υ² + 2 γ τυ + 2 δ τ + 2 ε υ + ζ = 0
  129.         // with
  130.         // α =  xU²  / a² +  yU²  / b² +  zU²  / c² > 0
  131.         // β =  xV²  / a² +  yV²  / b² +  zV²  / c² > 0
  132.         // γ = xU xV / a² + yU yV / b² + zU zV / c²
  133.         // δ = xP xU / a² + yP yU / b² + zP zU / c²
  134.         // ε = xP xV / a² + yP yV / b² + zP zV / c²
  135.         // ζ =  xP²  / a² +  yP²  / b² +  zP²  / c² - 1
  136.         // this is the equation of a conic (here an ellipse)
  137.         // Of course, we note that if the point P belongs to the ellipsoid
  138.         // then ζ = 0 and the equation holds at point P since τ = 0 and υ = 0
  139.         final Vector3D u     = planeNormal.orthogonal();
  140.         final Vector3D v     = Vector3D.crossProduct(planeNormal, u).normalize();
  141.         final double xUOa    = u.getX() / a;
  142.         final double yUOb    = u.getY() / b;
  143.         final double zUOc    = u.getZ() / c;
  144.         final double xVOa    = v.getX() / a;
  145.         final double yVOb    = v.getY() / b;
  146.         final double zVOc    = v.getZ() / c;
  147.         final double xPOa    = planePoint.getX() / a;
  148.         final double yPOb    = planePoint.getY() / b;
  149.         final double zPOc    = planePoint.getZ() / c;
  150.         final double alpha   = xUOa * xUOa + yUOb * yUOb + zUOc * zUOc;
  151.         final double beta    = xVOa * xVOa + yVOb * yVOb + zVOc * zVOc;
  152.         final double gamma   = MathArrays.linearCombination(xUOa, xVOa, yUOb, yVOb, zUOc, zVOc);
  153.         final double delta   = MathArrays.linearCombination(xPOa, xUOa, yPOb, yUOb, zPOc, zUOc);
  154.         final double epsilon = MathArrays.linearCombination(xPOa, xVOa, yPOb, yVOb, zPOc, zVOc);
  155.         final double zeta    = MathArrays.linearCombination(xPOa, xPOa, yPOb, yPOb, zPOc, zPOc, 1, -1);

  156.         // reduce the general equation α τ² + β υ² + 2 γ τυ + 2 δ τ + 2 ε υ + ζ = 0
  157.         // to canonical form (λ/l)² + (μ/m)² = 1
  158.         // using a coordinates change
  159.         //       τ = τC + λ cosθ - μ sinθ
  160.         //       υ = υC + λ sinθ + μ cosθ
  161.         // or equivalently
  162.         //       λ =   (τ - τC) cosθ + (υ - υC) sinθ
  163.         //       μ = - (τ - τC) sinθ + (υ - υC) cosθ
  164.         // τC and υC are the coordinates of the 2D ellipse center with respect to P
  165.         // 2l and 2m and are the axes lengths (major or minor depending on which one is greatest)
  166.         // θ is the angle of the 2D ellipse axis corresponding to axis with length 2l

  167.         // choose θ in order to cancel the coupling term in λμ
  168.         // expanding the general equation, we get: A λ² + B μ² + C λμ + D λ + E μ + F = 0
  169.         // with C = 2[(β - α) cosθ sinθ + γ (cos²θ - sin²θ)]
  170.         // hence the term is cancelled when θ = arctan(t), with γ t² + (α - β) t - γ = 0
  171.         // As the solutions of the quadratic equation obey t₁t₂ = -1, they correspond to
  172.         // angles θ in quadrature to each other. Selecting one solution or the other simply
  173.         // exchanges the principal axes. As we don't care about which axis we want as the
  174.         // first one, we select an arbitrary solution
  175.         final double tanTheta;
  176.         if (FastMath.abs(gamma) < Precision.SAFE_MIN) {
  177.             tanTheta = 0.0;
  178.         } else {
  179.             final double bMA = beta - alpha;
  180.             tanTheta = (bMA >= 0) ?
  181.                        (-2 * gamma / (bMA + FastMath.sqrt(bMA * bMA + 4 * gamma * gamma))) :
  182.                        (-2 * gamma / (bMA - FastMath.sqrt(bMA * bMA + 4 * gamma * gamma)));
  183.         }
  184.         final double tan2   = tanTheta * tanTheta;
  185.         final double cos2   = 1 / (1 + tan2);
  186.         final double sin2   = tan2 * cos2;
  187.         final double cosSin = tanTheta * cos2;
  188.         final double cos    = FastMath.sqrt(cos2);
  189.         final double sin    = tanTheta * cos;

  190.         // choose τC and υC in order to cancel the linear terms in λ and μ
  191.         // expanding the general equation, we get: A λ² + B μ² + C λμ + D λ + E μ + F = 0
  192.         // with D = 2[ (α τC + γ υC + δ) cosθ + (γ τC + β υC + ε) sinθ]
  193.         //      E = 2[-(α τC + γ υC + δ) sinθ + (γ τC + β υC + ε) cosθ]
  194.         // θ can be eliminated by combining the equations
  195.         // D cosθ - E sinθ = 2[α τC + γ υC + δ]
  196.         // E cosθ + D sinθ = 2[γ τC + β υC + ε]
  197.         // hence the terms D and E are both cancelled (regardless of θ) when
  198.         //     τC = (β δ - γ ε) / (γ² - α β)
  199.         //     υC = (α ε - γ δ) / (γ² - α β)
  200.         final double denom = MathArrays.linearCombination(gamma, gamma,   -alpha, beta);
  201.         final double tauC  = MathArrays.linearCombination(beta,  delta,   -gamma, epsilon) / denom;
  202.         final double nuC   = MathArrays.linearCombination(alpha, epsilon, -gamma, delta)   / denom;

  203.         // compute l and m
  204.         // expanding the general equation, we get: A λ² + B μ² + C λμ + D λ + E μ + F = 0
  205.         // with A = α cos²θ + β sin²θ + 2 γ cosθ sinθ
  206.         //      B = α sin²θ + β cos²θ - 2 γ cosθ sinθ
  207.         //      F = α τC² + β υC² + 2 γ τC υC + 2 δ τC + 2 ε υC + ζ
  208.         // hence we compute directly l = √(-F/A) and m = √(-F/B)
  209.         final double twogcs = 2 * gamma * cosSin;
  210.         final double bigA   = alpha * cos2 + beta * sin2 + twogcs;
  211.         final double bigB   = alpha * sin2 + beta * cos2 - twogcs;
  212.         final double bigF   = (alpha * tauC + 2 * (gamma * nuC + delta)) * tauC +
  213.                               (beta * nuC + 2 * epsilon) * nuC + zeta;
  214.         final double l      = FastMath.sqrt(-bigF / bigA);
  215.         final double m      = FastMath.sqrt(-bigF / bigB);
  216.         if (Double.isNaN(l + m)) {
  217.             // the plane does not intersect the ellipsoid
  218.             return null;
  219.         }

  220.         if (l > m) {
  221.             return new Ellipse(new Vector3D(1, planePoint, tauC, u, nuC, v),
  222.                                new Vector3D( cos, u, sin, v),
  223.                                new Vector3D(-sin, u, cos, v),
  224.                                l, m, frame);
  225.         } else {
  226.             return new Ellipse(new Vector3D(1, planePoint, tauC, u, nuC, v),
  227.                                new Vector3D(sin, u, -cos, v),
  228.                                new Vector3D(cos, u,  sin, v),
  229.                                m, l, frame);
  230.         }

  231.     }

  232.     /** Compute the 2D ellipse at the intersection of the 3D ellipsoid and a plane.
  233.      * @param planePoint point belonging to the plane, in the ellipsoid frame
  234.      * @param planeNormal normal of the plane, in the ellipsoid frame
  235.      * @return plane section or null if there are no intersections
  236.      * @exception MathRuntimeException if the norm of planeNormal is null
  237.      * @param <T> the type of the field elements
  238.      * @since 12.0
  239.      */
  240.     public <T extends CalculusFieldElement<T>> FieldEllipse<T> getPlaneSection(final FieldVector3D<T> planePoint, final FieldVector3D<T> planeNormal)
  241.         throws MathRuntimeException {

  242.         final T zero = planePoint.getX().getField().getZero();
  243.         final T one  = planePoint.getX().getField().getOne();

  244.         // we define the points Q in the plane using two free variables τ and υ as:
  245.         // Q = P + τ u + υ v
  246.         // where u and v are two unit vectors belonging to the plane
  247.         // Q belongs to the 3D ellipsoid so:
  248.         // (xQ / a)² + (yQ / b)² + (zQ / c)² = 1
  249.         // combining both equations, we get:
  250.         //   (xP² + 2 xP (τ xU + υ xV) + (τ xU + υ xV)²) / a²
  251.         // + (yP² + 2 yP (τ yU + υ yV) + (τ yU + υ yV)²) / b²
  252.         // + (zP² + 2 zP (τ zU + υ zV) + (τ zU + υ zV)²) / c²
  253.         // = 1
  254.         // which can be rewritten:
  255.         // α τ² + β υ² + 2 γ τυ + 2 δ τ + 2 ε υ + ζ = 0
  256.         // with
  257.         // α =  xU²  / a² +  yU²  / b² +  zU²  / c² > 0
  258.         // β =  xV²  / a² +  yV²  / b² +  zV²  / c² > 0
  259.         // γ = xU xV / a² + yU yV / b² + zU zV / c²
  260.         // δ = xP xU / a² + yP yU / b² + zP zU / c²
  261.         // ε = xP xV / a² + yP yV / b² + zP zV / c²
  262.         // ζ =  xP²  / a² +  yP²  / b² +  zP²  / c² - 1
  263.         // this is the equation of a conic (here an ellipse)
  264.         // Of course, we note that if the point P belongs to the ellipsoid
  265.         // then ζ = 0 and the equation holds at point P since τ = 0 and υ = 0
  266.         final FieldVector3D<T> u     = planeNormal.orthogonal();
  267.         final FieldVector3D<T> v     = FieldVector3D.crossProduct(planeNormal, u).normalize();
  268.         final T xUOa    = u.getX().divide(a);
  269.         final T yUOb    = u.getY().divide(b);
  270.         final T zUOc    = u.getZ().divide(c);
  271.         final T xVOa    = v.getX().divide(a);
  272.         final T yVOb    = v.getY().divide(b);
  273.         final T zVOc    = v.getZ().divide(c);
  274.         final T xPOa    = planePoint.getX().divide(a);
  275.         final T yPOb    = planePoint.getY().divide(b);
  276.         final T zPOc    = planePoint.getZ().divide(c);
  277.         final T alpha   = xUOa.multiply(xUOa).add(yUOb.multiply(yUOb)).add(zUOc.multiply(zUOc));
  278.         final T beta    = xVOa.multiply(xVOa).add(yVOb.multiply(yVOb)).add(zVOc.multiply(zVOc));
  279.         final T gamma   = alpha.linearCombination(xUOa, xVOa, yUOb, yVOb, zUOc, zVOc);
  280.         final T delta   = alpha.linearCombination(xPOa, xUOa, yPOb, yUOb, zPOc, zUOc);
  281.         final T epsilon = alpha.linearCombination(xPOa, xVOa, yPOb, yVOb, zPOc, zVOc);
  282.         final T zeta    = alpha.linearCombination(xPOa, xPOa, yPOb, yPOb, zPOc, zPOc, one, one.negate());

  283.         // reduce the general equation α τ² + β υ² + 2 γ τυ + 2 δ τ + 2 ε υ + ζ = 0
  284.         // to canonical form (λ/l)² + (μ/m)² = 1
  285.         // using a coordinates change
  286.         //       τ = τC + λ cosθ - μ sinθ
  287.         //       υ = υC + λ sinθ + μ cosθ
  288.         // or equivalently
  289.         //       λ =   (τ - τC) cosθ + (υ - υC) sinθ
  290.         //       μ = - (τ - τC) sinθ + (υ - υC) cosθ
  291.         // τC and υC are the coordinates of the 2D ellipse center with respect to P
  292.         // 2l and 2m and are the axes lengths (major or minor depending on which one is greatest)
  293.         // θ is the angle of the 2D ellipse axis corresponding to axis with length 2l

  294.         // choose θ in order to cancel the coupling term in λμ
  295.         // expanding the general equation, we get: A λ² + B μ² + C λμ + D λ + E μ + F = 0
  296.         // with C = 2[(β - α) cosθ sinθ + γ (cos²θ - sin²θ)]
  297.         // hence the term is cancelled when θ = arctan(t), with γ t² + (α - β) t - γ = 0
  298.         // As the solutions of the quadratic equation obey t₁t₂ = -1, they correspond to
  299.         // angles θ in quadrature to each other. Selecting one solution or the other simply
  300.         // exchanges the principal axes. As we don't care about which axis we want as the
  301.         // first one, we select an arbitrary solution
  302.         final T tanTheta;
  303.         if (FastMath.abs(gamma.getReal()) < Precision.SAFE_MIN) {
  304.             tanTheta = zero;
  305.         } else {
  306.             final T bMA = beta.subtract(alpha);
  307.             tanTheta = (bMA.getReal() >= 0) ?
  308.                        gamma.multiply(-2).divide(bMA.add(FastMath.sqrt(bMA.multiply(bMA).add(gamma.multiply(gamma).multiply(4))))) :
  309.                        gamma.multiply(-2).divide(bMA.subtract(FastMath.sqrt(bMA.multiply(bMA).add(gamma.multiply(gamma).multiply(4)))));
  310.         }
  311.         final T tan2   = tanTheta.multiply(tanTheta);
  312.         final T cos2   = tan2.add(1).reciprocal();
  313.         final T sin2   = tan2.multiply(cos2);
  314.         final T cosSin = tanTheta.multiply(cos2);
  315.         final T cos    = FastMath.sqrt(cos2);
  316.         final T sin    = tanTheta.multiply(cos);

  317.         // choose τC and υC in order to cancel the linear terms in λ and μ
  318.         // expanding the general equation, we get: A λ² + B μ² + C λμ + D λ + E μ + F = 0
  319.         // with D = 2[ (α τC + γ υC + δ) cosθ + (γ τC + β υC + ε) sinθ]
  320.         //      E = 2[-(α τC + γ υC + δ) sinθ + (γ τC + β υC + ε) cosθ]
  321.         // θ can be eliminated by combining the equations
  322.         // D cosθ - E sinθ = 2[α τC + γ υC + δ]
  323.         // E cosθ + D sinθ = 2[γ τC + β υC + ε]
  324.         // hence the terms D and E are both cancelled (regardless of θ) when
  325.         //     τC = (β δ - γ ε) / (γ² - α β)
  326.         //     υC = (α ε - γ δ) / (γ² - α β)
  327.         final T invDenom = gamma.linearCombination(gamma, gamma,   alpha.negate(), beta).reciprocal();
  328.         final T tauC     = gamma.linearCombination(beta,  delta,   gamma.negate(), epsilon).multiply(invDenom);
  329.         final T nuC      = gamma.linearCombination(alpha, epsilon, gamma.negate(), delta).multiply(invDenom);

  330.         // compute l and m
  331.         // expanding the general equation, we get: A λ² + B μ² + C λμ + D λ + E μ + F = 0
  332.         // with A = α cos²θ + β sin²θ + 2 γ cosθ sinθ
  333.         //      B = α sin²θ + β cos²θ - 2 γ cosθ sinθ
  334.         //      F = α τC² + β υC² + 2 γ τC υC + 2 δ τC + 2 ε υC + ζ
  335.         // hence we compute directly l = √(-F/A) and m = √(-F/B)
  336.         final T twogcs = gamma.multiply(cosSin).multiply(2);
  337.         final T bigA   = alpha.multiply(cos2).add(beta.multiply(sin2)).add(twogcs);
  338.         final T bigB   = alpha.multiply(sin2).add(beta.multiply(cos2)).subtract(twogcs);
  339.         final T bigFN  = alpha.multiply(tauC).add(gamma.multiply(nuC).add(delta).multiply(2)).multiply(tauC).
  340.                          add(beta.multiply(nuC).add(epsilon.multiply(2)).multiply(nuC)).
  341.                          add(zeta).
  342.                          negate();
  343.         final T l      = FastMath.sqrt(bigFN.divide(bigA));
  344.         final T m      = FastMath.sqrt(bigFN.divide(bigB));
  345.         if (l.add(m).isNaN()) {
  346.             // the plane does not intersect the ellipsoid
  347.             return null;
  348.         }

  349.         if (l.subtract(m).getReal() > 0) {
  350.             return new FieldEllipse<>(new FieldVector3D<>(tauC.getField().getOne(), planePoint, tauC, u, nuC, v),
  351.                                       new FieldVector3D<>(cos,          u, sin, v),
  352.                                       new FieldVector3D<>(sin.negate(), u, cos, v),
  353.                                       l, m, frame);
  354.         } else {
  355.             return new FieldEllipse<>(new FieldVector3D<>(tauC.getField().getOne(), planePoint, tauC, u, nuC, v),
  356.                                       new FieldVector3D<>(sin, u, cos.negate(), v),
  357.                                       new FieldVector3D<>(cos, u, sin,          v),
  358.                                       m, l, frame);
  359.         }

  360.     }

  361.     /** Find a point on ellipsoid limb, as seen by an external observer.
  362.      * @param observer observer position in ellipsoid frame
  363.      * @param outside point outside ellipsoid in ellipsoid frame, defining the phase around limb
  364.      * @return point on ellipsoid limb
  365.      * @exception MathRuntimeException if ellipsoid center, observer and outside
  366.      * points are aligned
  367.      * @since 7.1
  368.      */
  369.     public Vector3D pointOnLimb(final Vector3D observer, final Vector3D outside)
  370.         throws MathRuntimeException {

  371.         // There is no limb if we are inside the ellipsoid
  372.         if (isInside(observer)) {
  373.             throw new OrekitException(OrekitMessages.POINT_INSIDE_ELLIPSOID);
  374.         }
  375.         // Cut the ellipsoid, to find an elliptical plane section
  376.         final Vector3D normal  = Vector3D.crossProduct(observer, outside);
  377.         final Ellipse  section = getPlaneSection(Vector3D.ZERO, normal);

  378.         // the point on limb is tangential to the ellipse
  379.         // if T(xt, yt) is an ellipse point at which the tangent is drawn
  380.         // if O(xo, yo) is a point outside of the ellipse belonging to the tangent at T,
  381.         // then the two following equations holds:
  382.         // (1) a² yt²   + b² xt²   = a² b²  (T belongs to the ellipse)
  383.         // (2) a² yt yo + b² xt xo = a² b²  (TP is tangent to the ellipse)
  384.         // using the second equation to eliminate yt from the first equation, we get
  385.         // b² (a² - xt xo)² + a² xt² yo² = a⁴ yo²
  386.         // (3) (a² yo² + b² xo²) xt² - 2 a² b² xo xt + a⁴ (b² - yo²) = 0
  387.         // which can easily be solved for xt

  388.         // To avoid numerical errors, the x and y coordinates in the ellipse plane are normalized using:
  389.         // x' = x / a and y' = y / b
  390.         //
  391.         // This gives:
  392.         // (1) y't² + x't² = 1
  393.         // (2) y't y'o + x't x'o = 1
  394.         //
  395.         // And finally:
  396.         // (3) (x'o² + y'o²) x't² - 2 x't x'o + 1 - y'o² = 0
  397.         //
  398.         // Solving for x't, we get the reduced discriminant:
  399.         // delta' = beta'² - alpha' * gamma'
  400.         //
  401.         // With:
  402.         // beta' = x'o
  403.         // alpha' = x'o² + y'o²
  404.         // gamma' = 1 - y'o²
  405.         //
  406.         // Simplifying to  cancel a term of x'o².
  407.         // delta' = y'o² (x'o² + y'o² - 1) = y'o² (alpha' - 1)
  408.         //
  409.         // After solving for xt1, xt2 using (3) the values are substituted into (2) to
  410.         // compute yt1, yt2. Then terms of x'o may be canceled from the expressions for
  411.         // yt1 and yt2. Additionally a point discontinuity is removed at y'o=0 from both
  412.         // yt1 and yt2.
  413.         //
  414.         // y't1 = (y'o - x'o d) / (x'o² + y'o²)
  415.         // y't2 = (x'o y'o + d) / (x + sqrt(delta'))
  416.         //
  417.         // where:
  418.         // d = sign(y'o) sqrt(alpha' - 1)

  419.         // Get the point in ellipse plane frame (2D)
  420.         final Vector2D observer2D = section.toPlane(observer);

  421.         // Normalize and compute intermediary terms
  422.         final double ap = section.getA();
  423.         final double bp = section.getB();
  424.         final double xpo = observer2D.getX() / ap;
  425.         final double ypo = observer2D.getY() / bp;
  426.         final double xpo2 = xpo * xpo;
  427.         final double ypo2 = ypo * ypo;
  428.         final double   alphap      = ypo2 + xpo2;
  429.         final double   gammap      = 1. - ypo2;

  430.         // Compute the roots
  431.         // We know there are two solutions as we already checked the point is outside ellipsoid
  432.         final double sqrt = FastMath.sqrt(alphap - 1);
  433.         final double sqrtp = FastMath.abs(ypo) * sqrt;
  434.         final double sqrtSigned = FastMath.copySign(sqrt, ypo);

  435.         // Compute the roots (ordered by value)
  436.         final double   xpt1;
  437.         final double   xpt2;
  438.         final double   ypt1;
  439.         final double   ypt2;
  440.         if (xpo > 0) {
  441.             final double s = xpo + sqrtp;
  442.             // xpt1 = (beta' + sqrt(delta')) / alpha' (with beta' = x'o)
  443.             xpt1 = s / alphap;
  444.             // x't2 = gamma' / (beta' + sqrt(delta')) since x't1 * x't2 = gamma' / alpha'
  445.             xpt2 = gammap / s;
  446.             // Get the corresponding values of y't
  447.             ypt1 = (ypo - xpo * sqrtSigned) / alphap;
  448.             ypt2 = (xpo * ypo + sqrtSigned) / s;
  449.         } else {
  450.             final double s = xpo - sqrtp;
  451.             // x't1 and x't2 are reverted compared to previous solution
  452.             xpt1 = gammap / s;
  453.             xpt2 = s / alphap;
  454.             // Get the corresponding values of y't
  455.             ypt2 = (ypo + xpo * sqrtSigned) / alphap;
  456.             ypt1 = (xpo * ypo - sqrtSigned) / s;
  457.         }

  458.         // De-normalize and express the two solutions in 3D
  459.         final Vector3D tp1 = section.toSpace(new Vector2D(ap * xpt1, bp * ypt1));
  460.         final Vector3D tp2 = section.toSpace(new Vector2D(ap * xpt2, bp * ypt2));

  461.         // Return the limb point in the direction of the outside point
  462.         return Vector3D.distance(tp1, outside) <= Vector3D.distance(tp2, outside) ? tp1 : tp2;

  463.     }

  464.     /** Find a point on ellipsoid limb, as seen by an external observer.
  465.      * @param observer observer position in ellipsoid frame
  466.      * @param outside point outside ellipsoid in ellipsoid frame, defining the phase around limb
  467.      * @return point on ellipsoid limb
  468.      * @exception MathRuntimeException if ellipsoid center, observer and outside
  469.      * points are aligned
  470.      * @param <T> the type of the field elements
  471.      * @since 12.0
  472.      */
  473.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> pointOnLimb(final FieldVector3D<T> observer, final FieldVector3D<T> outside)
  474.         throws MathRuntimeException {

  475.         // There is no limb if we are inside the ellipsoid
  476.         if (isInside(observer)) {
  477.             throw new OrekitException(OrekitMessages.POINT_INSIDE_ELLIPSOID);
  478.         }
  479.         // Cut the ellipsoid, to find an elliptical plane section
  480.         final FieldVector3D<T> normal  = FieldVector3D.crossProduct(observer, outside);
  481.         final FieldEllipse<T>  section = getPlaneSection(FieldVector3D.getZero(observer.getX().getField()), normal);

  482.         // the point on limb is tangential to the ellipse
  483.         // if T(xt, yt) is an ellipse point at which the tangent is drawn
  484.         // if O(xo, yo) is a point outside of the ellipse belonging to the tangent at T,
  485.         // then the two following equations holds:
  486.         // (1) a² yt²   + b² xt²   = a² b²  (T belongs to the ellipse)
  487.         // (2) a² yt yo + b² xt xo = a² b²  (TP is tangent to the ellipse)
  488.         // using the second equation to eliminate yt from the first equation, we get
  489.         // b² (a² - xt xo)² + a² xt² yo² = a⁴ yo²
  490.         // (3) (a² yo² + b² xo²) xt² - 2 a² b² xo xt + a⁴ (b² - yo²) = 0
  491.         // which can easily be solved for xt

  492.         // To avoid numerical errors, the x and y coordinates in the ellipse plane are normalized using:
  493.         // x' = x / a and y' = y / b
  494.         //
  495.         // This gives:
  496.         // (1) y't² + x't² = 1
  497.         // (2) y't y'o + x't x'o = 1
  498.         //
  499.         // And finally:
  500.         // (3) (x'o² + y'o²) x't² - 2 x't x'o + 1 - y'o² = 0
  501.         //
  502.         // Solving for x't, we get the reduced discriminant:
  503.         // delta' = beta'² - alpha' * gamma'
  504.         //
  505.         // With:
  506.         // beta' = x'o
  507.         // alpha' = x'o² + y'o²
  508.         // gamma' = 1 - y'o²
  509.         //
  510.         // Simplifying to  cancel a term of x'o².
  511.         // delta' = y'o² (x'o² + y'o² - 1) = y'o² (alpha' - 1)
  512.         //
  513.         // After solving for xt1, xt2 using (3) the values are substituted into (2) to
  514.         // compute yt1, yt2. Then terms of x'o may be canceled from the expressions for
  515.         // yt1 and yt2. Additionally a point discontinuity is removed at y'o=0 from both
  516.         // yt1 and yt2.
  517.         //
  518.         // y't1 = (y'o - x'o d) / (x'o² + y'o²)
  519.         // y't2 = (x'o y'o + d) / (x + sqrt(delta'))
  520.         //
  521.         // where:
  522.         // d = sign(y'o) sqrt(alpha' - 1)

  523.         // Get the point in ellipse plane frame (2D)
  524.         final FieldVector2D<T> observer2D = section.toPlane(observer);

  525.         // Normalize and compute intermediary terms
  526.         final T ap     = section.getA();
  527.         final T bp     = section.getB();
  528.         final T xpo    = observer2D.getX().divide(ap);
  529.         final T ypo    = observer2D.getY().divide(bp);
  530.         final T xpo2   = xpo.multiply(xpo);
  531.         final T ypo2   = ypo.multiply(ypo);
  532.         final T alphap = ypo2.add(xpo2);
  533.         final T gammap = ypo2.negate().add(1);

  534.         // Compute the roots
  535.         // We know there are two solutions as we already checked the point is outside ellipsoid
  536.         final T sqrt = FastMath.sqrt(alphap.subtract(1));
  537.         final T sqrtp = FastMath.abs(ypo).multiply(sqrt);
  538.         final T sqrtSigned = FastMath.copySign(sqrt, ypo);

  539.         // Compute the roots (ordered by value)
  540.         final T   xpt1;
  541.         final T   xpt2;
  542.         final T   ypt1;
  543.         final T   ypt2;
  544.         if (xpo.getReal() > 0) {
  545.             final T s = xpo.add(sqrtp);
  546.             // xpt1 = (beta' + sqrt(delta')) / alpha' (with beta' = x'o)
  547.             xpt1 = s.divide(alphap);
  548.             // x't2 = gamma' / (beta' + sqrt(delta')) since x't1 * x't2 = gamma' / alpha'
  549.             xpt2 = gammap.divide(s);
  550.             // Get the corresponding values of y't
  551.             ypt1 = ypo.subtract(xpo.multiply(sqrtSigned)).divide(alphap);
  552.             ypt2 = xpo.multiply(ypo).add(sqrtSigned).divide(s);
  553.         } else {
  554.             final T s = xpo.subtract(sqrtp);
  555.             // x't1 and x't2 are reverted compared to previous solution
  556.             xpt1 = gammap.divide(s);
  557.             xpt2 = s.divide(alphap);
  558.             // Get the corresponding values of y't
  559.             ypt2 = ypo.add(xpo.multiply(sqrtSigned)).divide(alphap);
  560.             ypt1 = xpo.multiply(ypo).subtract(sqrtSigned).divide(s);
  561.         }

  562.         // De-normalize and express the two solutions in 3D
  563.         final FieldVector3D<T> tp1 = section.toSpace(new FieldVector2D<>(ap.multiply(xpt1), bp.multiply(ypt1)));
  564.         final FieldVector3D<T> tp2 = section.toSpace(new FieldVector2D<>(ap.multiply(xpt2), bp.multiply(ypt2)));

  565.         // Return the limb point in the direction of the outside point
  566.         return FieldVector3D.distance(tp1, outside).subtract(FieldVector3D.distance(tp2, outside)).getReal() <= 0 ? tp1 : tp2;

  567.     }

  568. }