DSSTSolarRadiationPressure.java

  1. /* Copyright 2002-2018 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.MathUtils;
  22. import org.hipparchus.util.Precision;
  23. import org.orekit.errors.OrekitException;
  24. import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
  25. import org.orekit.forces.radiation.RadiationSensitive;
  26. import org.orekit.forces.radiation.SolarRadiationPressure;
  27. import org.orekit.propagation.SpacecraftState;
  28. import org.orekit.propagation.events.EventDetector;
  29. import org.orekit.utils.ExtendedPVCoordinatesProvider;
  30. import org.orekit.utils.PVCoordinatesProvider;

  32. /** Solar radiation pressure contribution to the
  33.  *  {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
  34.  *  <p>
  35.  *  The solar radiation pressure acceleration is computed through the
  36.  *  acceleration model of
  37.  *  {@link org.orekit.forces.radiation.SolarRadiationPressure SolarRadiationPressure}.
  38.  *  </p>
  39.  *
  40.  *  @author Pascal Parraud
  41.  */
  42. public class DSSTSolarRadiationPressure extends AbstractGaussianContribution {

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

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

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

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

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

  59.     /** Sun model. */
  60.     private final PVCoordinatesProvider sun;

  62.     /** Central Body radius. */
  63.     private final double                ae;

  65.     /** Spacecraft model for radiation acceleration computation. */
  66.     private final RadiationSensitive spacecraft;


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

  96.     /**
  97.      * Simple constructor with default reference values and spherical spacecraft.
  98.      * <p>
  99.      * When this constructor is used, the reference values are:
  100.      * </p>
  101.      * <ul>
  102.      * <li>d<sub>ref</sub> = 149597870000.0 m</li>
  103.      * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  104.      * </ul>
  105.      * <p>
  106.      * The spacecraft has a spherical shape.
  107.      * </p>
  108.      *
  109.      * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
  110.      * @param area cross sectional area of satellite
  111.      * @param sun Sun model
  112.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  113.      * @since 9.2
  114.      */
  115.     public DSSTSolarRadiationPressure(final double cr, final double area,
  116.                                       final ExtendedPVCoordinatesProvider sun,
  117.                                       final double equatorialRadius) {
  118.         this(D_REF, P_REF, cr, area, sun, equatorialRadius);
  119.     }

  121.     /**
  122.      * Simple constructor with default reference values, but custom spacecraft.
  123.      * <p>
  124.      * When this constructor is used, the reference values are:
  125.      * </p>
  126.      * <ul>
  127.      * <li>d<sub>ref</sub> = 149597870000.0 m</li>
  128.      * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  129.      * </ul>
  130.      *
  131.      * @param sun Sun model
  132.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  133.      * @param spacecraft spacecraft model
  134.      * @deprecated as of 9.2 replaced by {{@link #DSSTSolarRadiationPressure(ExtendedPVCoordinatesProvider,
  135.      * double, RadiationSensitive)}
  136.      */
  137.     public DSSTSolarRadiationPressure(final PVCoordinatesProvider sun,
  138.                                       final double equatorialRadius,
  139.                                       final RadiationSensitive spacecraft) {
  140.         this(D_REF, P_REF, sun, equatorialRadius, spacecraft);
  141.     }

  143.     /**
  144.      * Simple constructor with default reference values, but custom spacecraft.
  145.      * <p>
  146.      * When this constructor is used, the reference values are:
  147.      * </p>
  148.      * <ul>
  149.      * <li>d<sub>ref</sub> = 149597870000.0 m</li>
  150.      * <li>p<sub>ref</sub> = 4.56 10<sup>-6</sup> N/m²</li>
  151.      * </ul>
  152.      *
  153.      * @param sun Sun model
  154.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  155.      * @param spacecraft spacecraft model
  156.      * @since 9.2
  157.      */
  158.     public DSSTSolarRadiationPressure(final ExtendedPVCoordinatesProvider sun,
  159.                                       final double equatorialRadius,
  160.                                       final RadiationSensitive spacecraft) {
  161.         this(D_REF, P_REF, sun, equatorialRadius, spacecraft);
  162.     }

  164.     /**
  165.      * Constructor with customizable reference values but spherical spacecraft.
  166.      * <p>
  167.      * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
  168.      * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
  169.      * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  170.      * N/m² solar radiation pressure.
  171.      * </p>
  172.      *
  173.      * @param dRef reference distance for the solar radiation pressure (m)
  174.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  175.      * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
  176.      * @param area cross sectional area of satellite
  177.      * @param sun Sun model
  178.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  179.      * @deprecated as of 9.2 replaced by {@link #DSSTSolarRadiationPressure(double, double,
  180.      * double, double, ExtendedPVCoordinatesProvider, double)
  181.      */
  182.     @Deprecated
  183.     public DSSTSolarRadiationPressure(final double dRef, final double pRef,
  184.                                       final double cr, final double area,
  185.                                       final PVCoordinatesProvider sun,
  186.                                       final double equatorialRadius) {

  188.         // cR being the DSST SRP coef and assuming a spherical spacecraft,
  189.         // the conversion is:
  190.         // cR = 1 + (1 - kA) * (1 - kR) * 4 / 9
  191.         // with kA arbitrary sets to 0
  192.         this(dRef, pRef, sun, equatorialRadius, new IsotropicRadiationSingleCoefficient(area, cr));
  193.     }

  195.     /**
  196.      * Constructor with customizable reference values but spherical spacecraft.
  197.      * <p>
  198.      * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
  199.      * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
  200.      * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  201.      * N/m² solar radiation pressure.
  202.      * </p>
  203.      *
  204.      * @param dRef reference distance for the solar radiation pressure (m)
  205.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  206.      * @param cr satellite radiation pressure coefficient (assuming total specular reflection)
  207.      * @param area cross sectional area of satellite
  208.      * @param sun Sun model
  209.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  210.      * @since 9.2
  211.      */
  212.     public DSSTSolarRadiationPressure(final double dRef, final double pRef,
  213.                                       final double cr, final double area,
  214.                                       final ExtendedPVCoordinatesProvider sun,
  215.                                       final double equatorialRadius) {

  217.         // cR being the DSST SRP coef and assuming a spherical spacecraft,
  218.         // the conversion is:
  219.         // cR = 1 + (1 - kA) * (1 - kR) * 4 / 9
  220.         // with kA arbitrary sets to 0
  221.         this(dRef, pRef, sun, equatorialRadius, new IsotropicRadiationSingleCoefficient(area, cr));
  222.     }

  224.     /**
  225.      * Complete constructor.
  226.      * <p>
  227.      * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
  228.      * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
  229.      * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  230.      * N/m² solar radiation pressure.
  231.      * </p>
  232.      *
  233.      * @param dRef reference distance for the solar radiation pressure (m)
  234.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  235.      * @param sun Sun model
  236.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  237.      * @param spacecraft spacecraft model
  238.      * @deprecated as of 9.2 replaced by {@link #DSSTSolarRadiationPressure(double, double,
  239.      * ExtendedPVCoordinatesProvider, double, RadiationSensitive)
  240.      */
  241.     @Deprecated
  242.     public DSSTSolarRadiationPressure(final double dRef, final double pRef,
  243.                                       final PVCoordinatesProvider sun, final double equatorialRadius,
  244.                                       final RadiationSensitive spacecraft) {

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

  250.         this.sun  = sun;
  251.         this.ae   = equatorialRadius;
  252.         this.spacecraft = spacecraft;
  253.     }

  255.     /**
  256.      * Complete constructor.
  257.      * <p>
  258.      * Note that reference solar radiation pressure <code>pRef</code> in N/m² is linked
  259.      * to solar flux SF in W/m² using formula pRef = SF/c where c is the speed of light
  260.      * (299792458 m/s). So at 1UA a 1367 W/m² solar flux is a 4.56 10<sup>-6</sup>
  261.      * N/m² solar radiation pressure.
  262.      * </p>
  263.      *
  264.      * @param dRef reference distance for the solar radiation pressure (m)
  265.      * @param pRef reference solar radiation pressure at dRef (N/m²)
  266.      * @param sun Sun model
  267.      * @param equatorialRadius central body equatorial radius (for shadow computation)
  268.      * @param spacecraft spacecraft model
  269.      * @since 9.2
  270.      */
  271.     public DSSTSolarRadiationPressure(final double dRef, final double pRef,
  272.                                       final ExtendedPVCoordinatesProvider sun,
  273.                                       final double equatorialRadius,
  274.                                       final RadiationSensitive spacecraft) {

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

  280.         this.sun  = sun;
  281.         this.ae   = equatorialRadius;
  282.         this.spacecraft = spacecraft;
  283.     }

  285.     /** Get spacecraft shape.
  286.      * @return the spacecraft shape.
  287.      */
  288.     public RadiationSensitive getSpacecraft() {
  289.         return spacecraft;
  290.     }

  292.     /** {@inheritDoc} */
  293.     public EventDetector[] getEventsDetectors() {
  294.         return null;
  295.     }

  297.     /** {@inheritDoc} */
  298.     protected double[] getLLimits(final SpacecraftState state) throws OrekitException {
  299.         // Default bounds without shadow [-PI, PI]
  300.         final double[] ll = {-FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0),
  301.                              FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0)};

  303.         // Direction cosines of the Sun in the equinoctial frame
  304.         final Vector3D sunDir = sun.getPVCoordinates(state.getDate(), state.getFrame()).getPosition().normalize();
  305.         final double alpha = sunDir.dotProduct(f);
  306.         final double beta  = sunDir.dotProduct(g);
  307.         final double gamma = sunDir.dotProduct(w);

  309.         // Compute limits only if the perigee is close enough from the central body to be in the shadow
  310.         if (FastMath.abs(gamma * a * (1. - ecc)) < ae) {

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

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

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

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

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

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

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

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

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

  465.     /**
  466.      * Compute the real roots of a cubic equation.
  467.      *
  468.      * <pre>
  469.      * a[0] * x³ + a[1] * x² + a[2] * x + a[3] = 0.
  470.      * </pre>
  471.      *
  472.      * @param a the 4 coefficients
  473.      * @param y the real roots sorted in descending order
  474.      * @return the number of real roots
  475.      */
  476.     private int realCubicRoots(final double[] a, final double[] y) {
  477.         if (Precision.equals(a[0], 0.)) {
  478.             // Treat the degenerate cubic as quadratic
  479.             final double[] aa = new double[a.length - 1];
  480.             System.arraycopy(a, 1, aa, 0, aa.length);
  481.             return realQuadraticRoots(aa, y);
  482.         }

  484.         // Transform coefficients
  485.         final double b  = -a[1] / (3. * a[0]);
  486.         final double c  =  a[2] / a[0];
  487.         final double d  =  a[3] / a[0];
  488.         final double b2 =  b * b;
  489.         final double p  =  b2 - c / 3.;
  490.         final double q  =  b * (b2 - c * 0.5) - d * 0.5;

  492.         // Compute discriminant
  493.         final double disc = p * p * p - q * q;

  495.         if (disc < 0.) {
  496.             // One real root
  497.             final double alpha  = q + FastMath.copySign(FastMath.sqrt(-disc), q);
  498.             final double cbrtAl = FastMath.cbrt(alpha);
  499.             final double cbrtBe = p / cbrtAl;

  501.             if (p < 0.) {
  502.                 y[0] = b + 2. * q / (cbrtAl * cbrtAl + cbrtBe * cbrtBe - p);
  503.             } else if (p > 0.) {
  504.                 y[0] = b + cbrtAl + cbrtBe;
  505.             } else {
  506.                 y[0] = b + cbrtAl;
  507.             }

  509.             return 1;

  511.         } else if (disc > 0.) {
  512.             // Three distinct real roots
  513.             final double phi = FastMath.atan2(FastMath.sqrt(disc), q) / 3.;
  514.             final double sqP = 2.0 * FastMath.sqrt(p);

  516.             y[0] = b + sqP * FastMath.cos(phi);
  517.             y[1] = b - sqP * FastMath.cos(FastMath.PI / 3. + phi);
  518.             y[2] = b - sqP * FastMath.cos(FastMath.PI / 3. - phi);

  520.             return 3;

  522.         } else {
  523.             // One distinct and two equals real roots
  524.             final double cbrtQ = FastMath.cbrt(q);
  525.             final double root1 = b + 2. * cbrtQ;
  526.             final double root2 = b - cbrtQ;
  527.             if (q < 0.) {
  528.                 y[0] = root2;
  529.                 y[1] = root2;
  530.                 y[2] = root1;
  531.             } else {
  532.                 y[0] = root1;
  533.                 y[1] = root2;
  534.                 y[2] = root2;
  535.             }

  537.             return 3;
  538.         }
  539.     }

  541.     /**
  542.      * Compute the real roots of a quadratic equation.
  543.      *
  544.      * <pre>
  545.      * a[0] * x² + a[1] * x + a[2] = 0.
  546.      * </pre>
  547.      *
  548.      * @param a the 3 coefficients
  549.      * @param y the real roots sorted in descending order
  550.      * @return the number of real roots
  551.      */
  552.     private int realQuadraticRoots(final double[] a, final double[] y) {
  553.         if (Precision.equals(a[0], 0.)) {
  554.             // Degenerate quadratic
  555.             if (Precision.equals(a[1], 0.)) {
  556.                 // Degenerate linear equation: no real roots
  557.                 return 0;
  558.             }
  559.             // Linear equation: one real root
  560.             y[0] = -a[2] / a[1];
  561.             return 1;
  562.         }

  564.         // Transform coefficients
  565.         final double b = -0.5 * a[1] / a[0];
  566.         final double c =  a[2] / a[0];

  568.         // Compute discriminant
  569.         final double d =  b * b - c;

  571.         if (d < 0.) {
  572.             // No real roots
  573.             return 0;
  574.         } else if (d > 0.) {
  575.             // Two distinct real roots
  576.             final double y0 = b + FastMath.copySign(FastMath.sqrt(d), b);
  577.             final double y1 = c / y0;
  578.             y[0] = FastMath.max(y0, y1);
  579.             y[1] = FastMath.min(y0, y1);
  580.             return 2;
  581.         } else {
  582.             // Discriminant is zero: two equal real roots
  583.             y[0] = b;
  584.             y[1] = b;
  585.             return 2;
  586.         }
  587.     }
  588. }
Created with JaCoCo 0.7.9.201702052155