DSSTSolarRadiationPressure.java

  1. /* Copyright 2002-2023 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 java.util.List;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  184.         this.sun  = sun;
  185.         this.ae   = centralBody.getEquatorialRadius();
  186.         this.spacecraft = spacecraft;
  187.     }

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

  196.     /** {@inheritDoc} */
  197.     protected List<ParameterDriver> getParametersDriversWithoutMu() {
  198.         return spacecraft.getRadiationParametersDrivers();
  199.     }

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

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

  208.         // Direction cosines of the Sun in the equinoctial frame
  209.         final Vector3D sunDir = sun.getPosition(state.getDate(), state.getFrame()).normalize();
  210.         final double alpha = sunDir.dotProduct(auxiliaryElements.getVectorF());
  211.         final double beta  = sunDir.dotProduct(auxiliaryElements.getVectorG());
  212.         final double gamma = sunDir.dotProduct(auxiliaryElements.getVectorW());

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

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

  290.     /** {@inheritDoc} */
  291.     protected <T extends CalculusFieldElement<T>> T[] getLLimits(final FieldSpacecraftState<T> state,
  292.                                                              final FieldAuxiliaryElements<T> auxiliaryElements) {

  294.         final Field<T> field = state.getDate().getField();
  295.         final T zero = field.getZero();
  296.         final T one  = field.getOne();
  297.         final T pi   = one.getPi();

  299.         // Default bounds without shadow [-PI, PI]
  300.         final T[] ll = MathArrays.buildArray(field, 2);
  301.         ll[0] = MathUtils.normalizeAngle(state.getLv(), zero).subtract(pi);
  302.         ll[1] = MathUtils.normalizeAngle(state.getLv(), zero).add(pi);

  304.         // Direction cosines of the Sun in the equinoctial frame
  305.         final FieldVector3D<T> sunDir = sun.getPosition(state.getDate(), state.getFrame()).normalize();
  306.         final T alpha = sunDir.dotProduct(auxiliaryElements.getVectorF());
  307.         final T beta  = sunDir.dotProduct(auxiliaryElements.getVectorG());
  308.         final T gamma = sunDir.dotProduct(auxiliaryElements.getVectorW());

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

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

  386.     /** Get the central body equatorial radius.
  387.      *  @return central body equatorial radius (m)
  388.      */
  389.     public double getEquatorialRadius() {
  390.         return ae;
  391.     }

  393.     /**
  394.      * Compute the real roots of a quartic equation.
  395.      *
  396.      * <pre>
  397.      * a[0] * x⁴ + a[1] * x³ + a[2] * x² + a[3] * x + a[4] = 0.
  398.      * </pre>
  399.      *
  400.      * @param a the 5 coefficients
  401.      * @param y the real roots
  402.      * @return the number of real roots
  403.      */
  404.     private int realQuarticRoots(final double[] a, final double[] y) {
  405.         /* Treat the degenerate quartic as cubic */
  406.         if (Precision.equals(a[0], 0.)) {
  407.             final double[] aa = new double[a.length - 1];
  408.             System.arraycopy(a, 1, aa, 0, aa.length);
  409.             return realCubicRoots(aa, y);
  410.         }

  412.         // Transform coefficients
  413.         final double b  = a[1] / a[0];
  414.         final double c  = a[2] / a[0];
  415.         final double d  = a[3] / a[0];
  416.         final double e  = a[4] / a[0];
  417.         final double bh = b * 0.5;

  419.         // Solve resolvant cubic
  420.         final double[] z3 = new double[3];
  421.         final int i3 = realCubicRoots(new double[] {1.0, -c, b * d - 4. * e, e * (4. * c - b * b) - d * d}, z3);
  422.         if (i3 == 0) {
  423.             return 0;
  424.         }

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

  429.         // Compute auxiliary quantities
  430.         final double zh = z * 0.5;
  431.         final double p  = FastMath.max(z + bh * bh - c, 0.);
  432.         final double q  = FastMath.max(zh * zh - e, 0.);
  433.         final double r  = bh * z - d;
  434.         final double pp = FastMath.sqrt(p);
  435.         final double qq = FastMath.copySign(FastMath.sqrt(q), r);

  437.         // Solve quadratic factors of quartic equation
  438.         final double[] y1 = new double[2];
  439.         final int n1 = realQuadraticRoots(new double[] {1.0, bh - pp, zh - qq}, y1);
  440.         final double[] y2 = new double[2];
  441.         final int n2 = realQuadraticRoots(new double[] {1.0, bh + pp, zh + qq}, y2);

  443.         if (n1 == 2) {
  444.             if (n2 == 2) {
  445.                 y[0] = y1[0];
  446.                 y[1] = y1[1];
  447.                 y[2] = y2[0];
  448.                 y[3] = y2[1];
  449.                 return 4;
  450.             } else {
  451.                 y[0] = y1[0];
  452.                 y[1] = y1[1];
  453.                 return 2;
  454.             }
  455.         } else {
  456.             if (n2 == 2) {
  457.                 y[0] = y2[0];
  458.                 y[1] = y2[1];
  459.                 return 2;
  460.             } else {
  461.                 return 0;
  462.             }
  463.         }
  464.     }

  466.     /**
  467.      * Compute the real roots of a quartic equation.
  468.      *
  469.      * <pre>
  470.      * a[0] * x⁴ + a[1] * x³ + a[2] * x² + a[3] * x + a[4] = 0.
  471.      * </pre>
  472.      * @param <T> the type of the field elements
  473.      *
  474.      * @param a the 5 coefficients
  475.      * @param y the real roots
  476.      * @param field field of elements
  477.      * @return the number of real roots
  478.      */
  479.     private <T extends CalculusFieldElement<T>> int realQuarticRoots(final T[] a, final T[] y,
  480.                                                                  final Field<T> field) {

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

  484.         // Treat the degenerate quartic as cubic
  485.         if (Precision.equals(a[0].getReal(), 0.)) {
  486.             final T[] aa = MathArrays.buildArray(field, a.length - 1);
  487.             System.arraycopy(a, 1, aa, 0, aa.length);
  488.             return realCubicRoots(aa, y, field);
  489.         }

  491.         // Transform coefficients
  492.         final T b  = a[1].divide(a[0]);
  493.         final T c  = a[2].divide(a[0]);
  494.         final T d  = a[3].divide(a[0]);
  495.         final T e  = a[4].divide(a[0]);
  496.         final T bh = b.multiply(0.5);

  498.         // Solve resolvant cubic
  499.         final T[] z3 = MathArrays.buildArray(field, 3);
  500.         final T[] i = MathArrays.buildArray(field, 4);
  501.         i[0] = zero.add(1.0);
  502.         i[1] = c.negate();
  503.         i[2] = b.multiply(d).subtract(e.multiply(4.0));
  504.         i[3] = e.multiply(c.multiply(4.).subtract(b.multiply(b))).subtract(d.multiply(d));
  505.         final int i3 = realCubicRoots(i, z3, field);
  506.         if (i3 == 0) {
  507.             return 0;
  508.         }

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

  513.         // Compute auxiliary quantities
  514.         final T zh = z.multiply(0.5);
  515.         final T p  = FastMath.max(z.add(bh.multiply(bh)).subtract(c), zero);
  516.         final T q  = FastMath.max(zh.multiply(zh).subtract(e), zero);
  517.         final T r  = bh.multiply(z).subtract(d);
  518.         final T pp = FastMath.sqrt(p);
  519.         final T qq = FastMath.copySign(FastMath.sqrt(q), r);

  521.         // Solve quadratic factors of quartic equation
  522.         final T[] y1 = MathArrays.buildArray(field, 2);
  523.         final T[] n = MathArrays.buildArray(field, 3);
  524.         n[0] = zero.add(1.0);
  525.         n[1] = bh.subtract(pp);
  526.         n[2] = zh.subtract(qq);
  527.         final int n1 = realQuadraticRoots(n, y1);
  528.         final T[] y2 = MathArrays.buildArray(field, 2);
  529.         final T[] nn = MathArrays.buildArray(field, 3);
  530.         nn[0] = zero.add(1.0);
  531.         nn[1] = bh.add(pp);
  532.         nn[2] = zh.add(qq);
  533.         final int n2 = realQuadraticRoots(nn, y2);

  535.         if (n1 == 2) {
  536.             if (n2 == 2) {
  537.                 y[0] = y1[0];
  538.                 y[1] = y1[1];
  539.                 y[2] = y2[0];
  540.                 y[3] = y2[1];
  541.                 return 4;
  542.             } else {
  543.                 y[0] = y1[0];
  544.                 y[1] = y1[1];
  545.                 return 2;
  546.             }
  547.         } else {
  548.             if (n2 == 2) {
  549.                 y[0] = y2[0];
  550.                 y[1] = y2[1];
  551.                 return 2;
  552.             } else {
  553.                 return 0;
  554.             }
  555.         }
  556.     }

  558.     /**
  559.      * Compute the real roots of a cubic equation.
  560.      *
  561.      * <pre>
  562.      * a[0] * x³ + a[1] * x² + a[2] * x + a[3] = 0.
  563.      * </pre>
  564.      *
  565.      * @param a the 4 coefficients
  566.      * @param y the real roots sorted in descending order
  567.      * @return the number of real roots
  568.      */
  569.     private int realCubicRoots(final double[] a, final double[] y) {
  570.         if (Precision.equals(a[0], 0.)) {
  571.             // Treat the degenerate cubic as quadratic
  572.             final double[] aa = new double[a.length - 1];
  573.             System.arraycopy(a, 1, aa, 0, aa.length);
  574.             return realQuadraticRoots(aa, y);
  575.         }

  577.         // Transform coefficients
  578.         final double b  = -a[1] / (3. * a[0]);
  579.         final double c  =  a[2] / a[0];
  580.         final double d  =  a[3] / a[0];
  581.         final double b2 =  b * b;
  582.         final double p  =  b2 - c / 3.;
  583.         final double q  =  b * (b2 - c * 0.5) - d * 0.5;

  585.         // Compute discriminant
  586.         final double disc = p * p * p - q * q;

  588.         if (disc < 0.) {
  589.             // One real root
  590.             final double alpha  = q + FastMath.copySign(FastMath.sqrt(-disc), q);
  591.             final double cbrtAl = FastMath.cbrt(alpha);
  592.             final double cbrtBe = p / cbrtAl;

  594.             if (p < 0.) {
  595.                 y[0] = b + 2. * q / (cbrtAl * cbrtAl + cbrtBe * cbrtBe - p);
  596.             } else if (p > 0.) {
  597.                 y[0] = b + cbrtAl + cbrtBe;
  598.             } else {
  599.                 y[0] = b + cbrtAl;
  600.             }

  602.             return 1;

  604.         } else if (disc > 0.) {
  605.             // Three distinct real roots
  606.             final double phi = FastMath.atan2(FastMath.sqrt(disc), q) / 3.;
  607.             final double sqP = 2.0 * FastMath.sqrt(p);

  609.             y[0] = b + sqP * FastMath.cos(phi);
  610.             y[1] = b - sqP * FastMath.cos(FastMath.PI / 3. + phi);
  611.             y[2] = b - sqP * FastMath.cos(FastMath.PI / 3. - phi);

  613.             return 3;

  615.         } else {
  616.             // One distinct and two equals real roots
  617.             final double cbrtQ = FastMath.cbrt(q);
  618.             final double root1 = b + 2. * cbrtQ;
  619.             final double root2 = b - cbrtQ;
  620.             if (q < 0.) {
  621.                 y[0] = root2;
  622.                 y[1] = root2;
  623.                 y[2] = root1;
  624.             } else {
  625.                 y[0] = root1;
  626.                 y[1] = root2;
  627.                 y[2] = root2;
  628.             }

  630.             return 3;
  631.         }
  632.     }

  634.     /**
  635.      * Compute the real roots of a cubic equation.
  636.      *
  637.      * <pre>
  638.      * a[0] * x³ + a[1] * x² + a[2] * x + a[3] = 0.
  639.      * </pre>
  640.      *
  641.      * @param <T> the type of the field elements
  642.      * @param a the 4 coefficients
  643.      * @param y the real roots sorted in descending order
  644.      * @param field field of elements
  645.      * @return the number of real roots
  646.      */
  647.     private <T extends CalculusFieldElement<T>> int realCubicRoots(final T[] a, final T[] y,
  648.                                                                final Field<T> field) {

  650.         if (Precision.equals(a[0].getReal(), 0.)) {
  651.             // Treat the degenerate cubic as quadratic
  652.             final T[] aa = MathArrays.buildArray(field, a.length - 1);
  653.             System.arraycopy(a, 1, aa, 0, aa.length);
  654.             return realQuadraticRoots(aa, y);
  655.         }

  657.         // Transform coefficients
  658.         final T b  =  a[1].divide(a[0].multiply(3.)).negate();
  659.         final T c  =  a[2].divide(a[0]);
  660.         final T d  =  a[3].divide(a[0]);
  661.         final T b2 =  b.multiply(b);
  662.         final T p  =  b2.subtract(c.divide(3.));
  663.         final T q  =  b.multiply(b2.subtract(c.multiply(0.5))).subtract(d.multiply(0.5));

  665.         // Compute discriminant
  666.         final T disc = p.multiply(p).multiply(p).subtract(q.multiply(q));

  668.         if (disc.getReal() < 0.) {
  669.             // One real root
  670.             final T alpha  = FastMath.copySign(FastMath.sqrt(disc.negate()), q).add(q);
  671.             final T cbrtAl = FastMath.cbrt(alpha);
  672.             final T cbrtBe = p.divide(cbrtAl);

  674.             if (p .getReal() < 0.) {
  675.                 y[0] = q.divide(cbrtAl.multiply(cbrtAl).add(cbrtBe.multiply(cbrtBe)).subtract(p)).multiply(2.).add(b);
  676.             } else if (p.getReal() > 0.) {
  677.                 y[0] = b.add(cbrtAl).add(cbrtBe);
  678.             } else {
  679.                 y[0] = b.add(cbrtAl);
  680.             }

  682.             return 1;

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

  689.             y[0] = b.add(sqP.multiply(FastMath.cos(phi)));
  690.             y[1] = b.subtract(sqP.multiply(FastMath.cos(phi.add(b.getPi().divide(3.)))));
  691.             y[2] = b.subtract(sqP.multiply(FastMath.cos(phi.negate().add(b.getPi().divide(3.)))));

  693.             return 3;

  695.         } else {
  696.             // One distinct and two equals real roots
  697.             final T cbrtQ = FastMath.cbrt(q);
  698.             final T root1 = b.add(cbrtQ.multiply(2.0));
  699.             final T root2 = b.subtract(cbrtQ);
  700.             if (q.getReal() < 0.) {
  701.                 y[0] = root2;
  702.                 y[1] = root2;
  703.                 y[2] = root1;
  704.             } else {
  705.                 y[0] = root1;
  706.                 y[1] = root2;
  707.                 y[2] = root2;
  708.             }

  710.             return 3;
  711.         }
  712.     }

  714.     /**
  715.      * Compute the real roots of a quadratic equation.
  716.      *
  717.      * <pre>
  718.      * a[0] * x² + a[1] * x + a[2] = 0.
  719.      * </pre>
  720.      *
  721.      * @param a the 3 coefficients
  722.      * @param y the real roots sorted in descending order
  723.      * @return the number of real roots
  724.      */
  725.     private int realQuadraticRoots(final double[] a, final double[] y) {
  726.         if (Precision.equals(a[0], 0.)) {
  727.             // Degenerate quadratic
  728.             if (Precision.equals(a[1], 0.)) {
  729.                 // Degenerate linear equation: no real roots
  730.                 return 0;
  731.             }
  732.             // Linear equation: one real root
  733.             y[0] = -a[2] / a[1];
  734.             return 1;
  735.         }

  737.         // Transform coefficients
  738.         final double b = -0.5 * a[1] / a[0];
  739.         final double c =  a[2] / a[0];

  741.         // Compute discriminant
  742.         final double d =  b * b - c;

  744.         if (d < 0.) {
  745.             // No real roots
  746.             return 0;
  747.         } else if (d > 0.) {
  748.             // Two distinct real roots
  749.             final double y0 = b + FastMath.copySign(FastMath.sqrt(d), b);
  750.             final double y1 = c / y0;
  751.             y[0] = FastMath.max(y0, y1);
  752.             y[1] = FastMath.min(y0, y1);
  753.             return 2;
  754.         } else {
  755.             // Discriminant is zero: two equal real roots
  756.             y[0] = b;
  757.             y[1] = b;
  758.             return 2;
  759.         }
  760.     }

  762.     /**
  763.      * Compute the real roots of a quadratic equation.
  764.      *
  765.      * <pre>
  766.      * a[0] * x² + a[1] * x + a[2] = 0.
  767.      * </pre>
  768.      *
  769.      * @param <T> the type of the field elements
  770.      * @param a the 3 coefficients
  771.      * @param y the real roots sorted in descending order
  772.      * @return the number of real roots
  773.      */
  774.     private <T extends CalculusFieldElement<T>> int realQuadraticRoots(final T[] a, final T[] y) {

  776.         if (Precision.equals(a[0].getReal(), 0.)) {
  777.             // Degenerate quadratic
  778.             if (Precision.equals(a[1].getReal(), 0.)) {
  779.                 // Degenerate linear equation: no real roots
  780.                 return 0;
  781.             }
  782.             // Linear equation: one real root
  783.             y[0] = a[2].divide(a[1]).negate();
  784.             return 1;
  785.         }

  787.         // Transform coefficients
  788.         final T b = a[1].divide(a[0]).multiply(0.5).negate();
  789.         final T c =  a[2].divide(a[0]);

  791.         // Compute discriminant
  792.         final T d =  b.multiply(b).subtract(c);

  794.         if (d.getReal() < 0.) {
  795.             // No real roots
  796.             return 0;
  797.         } else if (d.getReal() > 0.) {
  798.             // Two distinct real roots
  799.             final T y0 = b.add(FastMath.copySign(FastMath.sqrt(d), b));
  800.             final T y1 = c.divide(y0);
  801.             y[0] = FastMath.max(y0, y1);
  802.             y[1] = FastMath.min(y0, y1);
  803.             return 2;
  804.         } else {
  805.             // Discriminant is zero: two equal real roots
  806.             y[0] = b;
  807.             y[1] = b;
  808.             return 2;
  809.         }
  810.     }

  812. }
Created with JaCoCo 0.8.11.202310140853