PRO solve_realroots_quadratic, a, b, c, x1, x2, nroots

; Description: This module solves for the real root(s) of a quadratic equation, of which there
;              may be zero, one or two. A quadratic equation has the form:
;
;              a*x^2 + b*x + c = 0
;
;              where a, b and c are the polynomial coefficients, which also happen to be the
;              input parameters "a", "b" and "c" respectively. The module returns the distinct 
;              real roots "x1" and "x2" in ascending order, and the number "nroots" of distinct
;              real roots of the equation. If there are no real roots, then "x1" and "x2" are 
;              set to "0.0", and if there is only one distinct real root, then "x2" is set to 
;              "0.0".
;                Note that where at least one real root exists, the module makes use of the 
;              quadratic formula to calculate the real root with the greatest magnitude, and 
;              then employs a useful relationship between the real roots to obtain the other 
;              root. The quadratic formula is:
;
;              x = ( -b +- sqrt( b^2 - (4*a*c) )) / (2 * a)
;
;              The useful relationship is:
;
;              x2 = c / (a * x1)
;
;              This approach avoids the risk of "catastrophic cancellation", where substantial 
;              loss of precision for one of the roots can occur when the magnitudes of the 
;              terms "b" and "sqrt( b^2 - (4*a*c) )" in the quadratic equation are similar.
;                To allow many similar equations to be solved using a single call to this
;              module, any one of the input parameters "a", "b" and "c" 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 quadratic term in the
;                                    quadratic equation. This parameter must be non-zero (in
;                                    order for the equation to be a quadratic 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 linear term in the
;                                    quadratic equation.
;   c - FLOAT/DOUBLE SCALAR/VECTOR - The constant term in the quadratic equation.
;
; Output Parameters:
;
;   x1 - DOUBLE SCALAR/VECTOR - The smallest distinct real root of the quadratic equation. If
;                               the equation does not have any real roots, then "x1" is
;                               assigned a value of "0.0".
;   x2 - DOUBLE SCALAR/VECTOR - The largest distinct real root of the quadratic equation. If
;                               the equation only has one distinct real root, or no real roots,
;                               then "x2" is assigned a value of "0.0".
;   nroots - INTEGER SCALAR/VECTOR - The number of distinct real roots of the quadratic equation.
;                                    If the module fails, then "nroots" is set to "-1". If the 
;                                    module is successful, then "nroots" will take one of the 
;                                    values "0", "1" or "2".
;
; N.B: 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.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
;   27/04/2009 - Module created (dmb)

;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 quadratic term in the quadratic equation is not
  ;zero

  ;Convert the quadratic equation to a monic quadratic equation by dividing all the polynomial
  ;coefficients by the coefficient of the quadratic term

  ;Calculate the discriminant of the quadratic equation

  ;If the quadratic equation does not have any real roots, then set "nroots" to "0"

  ;If the quadratic equation has two distinct real roots, then calculate the values of the real
  ;roots, and sort them into ascending order

    ;If "b" is non-zero, then use the quadratic formula to calculate the real root with the
    ;greatest magnitude, and subsequently calculate the value of the remaining root using the 
    ;useful relationship between the roots.

      ;Calculate the values of the two distinct real roots

      ;Sort the roots into ascending order

    ;If "b" is zero, then use the quadratic formula to calculate the values of both real roots

      ;Calculate the values of the two distinct real roots

  ;If the quadratic equation has one distinct real root, then calculate the value of the real
  ;root

    ;Calculate the value of the one distinct real root

  ;Return with the roots of the quadratic equation

;If the polynomial coefficient of the quadratic term is a vector

  ;Check that "a" is of the correct number type, and determine its length

  ;Check that the remaining parameters are scalars

  ;Check that "a" has no zero values

  ;Convert the quadratic equations to monic quadratic equations by dividing all the polynomial
  ;coefficients by the coefficient of the quadratic term

;If the polynomial coefficient of the linear term is a vector

  ;Check that "b" is of the correct number type, and determine its length

  ;Check that the remaining parameters are scalars

  ;Check that "a" is not zero

  ;Convert the quadratic equations to monic quadratic equations by dividing all the polynomial
  ;coefficients by the coefficient of the quadratic term

;If the constant term is a vector

  ;Check that "c" is of the correct number type, and determine its length

  ;Check that the remaining parameters are scalars

  ;Check that "a" is not zero

  ;Convert the quadratic equations to monic quadratic equations by dividing all the polynomial
  ;coefficients by the coefficient of the quadratic term

;Set up the output parameters

;Calculate the discriminant of the quadratic equation

;Determine where the quadratic equation(s) have one distinct real root, and where they 
;have two distinct real roots. In the cases where the quadratic equation(s) have two 
;distinct real roots, distinguish between the cases where "b" is non-zero and where "b" 
;is equal to zero.

;If there is at least one quadratic equation with one distinct real root, then for each
;of the relevant quadratic equations calculate the value of the real root

  ;Calculate the value of the real root

;If there is at least one quadratic equation with two distinct real roots and with "b"
;non-zero, then for each of the relevant quadratic equations calculate the values of the 
;real roots, and sort them into ascending order.

  ;Use the quadratic formula to calculate the real root with the greatest magnitude, and
  ;subsequently calculate the value of the remaining root using the useful relationship 
  ;between the roots.

  ;Sort the roots into ascending order

;If there is at least one quadratic equation with two distinct real roots and with "b"
;equal to zero, then for each of the relevant quadratic equations calculate the values 
;of the real roots.

  ;Calculate the values of the two distinct real roots