Mantid
Searching...
No Matches
mathSupport.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright &copy; 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7#include <algorithm>
8#include <cmath>
9#include <complex>
10#include <iterator>
11#include <vector>
12
14#include "MantidKernel/Logger.h"
15
16namespace Mantid {
17
18namespace {
19Kernel::Logger logger("mathSupport");
20}
21
22template <typename InputIter>
23int solveQuadratic(const InputIter Coef, std::pair<std::complex<double>, std::complex<double>> &OutAns)
31{
32 const double a = (*Coef);
33 const double b = *(Coef + 1);
34 const double c = *(Coef + 2);
35
36 if (a == 0.0) {
37 if (b == 0.0) {
38 OutAns.first = std::complex<double>(0.0, 0.0);
39 OutAns.second = std::complex<double>(0.0, 0.0);
40 return 0;
41 } else {
42 OutAns.first = std::complex<double>(-c / b, 0.0);
43 OutAns.second = OutAns.first;
44 return 1;
45 }
46 }
47 const double complex_part_sq = b * b - 4 * a * c;
48 if (complex_part_sq == 0.) { // degenerate case
49 OutAns.first = std::complex<double>(-0.5 * b / a, 0.0);
50 OutAns.second = OutAns.first;
51 return 1;
52 } else if (complex_part_sq > 0.) { /* Real Roots */
53 const double complex_part = sqrt(complex_part_sq);
54 const double q = (b >= 0) ? -0.5 * (b + complex_part) : -0.5 * (b - complex_part);
55 OutAns.first = std::complex<double>(q / a, 0.0);
56 OutAns.second = std::complex<double>(c / q, 0.0);
57 return 2;
58 } else {
59 const double complex_part = sqrt(-complex_part_sq);
60 std::complex<double> CQ(-0.5 * b, (b >= 0 ? -0.5 * complex_part : 0.5 * complex_part));
61 OutAns.first = CQ / a;
62 OutAns.second = c / CQ;
63 return 2;
64 }
65}
66
67template <typename CInputIter>
68int solveCubic(const CInputIter Coef, std::complex<double> &AnsA, std::complex<double> &AnsB,
69 std::complex<double> &AnsC)
78{
79 using Cpair = std::complex<double>;
80 double q, r; /* solution parameters */
81 double termR, discrim;
82 double r13;
83 std::pair<std::complex<double>, std::complex<double>> SQ;
84
85 if ((*Coef) == 0) {
86 const int xi = solveQuadratic(Coef + 1, SQ);
87 AnsA = SQ.first;
88 AnsB = SQ.second;
89 AnsC = SQ.second;
90 return xi;
91 }
92 if (*(Coef + 3) == 0) {
93 const int xi = solveQuadratic(Coef, SQ);
94 logger.debug() << "Xi == " << xi << '\n';
95 AnsA = SQ.first;
96 AnsB = (xi == 1) ? SQ.first : SQ.second;
97 AnsC = std::complex<double>(0.0, 0.0);
98 return (AnsC != AnsA) ? xi + 1 : xi;
99 }
100 const double a = (*Coef);
101 const double b = *(Coef + 1) / a;
102 const double c = *(Coef + 2) / a;
103 const double d = *(Coef + 3) / a;
104
105 q = (3.0 * c - (b * b)) / 9.0; // -q
106 r = -27.0 * d + b * (9.0 * c - 2.0 * b * b); // -r
107 r /= 54.0;
108
109 discrim = q * q * q + r * r; // r^2-qq^3
110 /* The first root is always real. */
111 termR = (b / 3.0);
112
113 if (discrim > 1e-13) /* one root real, two are complex */
114 {
115 double s = r + sqrt(discrim);
116 s = ((s < 0) ? -pow(-s, (1.0 / 3.0)) : pow(s, (1.0 / 3.0)));
117 double t = r - sqrt(discrim);
118 t = ((t < 0) ? -pow(-t, (1.0 / 3.0)) : pow(t, (1.0 / 3.0)));
119 AnsA = Cpair(-termR + s + t, 0.0);
120 // second real point.
121 termR += (s + t) / 2.0;
122 double termI = sqrt(3.0) * (-t + s) / 2;
123 AnsB = Cpair(-termR, termI);
124 AnsC = Cpair(-termR, -termI);
125 return 3;
126 }
127
128 /* The remaining options are all real */
129
130 if (discrim < 1e-13) // All roots real
131 {
132 q = -q;
133 double q3 = q * q * q;
134 q3 = acos(-r / sqrt(q3));
135 r13 = -2.0 * sqrt(q);
136 AnsA = Cpair(-termR + r13 * cos(q3 / 3.0), 0.0);
137 AnsB = Cpair(-termR + r13 * cos((q3 + 2.0 * M_PI) / 3.0), 0.0);
138 AnsC = Cpair(-termR + r13 * cos((q3 - 2.0 * M_PI) / 3.0), 0.0);
139 return 3;
140 }
141
142 // Only option left is that all roots are real and unequal
143 // (to get here, q*q*q=r*r)
144
145 r13 = ((r < 0) ? -pow(-r, (1.0 / 3.0)) : pow(r, (1.0 / 3.0)));
146 AnsA = Cpair(-termR + 2.0 * r13, 0.0);
147 AnsB = Cpair(-(r13 + termR), 0.0);
148 AnsC = Cpair(-(r13 + termR), 0.0);
149 return 2;
150}
151
153
154template MANTID_GEOMETRY_DLL int solveQuadratic(const double *,
155 std::pair<std::complex<double>, std::complex<double>> &);
156template MANTID_GEOMETRY_DLL int solveQuadratic(double *, std::pair<std::complex<double>, std::complex<double>> &);
158 std::pair<std::complex<double>, std::complex<double>> &);
159template MANTID_GEOMETRY_DLL int solveCubic(const double *, std::complex<double> &, std::complex<double> &,
160 std::complex<double> &);
161template MANTID_GEOMETRY_DLL int solveCubic(double *, std::complex<double> &, std::complex<double> &,
162 std::complex<double> &);
163
164template MANTID_GEOMETRY_DLL int solveCubic(const std::vector<double>::iterator, std::complex<double> &,
165 std::complex<double> &, std::complex<double> &);
166template MANTID_GEOMETRY_DLL int solveCubic(const std::vector<double>::const_iterator, std::complex<double> &,
167 std::complex<double> &, std::complex<double> &);
168
170} // namespace Mantid
Helper class which provides the Collimation Length for SANS instruments.
MANTID_GEOMETRY_DLL int solveCubic(InputIter, std::complex< double > &, std::complex< double > &, std::complex< double > &)
Solve a Cubic equation.
MANTID_GEOMETRY_DLL int solveQuadratic(InputIter, std::pair< std::complex< double >, std::complex< double > > &)