DSSTSolarRadiationPressure.java

  1. /* Copyright 2002-2020 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.propagation.semianalytical.dsst.forces;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.MathArrays;
  25. import org.hipparchus.util.MathUtils;
  26. import org.hipparchus.util.Precision;
  27. import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
  28. import org.orekit.forces.radiation.RadiationSensitive;
  29. import org.orekit.forces.radiation.SolarRadiationPressure;
  30. import org.orekit.propagation.FieldSpacecraftState;
  31. import org.orekit.propagation.SpacecraftState;
  32. import org.orekit.propagation.events.EventDetector;
  33. import org.orekit.propagation.events.FieldEventDetector;
  34. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  35. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  36. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  37. import org.orekit.utils.ParameterDriver;

  39. /** Solar radiation pressure contribution to the
  40.  *  {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
  41.  *  <p>
  42.  *  The solar radiation pressure acceleration is computed through the
  43.  *  acceleration model of
  44.  *  {@link org.orekit.forces.radiation.SolarRadiationPressure SolarRadiationPressure}.
  45.  *  </p>
  46.  *
  47.  *  @author Pascal Parraud
  48.  */
  49. public class DSSTSolarRadiationPressure extends AbstractGaussianContribution {

  51.     /** Reference distance for the solar radiation pressure (m). */
  52.     private static final double D_REF = 149597870000.0;

  54.     /** Reference solar radiation pressure at D_REF (N/m²). */
  55.     private static final double P_REF = 4.56e-6;

  57.     /** Threshold for the choice of the Gauss quadrature order. */
  58.     private static final double GAUSS_THRESHOLD = 1.0e-15;

  60.     /** Threshold for shadow equation. */
  61.     private static final double S_ZERO = 1.0e-6;

  63.     /** Prefix for the coefficient keys. */
  64.     private static final String PREFIX = "DSST-SRP-";

  66.     /** Sun model. */
  67.     private final ExtendedPVCoordinatesProvider sun;

  69.     /** Central Body radius. */
  70.     private final double                ae;

  72.     /** Spacecraft model for radiation acceleration computation. */
  73.     private final RadiationSensitive spacecraft;

  75.     /**
  76.      * Simple constructor with default reference values and spherical spacecraft.
  77.      * <p>
  78.      * When this constructor is used, the reference values are:
  79.      * </p>
  80.      * <ul>
  81.      * <li>d<sub>ref</sub> = 149597870000.0 m</li>
  82.      * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  83.      * </ul>
  84.      * <p>
  85.      * The spacecraft has a spherical shape.
  86.      * </p>
  87.      *
  88.      * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
  89.      * @param area cross sectional area of satellite
  90.      * @param sun Sun model
  91.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  92.      * @param mu central attraction coefficient
  93.      * @since 9.2
  94.      */
  95.     public DSSTSolarRadiationPressure(final double cr, final double area,
  96.                                       final ExtendedPVCoordinatesProvider sun,
  97.                                       final double equatorialRadius,
  98.                                       final double mu) {
  99.         this(D_REF, P_REF, cr, area, sun, equatorialRadius, mu);
  100.     }

  102.     /**
  103.      * Simple constructor with default reference values, but custom spacecraft.
  104.      * <p>
  105.      * When this constructor is used, the reference values are:
  106.      * </p>
  107.      * <ul>
  108.      * <li>d<sub>ref</sub> = 149597870000.0 m</li>
  109.      * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  110.      * </ul>
  111.      *
  112.      * @param sun Sun model
  113.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  114.      * @param spacecraft spacecraft model
  115.      * @param mu central attraction coefficient
  116.      * @since 9.2
  117.      */
  118.     public DSSTSolarRadiationPressure(final ExtendedPVCoordinatesProvider sun,
  119.                                       final double equatorialRadius,
  120.                                       final RadiationSensitive spacecraft,
  121.                                       final double mu) {
  122.         this(D_REF, P_REF, sun, equatorialRadius, spacecraft, mu);
  123.     }

  125.     /**
  126.      * Constructor with customizable reference values but spherical spacecraft.
  127.      * <p>
  128.      * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
  129.      * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
  130.      * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  131.      * N/m² solar radiation pressure.
  132.      * </p>
  133.      *
  134.      * @param dRef reference distance for the solar radiation pressure (m)
  135.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  136.      * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
  137.      * @param area cross sectional area of satellite
  138.      * @param sun Sun model
  139.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  140.      * @param mu central attraction coefficient
  141.      * @since 9.2
  142.      */
  143.     public DSSTSolarRadiationPressure(final double dRef, final double pRef,
  144.                                       final double cr, final double area,
  145.                                       final ExtendedPVCoordinatesProvider sun,
  146.                                       final double equatorialRadius,
  147.                                       final double mu) {

  149.         // cR being the DSST SRP coef and assuming a spherical spacecraft,
  150.         // the conversion is:
  151.         // cR = 1 + (1 - kA) * (1 - kR) * 4 / 9
  152.         // with kA arbitrary sets to 0
  153.         this(dRef, pRef, sun, equatorialRadius, new IsotropicRadiationSingleCoefficient(area, cr), mu);
  154.     }

  156.     /**
  157.      * Complete constructor.
  158.      * <p>
  159.      * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
  160.      * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
  161.      * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  162.      * N/m² solar radiation pressure.
  163.      * </p>
  164.      *
  165.      * @param dRef reference distance for the solar radiation pressure (m)
  166.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  167.      * @param sun Sun model
  168.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  169.      * @param spacecraft spacecraft model
  170.      * @param mu central attraction coefficient
  171.      * @since 9.2
  172.      */
  173.     public DSSTSolarRadiationPressure(final double dRef, final double pRef,
  174.                                       final ExtendedPVCoordinatesProvider sun,
  175.                                       final double equatorialRadius,
  176.                                       final RadiationSensitive spacecraft,
  177.                                       final double mu) {

  179.         //Call to the constructor from superclass using the numerical SRP model as ForceModel
  180.         super(PREFIX, GAUSS_THRESHOLD,
  181.               new SolarRadiationPressure(dRef, pRef, sun, equatorialRadius, spacecraft), mu);

  183.         this.sun  = sun;
  184.         this.ae   = equatorialRadius;
  185.         this.spacecraft = spacecraft;
  186.     }

  188.     /** Get spacecraft shape.
  189.      * @return the spacecraft shape.
  190.      */
  191.     public RadiationSensitive getSpacecraft() {
  192.         return spacecraft;
  193.     }

  195.     /** {@inheritDoc} */
  196.     public EventDetector[] getEventsDetectors() {
  197.         return null;
  198.     }

  200.     /** {@inheritDoc} */
  201.     @Override
  202.     public <T extends RealFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
  203.         return null;
  204.     }

  206.     /** {@inheritDoc} */
  207.     protected ParameterDriver[] getParametersDriversWithoutMu() {
  208.         return spacecraft.getRadiationParametersDrivers();
  209.     }

  211.     /** {@inheritDoc} */
  212.     protected double[] getLLimits(final SpacecraftState state, final AuxiliaryElements auxiliaryElements) {

  214.         // Default bounds without shadow [-PI, PI]
  215.         final double[] ll = {-FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0),
  216.                              FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0)};

  218.         // Direction cosines of the Sun in the equinoctial frame
  219.         final Vector3D sunDir = sun.getPVCoordinates(state.getDate(), state.getFrame()).getPosition().normalize();
  220.         final double alpha = sunDir.dotProduct(auxiliaryElements.getVectorF());
  221.         final double beta  = sunDir.dotProduct(auxiliaryElements.getVectorG());
  222.         final double gamma = sunDir.dotProduct(auxiliaryElements.getVectorW());

  224.         // Compute limits only if the perigee is close enough from the central body to be in the shadow
  225.         if (FastMath.abs(gamma * auxiliaryElements.getSma() * (1. - auxiliaryElements.getEcc())) < ae) {

  227.             // Compute the coefficients of the quartic equation in cos(L) 3.5-(2)
  228.             final double bet2 = beta * beta;
  229.             final double h2 = auxiliaryElements.getH() * auxiliaryElements.getH();
  230.             final double k2 = auxiliaryElements.getK() * auxiliaryElements.getK();
  231.             final double m  = ae / (auxiliaryElements.getSma() * auxiliaryElements.getB());
  232.             final double m2 = m * m;
  233.             final double m4 = m2 * m2;
  234.             final double bb = alpha * beta + m2 * auxiliaryElements.getH() * auxiliaryElements.getK();
  235.             final double b2 = bb * bb;
  236.             final double cc = alpha * alpha - bet2 + m2 * (k2 - h2);
  237.             final double dd = 1. - bet2 - m2 * (1. + h2);
  238.             final double[] a = new double[5];
  239.             a[0] = 4. * b2 + cc * cc;
  240.             a[1] = 8. * bb * m2 * auxiliaryElements.getH() + 4. * cc * m2 * auxiliaryElements.getK();
  241.             a[2] = -4. * b2 + 4. * m4 * h2 - 2. * cc * dd + 4. * m4 * k2;
  242.             a[3] = -8. * bb * m2 * auxiliaryElements.getH() - 4. * dd * m2 * auxiliaryElements.getK();
  243.             a[4] = -4. * m4 * h2 + dd * dd;
  244.             // Compute the real roots of the quartic equation 3.5-2
  245.             final double[] roots = new double[4];
  246.             final int nbRoots = realQuarticRoots(a, roots);
  247.             if (nbRoots > 0) {
  248.                 // Check for consistency
  249.                 boolean entryFound = false;
  250.                 boolean exitFound  = false;
  251.                 // Eliminate spurious roots
  252.                 for (int i = 0; i < nbRoots; i++) {
  253.                     final double cosL = roots[i];
  254.                     final double sL = FastMath.sqrt((1. - cosL) * (1. + cosL));
  255.                     // Check both angles: L and -L
  256.                     for (int j = -1; j <= 1; j += 2) {
  257.                         final double sinL = j * sL;
  258.                         final double cPhi = alpha * cosL + beta * sinL;
  259.                         // Is the angle on the shadow side of the central body (eq. 3.5-3) ?
  260.                         if (cPhi < 0.) {
  261.                             final double range = 1. + auxiliaryElements.getK() * cosL + auxiliaryElements.getH() * sinL;
  262.                             final double S  = 1. - m2 * range * range - cPhi * cPhi;
  263.                             // Is the shadow equation 3.5-1 satisfied ?
  264.                             if (FastMath.abs(S) < S_ZERO) {
  265.                                 // Is this the entry or exit angle ?
  266.                                 final double dSdL = m2 * range * (auxiliaryElements.getK() * sinL - auxiliaryElements.getH() * cosL) + cPhi * (alpha * sinL - beta * cosL);
  267.                                 if (dSdL > 0.) {
  268.                                     // Exit from shadow: 3.5-4
  269.                                     exitFound = true;
  270.                                     ll[0] = FastMath.atan2(sinL, cosL);
  271.                                 } else {
  272.                                     // Entry into shadow: 3.5-5
  273.                                     entryFound = true;
  274.                                     ll[1] = FastMath.atan2(sinL, cosL);
  275.                                 }
  276.                             }
  277.                         }
  278.                     }
  279.                 }
  280.                 // Must be one entry and one exit or none
  281.                 if (!(entryFound == exitFound)) {
  282.                     // entry or exit found but not both ! In this case, consider there is no eclipse...
  283.                     ll[0] = -FastMath.PI;
  284.                     ll[1] = FastMath.PI;
  285.                 }
  286.                 // Quadrature between L at exit and L at entry so Lexit must be lower than Lentry
  287.                 if (ll[0] > ll[1]) {
  288.                     // Keep the angles between [-2PI, 2PI]
  289.                     if (ll[1] < 0.) {
  290.                         ll[1] += 2. * FastMath.PI;
  291.                     } else {
  292.                         ll[0] -= 2. * FastMath.PI;
  293.                     }
  294.                 }
  295.             }
  296.         }
  297.         return ll;
  298.     }

  300.     /** {@inheritDoc} */
  301.     protected <T extends RealFieldElement<T>> T[] getLLimits(final FieldSpacecraftState<T> state,
  302.                                                              final FieldAuxiliaryElements<T> auxiliaryElements) {

  304.         final Field<T> field = state.getDate().getField();
  305.         final T zero = field.getZero();
  306.         final T one  = field.getOne();

  308.         // Default bounds without shadow [-PI, PI]
  309.         final T[] ll = MathArrays.buildArray(field, 2);
  310.         ll[0] = MathUtils.normalizeAngle(state.getLv(), zero).subtract(FastMath.PI);
  311.         ll[1] = MathUtils.normalizeAngle(state.getLv(), zero).add(FastMath.PI);

  313.         // Direction cosines of the Sun in the equinoctial frame
  314.         final FieldVector3D<T> sunDir = sun.getPVCoordinates(state.getDate(), state.getFrame()).getPosition().normalize();
  315.         final T alpha = sunDir.dotProduct(auxiliaryElements.getVectorF());
  316.         final T beta  = sunDir.dotProduct(auxiliaryElements.getVectorG());
  317.         final T gamma = sunDir.dotProduct(auxiliaryElements.getVectorW());

  319.         // Compute limits only if the perigee is close enough from the central body to be in the shadow
  320.         if (FastMath.abs(gamma.multiply(auxiliaryElements.getSma()).multiply(auxiliaryElements.getEcc().negate().add(one))).getReal() < ae) {

  322.             // Compute the coefficients of the quartic equation in cos(L) 3.5-(2)
  323.             final T bet2 = beta.multiply(beta);
  324.             final T h2 = auxiliaryElements.getH().multiply(auxiliaryElements.getH());
  325.             final T k2 = auxiliaryElements.getK().multiply(auxiliaryElements.getK());
  326.             final T m  = (auxiliaryElements.getSma().multiply(auxiliaryElements.getB())).divide(ae).reciprocal();
  327.             final T m2 = m.multiply(m);
  328.             final T m4 = m2.multiply(m2);
  329.             final T bb = alpha.multiply(beta).add(m2.multiply(auxiliaryElements.getH()).multiply(auxiliaryElements.getK()));
  330.             final T b2 = bb.multiply(bb);
  331.             final T cc = alpha.multiply(alpha).subtract(bet2).add(m2.multiply(k2.subtract(h2)));
  332.             final T dd = bet2.add(m2.multiply(h2.add(1.))).negate().add(one);
  333.             final T[] a = MathArrays.buildArray(field, 5);
  334.             a[0] = b2.multiply(4.).add(cc.multiply(cc));
  335.             a[1] = bb.multiply(8.).multiply(m2).multiply(auxiliaryElements.getH()).add(cc.multiply(4.).multiply(m2).multiply(auxiliaryElements.getK()));
  336.             a[2] = m4.multiply(h2).multiply(4.).subtract(cc.multiply(dd).multiply(2.)).add(m4.multiply(k2).multiply(4.)).subtract(b2.multiply(4.));
  337.             a[3] = auxiliaryElements.getH().multiply(m2).multiply(bb).multiply(8.).add(auxiliaryElements.getK().multiply(m2).multiply(dd).multiply(4.)).negate();
  338.             a[4] = dd.multiply(dd).subtract(m4.multiply(h2).multiply(4.));
  339.             // Compute the real roots of the quartic equation 3.5-2
  340.             final T[] roots = MathArrays.buildArray(field, 4);
  341.             final int nbRoots = realQuarticRoots(a, roots, field);
  342.             if (nbRoots > 0) {
  343.                 // Check for consistency
  344.                 boolean entryFound = false;
  345.                 boolean exitFound  = false;
  346.                 // Eliminate spurious roots
  347.                 for (int i = 0; i < nbRoots; i++) {
  348.                     final T cosL = roots[i];
  349.                     final T sL = FastMath.sqrt((cosL.negate().add(one)).multiply(cosL.add(one)));
  350.                     // Check both angles: L and -L
  351.                     for (int j = -1; j <= 1; j += 2) {
  352.                         final T sinL = sL.multiply(j);
  353.                         final T cPhi = cosL.multiply(alpha).add(sinL.multiply(beta));
  354.                         // Is the angle on the shadow side of the central body (eq. 3.5-3) ?
  355.                         if (cPhi.getReal() < 0.) {
  356.                             final T range = cosL.multiply(auxiliaryElements.getK()).add(sinL.multiply(auxiliaryElements.getH())).add(one);
  357.                             final T S  = (range.multiply(range).multiply(m2).add(cPhi.multiply(cPhi))).negate().add(1.);
  358.                             // Is the shadow equation 3.5-1 satisfied ?
  359.                             if (FastMath.abs(S).getReal() < S_ZERO) {
  360.                                 // Is this the entry or exit angle ?
  361.                                 final T dSdL = m2.multiply(range).multiply(auxiliaryElements.getK().multiply(sinL).subtract(auxiliaryElements.getH().multiply(cosL))).add(cPhi.multiply(alpha.multiply(sinL).subtract(beta.multiply(cosL))));
  362.                                 if (dSdL.getReal() > 0.) {
  363.                                     // Exit from shadow: 3.5-4
  364.                                     exitFound = true;
  365.                                     ll[0] = FastMath.atan2(sinL, cosL);
  366.                                 } else {
  367.                                     // Entry into shadow: 3.5-5
  368.                                     entryFound = true;
  369.                                     ll[1] = FastMath.atan2(sinL, cosL);
  370.                                 }
  371.                             }
  372.                         }
  373.                     }
  374.                 }
  375.                 // Must be one entry and one exit or none
  376.                 if (!(entryFound == exitFound)) {
  377.                     // entry or exit found but not both ! In this case, consider there is no eclipse...
  378.                     ll[0] = zero.add(-FastMath.PI);
  379.                     ll[1] = zero.add(FastMath.PI);
  380.                 }
  381.                 // Quadrature between L at exit and L at entry so Lexit must be lower than Lentry
  382.                 if (ll[0].getReal() > ll[1].getReal()) {
  383.                     // Keep the angles between [-2PI, 2PI]
  384.                     if (ll[1].getReal() < 0.) {
  385.                         ll[1] = ll[1].add(zero.add(2. * FastMath.PI));
  386.                     } else {
  387.                         ll[0] = ll[0].subtract(zero.add(2. * FastMath.PI));
  388.                     }
  389.                 }
  390.             }
  391.         }
  392.         return ll;
  393.     }

  395.     /** Get the central body equatorial radius.
  396.      *  @return central body equatorial radius (m)
  397.      */
  398.     public double getEquatorialRadius() {
  399.         return ae;
  400.     }

  402.     /**
  403.      * Compute the real roots of a quartic equation.
  404.      *
  405.      * <pre>
  406.      * a[0] * x⁴ + a[1] * x³ + a[2] * x² + a[3] * x + a[4] = 0.
  407.      * </pre>
  408.      *
  409.      * @param a the 5 coefficients
  410.      * @param y the real roots
  411.      * @return the number of real roots
  412.      */
  413.     private int realQuarticRoots(final double[] a, final double[] y) {
  414.         /* Treat the degenerate quartic as cubic */
  415.         if (Precision.equals(a[0], 0.)) {
  416.             final double[] aa = new double[a.length - 1];
  417.             System.arraycopy(a, 1, aa, 0, aa.length);
  418.             return realCubicRoots(aa, y);
  419.         }

  421.         // Transform coefficients
  422.         final double b  = a[1] / a[0];
  423.         final double c  = a[2] / a[0];
  424.         final double d  = a[3] / a[0];
  425.         final double e  = a[4] / a[0];
  426.         final double bh = b * 0.5;

  428.         // Solve resolvant cubic
  429.         final double[] z3 = new double[3];
  430.         final int i3 = realCubicRoots(new double[] {1.0, -c, b * d - 4. * e, e * (4. * c - b * b) - d * d}, z3);
  431.         if (i3 == 0) {
  432.             return 0;
  433.         }

  435.         // Largest real root of resolvant cubic
  436.         final double z = z3[0];

  438.         // Compute auxiliary quantities
  439.         final double zh = z * 0.5;
  440.         final double p  = FastMath.max(z + bh * bh - c, 0.);
  441.         final double q  = FastMath.max(zh * zh - e, 0.);
  442.         final double r  = bh * z - d;
  443.         final double pp = FastMath.sqrt(p);
  444.         final double qq = FastMath.copySign(FastMath.sqrt(q), r);

  446.         // Solve quadratic factors of quartic equation
  447.         final double[] y1 = new double[2];
  448.         final int n1 = realQuadraticRoots(new double[] {1.0, bh - pp, zh - qq}, y1);
  449.         final double[] y2 = new double[2];
  450.         final int n2 = realQuadraticRoots(new double[] {1.0, bh + pp, zh + qq}, y2);

  452.         if (n1 == 2) {
  453.             if (n2 == 2) {
  454.                 y[0] = y1[0];
  455.                 y[1] = y1[1];
  456.                 y[2] = y2[0];
  457.                 y[3] = y2[1];
  458.                 return 4;
  459.             } else {
  460.                 y[0] = y1[0];
  461.                 y[1] = y1[1];
  462.                 return 2;
  463.             }
  464.         } else {
  465.             if (n2 == 2) {
  466.                 y[0] = y2[0];
  467.                 y[1] = y2[1];
  468.                 return 2;
  469.             } else {
  470.                 return 0;
  471.             }
  472.         }
  473.     }

  475.     /**
  476.      * Compute the real roots of a quartic equation.
  477.      *
  478.      * <pre>
  479.      * a[0] * x⁴ + a[1] * x³ + a[2] * x² + a[3] * x + a[4] = 0.
  480.      * </pre>
  481.      * @param <T> the type of the field elements
  482.      *
  483.      * @param a the 5 coefficients
  484.      * @param y the real roots
  485.      * @param field field of elements
  486.      * @return the number of real roots
  487.      */
  488.     private <T extends RealFieldElement<T>> int realQuarticRoots(final T[] a, final T[] y,
  489.                                                                  final Field<T> field) {

  491.         final T zero = field.getZero();

  493.         // Treat the degenerate quartic as cubic
  494.         if (Precision.equals(a[0].getReal(), 0.)) {
  495.             final T[] aa = MathArrays.buildArray(field, a.length - 1);
  496.             System.arraycopy(a, 1, aa, 0, aa.length);
  497.             return realCubicRoots(aa, y, field);
  498.         }

  500.         // Transform coefficients
  501.         final T b  = a[1].divide(a[0]);
  502.         final T c  = a[2].divide(a[0]);
  503.         final T d  = a[3].divide(a[0]);
  504.         final T e  = a[4].divide(a[0]);
  505.         final T bh = b.multiply(0.5);

  507.         // Solve resolvant cubic
  508.         final T[] z3 = MathArrays.buildArray(field, 3);
  509.         final T[] i = MathArrays.buildArray(field, 4);
  510.         i[0] = zero.add(1.0);
  511.         i[1] = c.negate();
  512.         i[2] = b.multiply(d).subtract(e.multiply(4.0));
  513.         i[3] = e.multiply(c.multiply(4.).subtract(b.multiply(b))).subtract(d.multiply(d));
  514.         final int i3 = realCubicRoots(i, z3, field);
  515.         if (i3 == 0) {
  516.             return 0;
  517.         }

  519.         // Largest real root of resolvant cubic
  520.         final T z = z3[0];

  522.         // Compute auxiliary quantities
  523.         final T zh = z.multiply(0.5);
  524.         final T p  = FastMath.max(z.add(bh.multiply(bh)).subtract(c), zero);
  525.         final T q  = FastMath.max(zh.multiply(zh).subtract(e), zero);
  526.         final T r  = bh.multiply(z).subtract(d);
  527.         final T pp = FastMath.sqrt(p);
  528.         final T qq = FastMath.copySign(FastMath.sqrt(q), r);

  530.         // Solve quadratic factors of quartic equation
  531.         final T[] y1 = MathArrays.buildArray(field, 2);
  532.         final T[] n = MathArrays.buildArray(field, 3);
  533.         n[0] = zero.add(1.0);
  534.         n[1] = bh.subtract(pp);
  535.         n[2] = zh.subtract(qq);
  536.         final int n1 = realQuadraticRoots(n, y1);
  537.         final T[] y2 = MathArrays.buildArray(field, 2);
  538.         final T[] nn = MathArrays.buildArray(field, 3);
  539.         nn[0] = zero.add(1.0);
  540.         nn[1] = bh.add(pp);
  541.         nn[2] = zh.add(qq);
  542.         final int n2 = realQuadraticRoots(nn, y2);

  544.         if (n1 == 2) {
  545.             if (n2 == 2) {
  546.                 y[0] = y1[0];
  547.                 y[1] = y1[1];
  548.                 y[2] = y2[0];
  549.                 y[3] = y2[1];
  550.                 return 4;
  551.             } else {
  552.                 y[0] = y1[0];
  553.                 y[1] = y1[1];
  554.                 return 2;
  555.             }
  556.         } else {
  557.             if (n2 == 2) {
  558.                 y[0] = y2[0];
  559.                 y[1] = y2[1];
  560.                 return 2;
  561.             } else {
  562.                 return 0;
  563.             }
  564.         }
  565.     }

  567.     /**
  568.      * Compute the real roots of a cubic equation.
  569.      *
  570.      * <pre>
  571.      * a[0] * x³ + a[1] * x² + a[2] * x + a[3] = 0.
  572.      * </pre>
  573.      *
  574.      * @param a the 4 coefficients
  575.      * @param y the real roots sorted in descending order
  576.      * @return the number of real roots
  577.      */
  578.     private int realCubicRoots(final double[] a, final double[] y) {
  579.         if (Precision.equals(a[0], 0.)) {
  580.             // Treat the degenerate cubic as quadratic
  581.             final double[] aa = new double[a.length - 1];
  582.             System.arraycopy(a, 1, aa, 0, aa.length);
  583.             return realQuadraticRoots(aa, y);
  584.         }

  586.         // Transform coefficients
  587.         final double b  = -a[1] / (3. * a[0]);
  588.         final double c  =  a[2] / a[0];
  589.         final double d  =  a[3] / a[0];
  590.         final double b2 =  b * b;
  591.         final double p  =  b2 - c / 3.;
  592.         final double q  =  b * (b2 - c * 0.5) - d * 0.5;

  594.         // Compute discriminant
  595.         final double disc = p * p * p - q * q;

  597.         if (disc < 0.) {
  598.             // One real root
  599.             final double alpha  = q + FastMath.copySign(FastMath.sqrt(-disc), q);
  600.             final double cbrtAl = FastMath.cbrt(alpha);
  601.             final double cbrtBe = p / cbrtAl;

  603.             if (p < 0.) {
  604.                 y[0] = b + 2. * q / (cbrtAl * cbrtAl + cbrtBe * cbrtBe - p);
  605.             } else if (p > 0.) {
  606.                 y[0] = b + cbrtAl + cbrtBe;
  607.             } else {
  608.                 y[0] = b + cbrtAl;
  609.             }

  611.             return 1;

  613.         } else if (disc > 0.) {
  614.             // Three distinct real roots
  615.             final double phi = FastMath.atan2(FastMath.sqrt(disc), q) / 3.;
  616.             final double sqP = 2.0 * FastMath.sqrt(p);

  618.             y[0] = b + sqP * FastMath.cos(phi);
  619.             y[1] = b - sqP * FastMath.cos(FastMath.PI / 3. + phi);
  620.             y[2] = b - sqP * FastMath.cos(FastMath.PI / 3. - phi);

  622.             return 3;

  624.         } else {
  625.             // One distinct and two equals real roots
  626.             final double cbrtQ = FastMath.cbrt(q);
  627.             final double root1 = b + 2. * cbrtQ;
  628.             final double root2 = b - cbrtQ;
  629.             if (q < 0.) {
  630.                 y[0] = root2;
  631.                 y[1] = root2;
  632.                 y[2] = root1;
  633.             } else {
  634.                 y[0] = root1;
  635.                 y[1] = root2;
  636.                 y[2] = root2;
  637.             }

  639.             return 3;
  640.         }
  641.     }

  643.     /**
  644.      * Compute the real roots of a cubic equation.
  645.      *
  646.      * <pre>
  647.      * a[0] * x³ + a[1] * x² + a[2] * x + a[3] = 0.
  648.      * </pre>
  649.      *
  650.      * @param <T> the type of the field elements
  651.      * @param a the 4 coefficients
  652.      * @param y the real roots sorted in descending order
  653.      * @param field field of elements
  654.      * @return the number of real roots
  655.      */
  656.     private <T extends RealFieldElement<T>> int realCubicRoots(final T[] a, final T[] y,
  657.                                                                final Field<T> field) {

  659.         if (Precision.equals(a[0].getReal(), 0.)) {
  660.             // Treat the degenerate cubic as quadratic
  661.             final T[] aa = MathArrays.buildArray(field, a.length - 1);
  662.             System.arraycopy(a, 1, aa, 0, aa.length);
  663.             return realQuadraticRoots(aa, y);
  664.         }

  666.         // Transform coefficients
  667.         final T b  =  a[1].divide(a[0].multiply(3.)).negate();
  668.         final T c  =  a[2].divide(a[0]);
  669.         final T d  =  a[3].divide(a[0]);
  670.         final T b2 =  b.multiply(b);
  671.         final T p  =  b2.subtract(c.divide(3.));
  672.         final T q  =  b.multiply(b2.subtract(c.multiply(0.5))).subtract(d.multiply(0.5));

  674.         // Compute discriminant
  675.         final T disc = p.multiply(p).multiply(p).subtract(q.multiply(q));

  677.         if (disc.getReal() < 0.) {
  678.             // One real root
  679.             final T alpha  = FastMath.copySign(FastMath.sqrt(disc.negate()), q).add(q);
  680.             final T cbrtAl = FastMath.cbrt(alpha);
  681.             final T cbrtBe = p.divide(cbrtAl);

  683.             if (p .getReal() < 0.) {
  684.                 y[0] = q.divide(cbrtAl.multiply(cbrtAl).add(cbrtBe.multiply(cbrtBe)).subtract(p)).multiply(2.).add(b);
  685.             } else if (p.getReal() > 0.) {
  686.                 y[0] = b.add(cbrtAl).add(cbrtBe);
  687.             } else {
  688.                 y[0] = b.add(cbrtAl);
  689.             }

  691.             return 1;

  693.         } else if (disc.getReal() > 0.) {
  694.             // Three distinct real roots
  695.             final T phi = FastMath.atan2(FastMath.sqrt(disc), q).divide(3.);
  696.             final T sqP = FastMath.sqrt(p).multiply(2.);

  698.             y[0] = b.add(sqP.multiply(FastMath.cos(phi)));
  699.             y[1] = b.subtract(sqP.multiply(FastMath.cos(phi.add(FastMath.PI / 3.))));
  700.             y[2] = b.subtract(sqP.multiply(FastMath.cos(phi.negate().add(FastMath.PI / 3.))));

  702.             return 3;

  704.         } else {
  705.             // One distinct and two equals real roots
  706.             final T cbrtQ = FastMath.cbrt(q);
  707.             final T root1 = b.add(cbrtQ.multiply(2.0));
  708.             final T root2 = b.subtract(cbrtQ);
  709.             if (q.getReal() < 0.) {
  710.                 y[0] = root2;
  711.                 y[1] = root2;
  712.                 y[2] = root1;
  713.             } else {
  714.                 y[0] = root1;
  715.                 y[1] = root2;
  716.                 y[2] = root2;
  717.             }

  719.             return 3;
  720.         }
  721.     }

  723.     /**
  724.      * Compute the real roots of a quadratic equation.
  725.      *
  726.      * <pre>
  727.      * a[0] * x² + a[1] * x + a[2] = 0.
  728.      * </pre>
  729.      *
  730.      * @param a the 3 coefficients
  731.      * @param y the real roots sorted in descending order
  732.      * @return the number of real roots
  733.      */
  734.     private int realQuadraticRoots(final double[] a, final double[] y) {
  735.         if (Precision.equals(a[0], 0.)) {
  736.             // Degenerate quadratic
  737.             if (Precision.equals(a[1], 0.)) {
  738.                 // Degenerate linear equation: no real roots
  739.                 return 0;
  740.             }
  741.             // Linear equation: one real root
  742.             y[0] = -a[2] / a[1];
  743.             return 1;
  744.         }

  746.         // Transform coefficients
  747.         final double b = -0.5 * a[1] / a[0];
  748.         final double c =  a[2] / a[0];

  750.         // Compute discriminant
  751.         final double d =  b * b - c;

  753.         if (d < 0.) {
  754.             // No real roots
  755.             return 0;
  756.         } else if (d > 0.) {
  757.             // Two distinct real roots
  758.             final double y0 = b + FastMath.copySign(FastMath.sqrt(d), b);
  759.             final double y1 = c / y0;
  760.             y[0] = FastMath.max(y0, y1);
  761.             y[1] = FastMath.min(y0, y1);
  762.             return 2;
  763.         } else {
  764.             // Discriminant is zero: two equal real roots
  765.             y[0] = b;
  766.             y[1] = b;
  767.             return 2;
  768.         }
  769.     }

  771.     /**
  772.      * Compute the real roots of a quadratic equation.
  773.      *
  774.      * <pre>
  775.      * a[0] * x² + a[1] * x + a[2] = 0.
  776.      * </pre>
  777.      *
  778.      * @param <T> the type of the field elements
  779.      * @param a the 3 coefficients
  780.      * @param y the real roots sorted in descending order
  781.      * @return the number of real roots
  782.      */
  783.     private <T extends RealFieldElement<T>> int realQuadraticRoots(final T[] a, final T[] y) {

  785.         if (Precision.equals(a[0].getReal(), 0.)) {
  786.             // Degenerate quadratic
  787.             if (Precision.equals(a[1].getReal(), 0.)) {
  788.                 // Degenerate linear equation: no real roots
  789.                 return 0;
  790.             }
  791.             // Linear equation: one real root
  792.             y[0] = a[2].divide(a[1]).negate();
  793.             return 1;
  794.         }

  796.         // Transform coefficients
  797.         final T b = a[1].divide(a[0]).multiply(0.5).negate();
  798.         final T c =  a[2].divide(a[0]);

  800.         // Compute discriminant
  801.         final T d =  b.multiply(b).subtract(c);

  803.         if (d.getReal() < 0.) {
  804.             // No real roots
  805.             return 0;
  806.         } else if (d.getReal() > 0.) {
  807.             // Two distinct real roots
  808.             final T y0 = b.add(FastMath.copySign(FastMath.sqrt(d), b));
  809.             final T y1 = c.divide(y0);
  810.             y[0] = FastMath.max(y0, y1);
  811.             y[1] = FastMath.min(y0, y1);
  812.             return 2;
  813.         } else {
  814.             // Discriminant is zero: two equal real roots
  815.             y[0] = b;
  816.             y[1] = b;
  817.             return 2;
  818.         }
  819.     }

  821. }
Created with JaCoCo 0.8.6.202009150832