DSSTSolarRadiationPressure.java

  1. /* Copyright 2002-2025 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.ExtendedPositionProvider;
  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 ExtendedPositionProvider 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 ExtendedPositionProvider 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 ExtendedPositionProvider 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 ExtendedPositionProvider 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 ExtendedPositionProvider 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.     @Override
  198.     protected List<ParameterDriver> getParametersDriversWithoutMu() {
  199.         return spacecraft.getRadiationParametersDrivers();
  200.     }

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

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

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

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

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

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

  297.         final Field<T> field = state.getDate().getField();
  298.         final T zero = field.getZero();
  299.         final T one  = field.getOne();
  300.         final T pi   = one.getPi();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  588.         // Compute discriminant
  589.         final double disc = p * p * p - q * q;

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

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

  605.             return 1;

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

  612.             y[0] = b + sqP * FastMath.cos(phi);
  613.             y[1] = b - sqP * FastMath.cos(FastMath.PI / 3. + phi);
  614.             y[2] = b - sqP * FastMath.cos(FastMath.PI / 3. - phi);

  616.             return 3;

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

  633.             return 3;
  634.         }
  635.     }

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

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

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

  668.         // Compute discriminant
  669.         final T disc = p.multiply(p).multiply(p).subtract(q.multiply(q));

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

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

  685.             return 1;

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

  692.             y[0] = b.add(sqP.multiply(FastMath.cos(phi)));
  693.             y[1] = b.subtract(sqP.multiply(FastMath.cos(phi.add(b.getPi().divide(3.)))));
  694.             y[2] = b.subtract(sqP.multiply(FastMath.cos(phi.negate().add(b.getPi().divide(3.)))));

  696.             return 3;

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

  713.             return 3;
  714.         }
  715.     }

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

  740.         // Transform coefficients
  741.         final double b = -0.5 * a[1] / a[0];
  742.         final double c =  a[2] / a[0];

  744.         // Compute discriminant
  745.         final double d =  b * b - c;

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

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

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

  790.         // Transform coefficients
  791.         final T b = a[1].divide(a[0]).multiply(0.5).negate();
  792.         final T c =  a[2].divide(a[0]);

  794.         // Compute discriminant
  795.         final T d =  b.square().subtract(c);

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

  815. }
Created with JaCoCo 0.8.12.202403310830