PRO solve_realroots_cubic, a, b, c, d, x1, x2, x3, nroots ; 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. ; If the module fails, then "nroots" is set to "-1". ; If the module is successful, then "nroots" will take one ; of the values "1", "2" or "3". ; ; 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/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 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 ;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 ;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 ;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 ;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 ;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