Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Public Member Functions | List of all members
Heed::Parabol Class Reference

#include <parabol.h>

Public Member Functions

double a () const
 
double b () const
 
double c () const
 
void put_a (const double fa)
 
void put_b (const double fb)
 
void put_c (const double fc)
 
 Parabol (void)
 Default constructor.
 
 Parabol (const Parabol &f)
 
 Parabol (double fa, double fb, double fc)
 
 Parabol (double x[3], double y[3])
 Constructor from three points.
 
 Parabol (double x[3], double y[3], int)
 
 Parabol (double x1, double x2, double x3, double y1, double y2, double y3)
 Constructor from 3 points.
 
double eval (const double x) const
 Evaluate the function.
 
int find_zero (double xzero[2]) const
 
double find_maxmin ()
 
double determinant () const
 

Detailed Description

Solution of quadratic equation. The main class is colled with a cut of the last character: Parabol to assure absence of coincidences with any other libraries and user files. This was inpired by a significant problem with link. The programs containig this class often don't want to be linked due to missing of references to functions from this class, if the file compiled from parabol.c is just included in library. The class parabol does not contain anything special, but this often happens. Why - unknown. The mentioned change of name did not solve this problem. The only solution is the inclusion of the object file parabol.o in the list of command supplied to linker. In general, this is very simple class, see definition.

Definition at line 30 of file parabol.h.

Constructor & Destructor Documentation

◆ Parabol() [1/6]

Heed::Parabol::Parabol ( void  )
inline

Default constructor.

Definition at line 52 of file parabol.h.

52: da(0.0), db(0.0), dc(0.0), s_det(0), s_dxzero(0) {}

◆ Parabol() [2/6]

Heed::Parabol::Parabol ( const Parabol &  f)

Definition at line 20 of file parabol.cpp.

20{ *this = f; }

◆ Parabol() [3/6]

Heed::Parabol::Parabol ( double  fa,
double  fb,
double  fc 
)
inline

Definition at line 54 of file parabol.h.

55 : da(fa), db(fb), dc(fc), s_det(0), s_dxzero(0) {}

◆ Parabol() [4/6]

Heed::Parabol::Parabol ( double  x[3],
double  y[3] 
)

Constructor from three points.

Definition at line 22 of file parabol.cpp.

22 : s_det(0), s_dxzero(0) {
23 mfunname("Parabol::Parabol(double x[3], double y[3])");
24
25 check_econd12a(x[0], ==, x[1], "x[2]=" << x[2] << " y[0]=" << y[0] << " y[1]="
26 << y[1] << " y[2]=" << y[2] << '\n',
27 mcerr);
28 check_econd12a(x[0], ==, x[2], "x[1]=" << x[1] << " y[0]=" << y[0] << " y[1]="
29 << y[1] << " y[2]=" << y[2] << '\n',
30 mcerr);
31 check_econd12a(x[1], ==, x[2], "x[0]=" << x[0] << " y[0]=" << y[0] << " y[1]="
32 << y[1] << " y[2]=" << y[2] << '\n',
33 mcerr);
34 DynArr<DoubleAc> mat(3, 3);
35 DynLinArr<DoubleAc> par(3);
36 DynLinArr<DoubleAc> f(3);
37 for (int i = 0; i < 3; ++i) {
38 f[i] = y[i];
39 mat.ac(i, 2) = 1.0;
40 mat.ac(i, 1) = x[i];
41 mat.ac(i, 0) = x[i] * x[i];
42 }
43 int ierr;
44 int szero;
45 DynArr<DoubleAc> mat_inv;
46 inverse_DynArr_prot(mat, mat_inv, szero, ierr);
47 // check_econd11a( ierr, != 0 , "should never happen\n", mcerr );
48 if (ierr == 0) {
49 par = mat_inv * f;
50 // Iprintdla_DoubleAc(mcout, par, 3);
51 // if(fabs(par[0]) == 0.0)
52 // da=0.0;
53 // else
54 // da=par[0];
55 da = par[0];
56 db = par[1];
57 dc = par[2];
58 } else {
59 da = 0.0;
60 DynLinArr<int> s_var(3);
61 s_var[0] = 0;
62 s_var[1] = 1;
63 s_var[2] = 1;
64 DynArr<DoubleAc> mat_inv1(3, 3);
65 // int ierr1;
66 // inverse_DynArr(mat, s_var, mat_inv, ierr, mat_inv1, ierr1);
67 inverse_DynArr_prot(mat, s_var, mat_inv, szero, ierr);
68 if (ierr != 0) {
69 // what if x1 and x2 are the same but the both differ from x0
70 // Then the following should help:
71 mat.ac(1, 1) = mat.ac(0, 1);
72 mat.ac(1, 2) = mat.ac(0, 2);
73 f[1] = f[0];
74 s_var[0] = 0;
75 s_var[1] = 1;
76 s_var[2] = 1;
77 inverse_DynArr_prot(mat, s_var, mat_inv, szero, ierr);
78 check_econd11a(ierr, != 0,
79 "should never happen\nmat=" << mat << "\ns_var=" << s_var
80 << "\nmat_inv=" << mat_inv,
81 mcerr);
82 }
83 par = mat_inv * f;
84 db = par[1];
85 dc = par[2];
86 }
87}
check_econd11a
#define check_econd11a(a, signb, add, stream)
Definition: FunNameStack.h:172
check_econd12a
#define check_econd12a(a, sign, b, add, stream)
Definition: FunNameStack.h:181
mfunname
#define mfunname(string)
Definition: FunNameStack.h:45
Heed::inverse_DynArr_prot
void inverse_DynArr_prot(const DynArr< DoubleAc > &mi, DynArr< DoubleAc > &mr, int &szero, int &serr, int s_stop)
Definition: inverse.cpp:17
mcerr
#define mcerr
Definition: prstream.h:128

◆ Parabol() [5/6]

Heed::Parabol::Parabol ( double  x[3],
double  y[3],
int   
)

Definition at line 165 of file parabol.cpp.

165 : s_det(0), s_dxzero(0) {
166 mfunname("Parabol::Parabol(double x[3], double y[3], int)");
167
168 check_econd12(x[0], ==, x[1], mcerr);
169 // check_econd12( x[0] , == , x[2] , mcerr);
170 // check_econd12( x[1] , == , x[2] , mcerr);
171
172 DynArr<DoubleAc> mat(3, 3);
173 DynLinArr<DoubleAc> par(3);
174 DynLinArr<DoubleAc> f(3);
175 for (int i = 0; i < 3; ++i) f[i] = y[i];
176 for (int i = 0; i < 2; ++i) {
177 mat.ac(i, 2) = 1.0;
178 mat.ac(i, 1) = x[i];
179 // Iprintdan(mcout, mat.ac(i,1));
180 mat.ac(i, 0) = x[i] * x[i];
181 }
182 mat.ac(2, 2) = 0.0;
183 mat.ac(2, 1) = 1.0;
184 mat.ac(2, 0) = 2.0 * x[2];
185 int ierr;
186 int szero;
187 DynArr<DoubleAc> mat_inv;
188 inverse_DynArr_prot(mat, mat_inv, szero, ierr);
189 check_econd11a(ierr, != 0, "should never happen\n", mcerr);
190 /*
191 if (ierr != 0) {
192 da = 0.0;
193 DynLinArr<int> s_var(3);
194 s_var[0] = 0;
195 s_var[1] = 1;
196 s_var[2] = 1;
197 inverse_DynArr(mat, mat_inv, ierr);
198 check_econd11a( ierr, != 0 , "should never happen\n", mcerr );
199 par = mat_inv * f;
200 db=par[1]; dc=par[2];
201 } else {
202 */
203 par = mat_inv * f;
204 if (fabs(par[0]) == 0.0) {
205 da = 0.0;
206 } else {
207 da = par[0];
208 db = par[1];
209 dc = par[2];
210 }
211 /*
212 HepMatrix mat(3,3,0);
213 HepVector par(3,0);
214 HepVector f(3,0);
215 for (int i = 0; i < 3; i++) f[i]=y[i];
216 for (int i = 0; i < 2; i++) {
217 mat[i][2] = 1.0;
218 mat[i][1] = x[i];
219 mat[i][0] = x[i] * x[i];
220 }
221 mat[2][2] = 0.0;
222 mat[2][1] = 1.0;
223 mat[2][0] = 2.0 * x[2];
224 par = solve(mat, f);
225 da = par[0]; db = par[1]; dc = par[2];
226 */
227}
check_econd12
#define check_econd12(a, sign, b, stream)
Definition: FunNameStack.h:163
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ Parabol() [6/6]

Heed::Parabol::Parabol ( double  x1,
double  x2,
double  x3,
double  y1,
double  y2,
double  y3 
)

Constructor from 3 points.

Definition at line 89 of file parabol.cpp.

91 : s_det(0), s_dxzero(0) {
92 mfunname("Parabol::Parabol(double x[3], double y[3])");
93
94 check_econd12a(x1, ==, x2, "x3=" << x3 << " y1=" << y1 << " y2=" << y2
95 << " y3=" << y3 << '\n',
96 mcerr);
97 check_econd12a(x1, ==, x3, "x2=" << x2 << " y1=" << y1 << " y2=" << y2
98 << " y3=" << y3 << '\n',
99 mcerr);
100 check_econd12a(x2, ==, x3, "x1=" << x1 << " y1=" << y1 << " y2=" << y2
101 << " y3=" << y3 << '\n',
102 mcerr);
103 DynArr<DoubleAc> mat(3, 3);
104 DynLinArr<DoubleAc> par(3);
105 DynLinArr<DoubleAc> f(3);
106 f[0] = y1;
107 mat.ac(0, 2) = 1.0;
108 mat.ac(0, 1) = x1;
109 mat.ac(0, 0) = x1 * x1;
110 f[1] = y2;
111 mat.ac(1, 2) = 1.0;
112 mat.ac(1, 1) = x2;
113 mat.ac(1, 0) = x2 * x2;
114 f[2] = y3;
115 mat.ac(2, 2) = 1.0;
116 mat.ac(2, 1) = x3;
117 mat.ac(2, 0) = x3 * x3;
118
119 int ierr;
120 int szero;
121 DynArr<DoubleAc> mat_inv;
122 inverse_DynArr_prot(mat, mat_inv, szero, ierr);
123 // check_econd11a( ierr, != 0 , "should never happen\n", mcerr );
124 if (ierr == 0) {
125 par = mat_inv * f;
126 // Iprintdla_DoubleAc(mcout, par, 3);
127 // if(fabs(par[0]) == 0.0)
128 // da=0.0;
129 // else
130 // da=par[0];
131 da = par[0];
132 db = par[1];
133 dc = par[2];
134 } else {
135 da = 0.0;
136 DynLinArr<int> s_var(3);
137 s_var[0] = 0;
138 s_var[1] = 1;
139 s_var[2] = 1;
140 DynArr<DoubleAc> mat_inv1(3, 3);
141 // int ierr1;
142 // inverse_DynArr(mat, s_var, mat_inv, ierr, mat_inv1, ierr1);
143 inverse_DynArr_prot(mat, s_var, mat_inv, szero, ierr);
144 if (ierr != 0) {
145 // what if x1 and x2 are the same but the both differ from x0
146 // Then the following should help:
147 mat.ac(1, 1) = mat.ac(0, 1);
148 mat.ac(1, 2) = mat.ac(0, 2);
149 f[1] = f[0];
150 s_var[0] = 0;
151 s_var[1] = 1;
152 s_var[2] = 1;
153 inverse_DynArr_prot(mat, s_var, mat_inv, szero, ierr);
154 check_econd11a(ierr, != 0,
155 "should never happen\nmat=" << mat << "\ns_var=" << s_var
156 << "\nmat_inv=" << mat_inv,
157 mcerr);
158 }
159 par = mat_inv * f;
160 db = par[1];
161 dc = par[2];
162 }
163}

Member Function Documentation

◆ a()

double Heed::Parabol::a ( ) const
inline

Definition at line 32 of file parabol.h.

32{ return da; }

Referenced by Heed::operator<<().

◆ b()

double Heed::Parabol::b ( ) const
inline

Definition at line 33 of file parabol.h.

33{ return db; }

Referenced by Heed::operator<<().

◆ c()

double Heed::Parabol::c ( ) const
inline

Definition at line 34 of file parabol.h.

34{ return dc; }

Referenced by Heed::operator<<().

◆ determinant()

double Heed::Parabol::determinant ( ) const
inline

Definition at line 72 of file parabol.h.

72 {
73 const Parabol& t = (*this);
74 if (s_det == 0) {
75 t.s_det = 1;
76 t.det = db * db - 4 * da * dc;
77 }
78 return det;
79 }
Heed::Parabol::Parabol
Parabol(void)
Default constructor.
Definition: parabol.h:52

Referenced by find_zero().

◆ eval()

double Heed::Parabol::eval ( const double  x) const
inline

Evaluate the function.

Definition at line 66 of file parabol.h.

66{ return da * x * x + db * x + dc; }

◆ find_maxmin()

double Heed::Parabol::find_maxmin ( void  )

Definition at line 274 of file parabol.cpp.

274 {
275 mfunname("double Parabol::find_maxmin(void)");
276 check_econd11(da, == 0, mcerr);
277 return -db / (2.0 * da);
278}
check_econd11
#define check_econd11(a, signb, stream)
Definition: FunNameStack.h:155

◆ find_zero()

int Heed::Parabol::find_zero ( double  xzero[2]) const

Definition at line 229 of file parabol.cpp.

229 {
230 mfunnamep("int Parabol::find_zero(double xzero[2]) const");
231 const Parabol& t = (*this);
232 if (s_dxzero == 0) {
233 // mcout<<"Parabol::find_zero: s_dxzero == 0\n";
234 t.s_dxzero = 1;
235 if (da == 0.0) {
236 if (db == 0.0) {
237 funnw.ehdr(mcerr);
238 mcerr << "can not find zero\n";
239 spexit(mcerr);
240 } else {
241 t.qdxzero = 1;
242 t.dxzero[0] = -dc / db;
243 }
244 } else {
245 if (determinant() < 0.0) {
246 t.qdxzero = 0;
247 t.dxzero[0] = 0;
248 t.dxzero[1] = 0;
249 } else if (determinant() == 0.0) {
250 t.qdxzero = 1;
251 t.dxzero[0] = -db / (2.0 * da);
252 } else {
253 const double sq = sqrt(determinant());
254 t.qdxzero = 2;
255 if (da > 0.0) {
256 t.dxzero[0] = (-db - sq) / (2.0 * da);
257 t.dxzero[1] = (-db + sq) / (2.0 * da);
258 } else {
259 t.dxzero[1] = (-db - sq) / (2.0 * da);
260 t.dxzero[0] = (-db + sq) / (2.0 * da);
261 }
262 // mcout<<"Parabol::find_zero: t.dxzero[0]="<<t.dxzero[0]
263 // <<" dxzero[0]="<<dxzero[0]<<'\n';
264 // mcout<<"Parabol::find_zero: t.dxzero[1]="<<t.dxzero[1]
265 // <<" dxzero[1]="<<dxzero[1]<<'\n';
266 }
267 }
268 }
269 xzero[0] = dxzero[0];
270 xzero[1] = dxzero[1];
271 return qdxzero;
272}
mfunnamep
#define mfunnamep(string)
Definition: FunNameStack.h:49
spexit
#define spexit(stream)
Definition: FunNameStack.h:256
Heed::Parabol::determinant
double determinant() const
Definition: parabol.h:72
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

Referenced by Heed::Cubic::find_maxmin(), and Heed::operator<<().

◆ put_a()

void Heed::Parabol::put_a ( const double  fa)
inline

Definition at line 35 of file parabol.h.

35 {
36 da = fa;
37 s_det = 0;
38 s_dxzero = 0;
39 }

◆ put_b()

void Heed::Parabol::put_b ( const double  fb)
inline

Definition at line 40 of file parabol.h.

40 {
41 db = fb;
42 s_det = 0;
43 s_dxzero = 0;
44 }

◆ put_c()

void Heed::Parabol::put_c ( const double  fc)
inline

Definition at line 45 of file parabol.h.

45 {
46 dc = fc;
47 s_det = 0;
48 s_dxzero = 0;
49 }
The documentation for this class was generated from the following files: