// Source file for the Oxbridge Stransverse Mass Library -- oxbridgekinetics. // See http://www.hep.phy.cam.ac.uk/~lester/mt2/index.html // Authors: Christopher Lester and Alan Barr #include namespace Mt2 { namespace quartic { //****************************************************************************** std::complex c8_sqrt (const std::complex x ); double c8_argument (const std::complex & x); double c8_magnitude (const std::complex & x); double r8_sign ( double x ); void r8poly3_root (const double a, const double b, const double c, const double d, std::complex & r1, std::complex & r2, std::complex & r3 ); void r8poly4_root (const double a, const double b, const double c, const double d, const double e, std::complex & r1, std::complex & r2, std::complex & r3, std::complex & r4 ) //****************************************************************************** // // Purpose: // // R8POLY4_ROOT returns the four roots of a quartic polynomial. // // Discussion: // // The polynomial has the form: // // A * X**4 + B * X**3 + C * X**2 + D * X + E = 0 // // Modified: // // 27 October 2005 // // Parameters: // // Input, double A, B, C, D, the coefficients of the polynomial. // A must not be zero. // // Output, complex *R1, *R2, *R3, *R4, the roots of the polynomial. // { double a3; double a4; double b3; double b4; double c3; double c4; double d3; double d4; std::complex p; std::complex q; std::complex r; const std::complex zero(0,0); if ( a == 0.0 ) { std::cout << "\n"; std::cout << "R8POLY4_ROOT - Fatal error!\n"; std::cout << " A must not be zero!\n"; throw "R8POLY4_ROOT - Fatal error! A must not be zero!"; } a4 = b / a; b4 = c / a; c4 = d / a; d4 = e / a; // // Set the coefficients of the resolvent cubic equation. // a3 = 1.0; b3 = -b4; c3 = a4 * c4 - 4.0 * d4; d3 = -a4 * a4 * d4 + 4.0 * b4 * d4 - c4 * c4; // // Find the roots of the resolvent cubic. // r8poly3_root ( a3, b3, c3, d3, r1, r2, r3 ); // // Choose one root of the cubic, here R1. // // Set R = sqrt ( 0.25 * A4**2 - B4 + R1 ) // // if r1 is real, and if r1 + 0.25 * a4 * a4 - b4 is positive, then r that is about to follow will be exactly real. r = c8_sqrt ( r1 + 0.25 * a4 * a4 - b4 ); // if ( r != zero ) if ((r.real()!=0.0) || (r.imag()!=0.0)) { p = c8_sqrt ( zero + 0.75 * a4 * a4 - r * r - 2.0 * b4 + 0.25 * ( 4.0 * a4 * b4 - 8.0 * c4 - a4 * a4 * a4 ) / r ); q = c8_sqrt ( zero + 0.75 * a4 * a4 - r * r - 2.0 * b4 - 0.25 * ( 4.0 * a4 * b4 - 8.0 * c4 - a4 * a4 * a4 ) / r ); } else { p = c8_sqrt ( zero + 0.75 * a4 * a4 - 2.0 * b4 + 2.0 * c8_sqrt ( r1 * r1 - 4.0 * d4 ) ); q = c8_sqrt ( zero + 0.75 * a4 * a4 - 2.0 * b4 - 2.0 * c8_sqrt ( r1 * r1 - 4.0 * d4 ) ); } // // Set the roots. // r1 = zero -0.25 * a4 + 0.5 * r + 0.5 * p; r2 = zero -0.25 * a4 + 0.5 * r - 0.5 * p; r3 = zero -0.25 * a4 - 0.5 * r + 0.5 * q; r4 = zero -0.25 * a4 - 0.5 * r - 0.5 * q; return; } //****************************************************************************** void r8poly3_root (const double a, const double b, const double c, const double d, std::complex & r1, std::complex & r2, std::complex & r3 ) //****************************************************************************** // // Purpose: // // R8POLY3_ROOT returns the three roots of a cubic polynomial. // // Discussion: // // The polynomial has the form // // A * X**3 + B * X**2 + C * X + D = 0 // // Modified: // // 25 October 2005 // // Parameters: // // Input, double A, B, C, D, the coefficients of the polynomial. // A must not be zero. // // Output, complex *R1, *R2, *R3, the roots of the polynomial, which // will include at least one real root. // { std::complex i; double pi = 3.141592653589793; double q; double r; double s1; double s2; double temp; double theta; const std::complex zero(0,0); if ( a == 0.0 ) { std::cout << "\n"; std::cout << "R8POLY3_ROOT - Fatal error!\n"; std::cout << " A must not be zero!\n"; throw "R8POLY3_ROOT - Fatal error! A must not be zero!"; } i = std::complex ( 0.0, 1.0 ); q = ( pow ( b / a, 2 ) - 3.0 * ( c / a ) ) / 9.0; r = ( 2.0 * pow ( b / a, 3 ) - 9.0 * ( b / a ) * ( c / a ) + 27.0 * ( d / a ) ) / 54.0; if ( r * r < q * q * q ) { theta = acos ( r / sqrt ( pow ( q, 3 ) ) ); r1 = -2.0 * sqrt ( q ) * cos ( theta / 3.0 ); r2 = -2.0 * sqrt ( q ) * cos ( ( theta + 2.0 * pi ) / 3.0 ); r3 = -2.0 * sqrt ( q ) * cos ( ( theta + 4.0 * pi ) / 3.0 ); } else if ( q * q * q <= r * r ) { temp = -r + sqrt ( r * r - q * q * q ); s1 = r8_sign ( temp ) * pow ( fabs ( temp ), 1.0 / 3.0 ); temp = -r - sqrt ( r * r - q * q * q ); s2 = r8_sign ( temp ) * pow ( fabs ( temp ), 1.0 / 3.0 ); r1 = s1 + s2; r2 = zero -0.5 * ( s1 + s2 ) + i * 0.5 * sqrt ( 3.0 ) * ( s1 - s2 ); r3 = zero -0.5 * ( s1 + s2 ) - i * 0.5 * sqrt ( 3.0 ) * ( s1 - s2 ); } r1 = r1 - b / ( 3.0 * a ); r2 = r2 - b / ( 3.0 * a ); r3 = r3 - b / ( 3.0 * a ); return; } //****************************************************************************** std::complex c8_sqrt (const std::complex x ) //****************************************************************************** // // Purpose: // // C8_SQRT returns the principal square root of a C8. // // Discussion: // // A C8 is a double precision complex value. // // Modified: // // 27 October 2005 // // Author: // // John Burkardt // // Parameters: // // Input, complex X, the number whose square root is desired. // // Output, complex C8_SQRT, the square root of X. // { double argument; double magnitude; std::complex value; argument = c8_argument ( x ); magnitude = c8_magnitude ( x ); if ( magnitude == 0.0 ) { value = std::complex ( 0.0, 0.0 ); } else { value = pow ( magnitude, 1.0 / 2.0 ) * std::complex ( cos ( argument / 2.0 ), sin ( argument / 2.0 ) ); } return value; } //********************************************************************* double r8_sign ( double x ) //********************************************************************* // // Purpose: // // R8_SIGN returns the sign of a double precision number. // // Modified: // // 18 October 2004 // // Author: // // John Burkardt // // Parameters: // // Input, double X, the number whose sign is desired. // // Output, double R8_SIGN, the sign of X. // { if ( x < 0.0 ) { return ( -1.0 ); } else { return ( 1.0 ); } } //****************************************************************************** double c8_argument (const std::complex & x) //****************************************************************************** // // Purpose: // // C8_ARGUMENT returns the argument of a C8. // // Discussion: // // A C8 is a double precision std::complex value. // // Modified: // // 21 October 2005 // // Author: // // John Burkardt // // Parameters: // // Input, complex X, the value whose argument is desired. // // Output, double C8_ARGUMENT, the argument of X. // { double argument; argument = atan2(x.imag(),x.real()); return argument; } //****************************************************************************** double c8_magnitude (const std::complex & x) //****************************************************************************** // // Purpose: // // C8_MAGNITUDE returns the magnitude of a C8. // // Discussion: // // A C8 is a double precision complex value. // // Modified: // // 22 October 2005 // // Author: // // John Burkardt // // Parameters: // // Input, complex X, the value whose norm is desired. // // Output, double C8_MAGNITUDE, the magnitude of X. // { return abs(x); //double magnitude; // //magnitude = sqrt(pow(x.real(),2) + pow(x.imag(),2)); // //return magnitude; } } // namespace quartic } // namespace Mt2