PRO solve_realroots_cubic, a, b, c, d, x1, x2, x3, nroots, status
; Description: This module solves for the real root(s) of a cubic equation, of which there
; is always at least one. A cubic equation has the form:
;
; a*x^3 + b*x^2 + c*x + d = 0
;
; where a, b, c and d are the polynomial coefficients, which also happen to be
; the input parameters "a", "b", "c" and "d" respectively. The module returns
; the distinct real roots "x1", "x2" and "x3" in ascending order, and the
; number "nroots" of distinct real roots of the equation. If there is only one
; distinct real root, then "x2" and "x3" are set to "0.0", and if there are
; only two distinct real roots, then "x3" is set to "0.0".
; To allow many similar equations to be solved using a single call to this
; module, any one of the input parameters "a", "b", "c" and "d" may be supplied
; as a vector rather than as a scalar. The remaining input parameters must be
; supplied as scalars. In this case, the output parameters are returned as
; vectors of the same length as the only vector input parameter.
;
; Input Parameters:
;
; a - FLOAT/DOUBLE SCALAR/VECTOR - The polynomial coefficient of the cubic term in the
; cubic equation. This parameter must be non-zero (in
; order for the equation to be a cubic equation). If
; this parameter is a vector, then all of the values must
; be non-zero.
; b - FLOAT/DOUBLE SCALAR/VECTOR - The polynomial coefficient of the quadratic term in the
; cubic equation.
; c - FLOAT/DOUBLE SCALAR/VECTOR - The polynomial coefficient of the linear term in the
; cubic equation.
; d - FLOAT/DOUBLE SCALAR/VECTOR - The constant term in the cubic equation.
;
; Output Parameters:
;
; x1 - DOUBLE SCALAR/VECTOR - The smallest distinct real root of the cubic equation.
; x2 - DOUBLE SCALAR/VECTOR - The next smallest distinct real root of the cubic equation.
; If the equation only has one distinct real root, then "x2"
; is assigned a value of "0.0".
; x3 - DOUBLE SCALAR/VECTOR - The largest distinct real root of the cubic equation. If
; the equation only has one or two distinct real roots, then
; "x3" is assigned a value of "0.0".
; nroots - INTEGER SCALAR/VECTOR - The number of distinct real roots of the cubic equation
; which may take one of the values "1", "2" or "3".
; status - INTEGER - If the module successfully determines the real root(s) of the cubic
; equation(s), then "status" is returned with a value of "1", otherwise
; it is returned with a value of "0".
;
; N.B: (i) If all of the input parameters are scalars, then all of the output parameters are
; scalars. If one of the input parameters is a vector, then all of the output
; parameters are vectors of the same length as the input vector. No more than ONE of
; the input parameters is allowed to be supplied as a vector.
; (ii) This implementation of the solution for the root(s) of a cubic equation correctly
; identifies when a cubic equation with three real roots has repeated roots, unlike
; the method described in Section 5.6 of "Numerical Recipes In Fortran 77: The Art
; Of Scientific Computing".
; (iii) No attempt has been made to identify and avoid problems with "catastrophic
; cancellation", when substantial loss of precision for the values of the roots
; can occur.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 21/10/2009 - Module created (dmb).
;Set the default output parameter values
;Determine the dimensions of the input parameters
;Check that at most only one of the input parameters is a vector
;If all of the input parameters are scalars
;Check that the input parameters are of the correct number type
;Check that the polynomial coefficient of the cubic term in the cubic equation is not zero
;Convert the cubic equation to a monic cubic equation by dividing all the polynomial
;coefficients by the coefficient of the cubic term
;Calculate the discriminant of the cubic equation
;If the cubic equation has only one real root (and two imaginary roots), then calculate
;the value of the real root
;Calculate the value of the real root
;If the cubic equation has three distinct real roots, then calculate the values of the
;roots, and sort them into ascending order
;Calculate the values of the three distinct real roots
;Sort the roots into ascending order
;If the cubic equation has three real roots with at least two of them being equal, then
;calculate the values of the roots, and sort them into ascending order
;If the cubic equation has two distinct real roots
;Calculate the values of the two distinct real roots
;If the cubic equation has three equal real roots
;Set "status" to "1"
;Return with the roots of the cubic equation
;If the polynomial coefficient of the cubic term is a vector
;Check that "a" is of the correct number type, and determine its length
;Check that the remaining parameters are scalars of the correct number type
;Check that "a" has no zero values
;Convert the cubic equations to monic cubic equations by dividing all the polynomial
;coefficients by the coefficient of the cubic term
;If the polynomial coefficient of the quadratic term is a vector
;Check that "b" is of the correct number type, and determine its length
;Check that the remaining parameters are scalars of the correct number type
;Check that "a" is not zero
;Convert the cubic equations to monic cubic equations by dividing all the polynomial
;coefficients by the coefficient of the cubic term
;If the polynomial coefficient of the linear term is a vector
;Check that "c" is of the correct number type, and determine its length
;Check that the remaining parameters are scalars of the correct number type
;Check that "a" is not zero
;Convert the cubic equations to monic cubic equations by dividing all the polynomial
;coefficients by the coefficient of the cubic term
;If the constant term is a vector
;Check that "d" is of the correct number type, and determine its length
;Check that the remaining parameters are scalars of the correct number type
;Check that "a" is not zero
;Convert the cubic equations to monic cubic equations by dividing all the polynomial
;coefficients by the coefficient of the cubic term
;Set up the output parameters
;Calculate the discriminant of the cubic equation
;Determine where the cubic equation(s) have only one real root (and two imaginary roots),
;where they have three distinct real roots, where they have three real roots with two of
;the roots being equal, and where they have three real roots with all three roots being
;equal.
;If there is at least one cubic equation with only one real root (and two imaginary roots),
;then for each of the relevant cubic equations calculate the value of the real root
;Calculate the value of the real root
;If there is at least one cubic equation with three distinct real roots, then for each of
;the relevant cubic equations calculate the values of the real roots, and sort them into
;ascending order
;Calculate the values of the three distinct real roots
;Sort the roots into ascending order
;If there is at least one cubic equation with three real roots where two of the roots are
;equal, then for each of the relevant cubic equations calculate the values of the distinct
;real roots, and sort them into ascending order
;Calculate the values of the two distinct real roots and sort the roots into ascending
;order
;If there is at least one cubic equation with three real roots where all three roots are
;equal, then for each of the relevant cubic equations calculate the value of the distinct
;real root
;Calculate the value of the distinct real root
;Set "status" to "1"