| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4QSS3
27 //
28 // G4QSS3 simulator
29
30 // Authors: Lucio Santi, Rodrigo Castro (Univ. Buenos Aires) - 2018-2021
31 // --------------------------------------------------------------------
32 #ifndef _G4QSS3_H_
33 #define _G4QSS3_H_ 1
34
35 #include "G4Types.hh"
36 #include "G4qss_misc.hh"
37
38 #include <cmath>
39
40 class G4QSS3
41 {
42 public:
43
44 G4QSS3(QSS_simulator);
45
46 inline QSS_simulator getSimulator() const { return this->simulator; }
47
48 inline G4int order() const { return 3; }
49
50 inline void full_definition(G4double coeff)
51 {
52 G4double* const x = simulator->q;
53 G4double* const dx = simulator->x;
54 G4double* const alg = simulator->alg;
55
56 dx[1] = x[12];
57 dx[2] = 0;
58 dx[3] = 0;
59
60 dx[5] = x[16];
61 dx[6] = 0;
62 dx[7] = 0;
63
64 dx[9] = x[20];
65 dx[10] = 0;
66 dx[11] = 0;
67
68 dx[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
69 dx[14] = 0;
70 dx[15] = 0;
71
72 dx[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
73 dx[18] = 0;
74 dx[19] = 0;
75
76 dx[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
77 dx[22] = 0;
78 dx[23] = 0;
79 }
80
81 inline void dependencies(G4int i, G4double coeff)
82 {
83 G4double* const x = simulator->q;
84 G4double* const der = simulator->x;
85 G4double* const alg = simulator->alg;
86
87 switch (i)
88 {
89 case 0:
90 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
91 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
92 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
93
94 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
95 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
96 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
97
98 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
99 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
100 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
101 return;
102 case 1:
103 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
104 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
105 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
106
107 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
108 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
109 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
110
111 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
112 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
113 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
114 return;
115 case 2:
116 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
117 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
118 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
119
120 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
121 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
122 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
123
124 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
125 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
126 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
127 return;
128 case 3:
129 der[1] = x[12];
130 der[2] = x[13] / 2;
131 der[3] = x[14] / 3;
132
133 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
134 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
135 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
136
137 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
138 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
139 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
140 return;
141 case 4:
142 der[5] = x[16];
143 der[6] = x[17] / 2;
144 der[7] = x[18] / 3;
145
146 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
147 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
148 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
149
150 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
151 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
152 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
153 return;
154 case 5:
155 der[9] = x[20];
156 der[10] = x[21] / 2;
157 der[11] = x[22] / 3;
158
159 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
160 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
161 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
162
163 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
164 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
165 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
166 return;
167 }
168 }
169
170 void recompute_next_times(G4int* inf, G4double t); // __attribute__((hot));
171
172 inline void recompute_all_state_times(G4double t)
173 {
174 G4double mpr;
175 G4double* const x = simulator->x;
176 G4double* const lqu = simulator->lqu;
177 G4double* const time = simulator->nextStateTime;
178
179 for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 4)
180 {
181 const G4int icf1 = icf0 + 1;
182
183 if (x[icf1] == 0)
184 {
185 time[var] = Qss_misc::INF;
186 }
187 else
188 {
189 mpr = lqu[var] / x[icf1];
190 if (mpr < 0) { mpr *= -1; }
191 time[var] = t + mpr;
192 }
193 }
194 }
195
196 inline void next_time(G4int i, G4double t)
197 {
198 const G4int cf3 = 4 * i + 3;
199 G4double* const x = simulator->x;
200 G4double* const lqu = simulator->lqu;
201 G4double* const time = simulator->nextStateTime;
202
203 if (likely(x[cf3])) {
204 time[i] = t + std::cbrt(lqu[i] / std::fabs(x[cf3]));
205 } else {
206 time[i] = Qss_misc::INF;
207 }
208 }
209
210 inline void update_quantized_state(G4int i)
211 {
212 const G4int cf0 = i * 4, cf1 = cf0 + 1, cf2 = cf1 + 1;
213 G4double* const q = simulator->q;
214 G4double* const x = simulator->x;
215
216 q[cf0] = x[cf0];
217 q[cf1] = x[cf1];
218 q[cf2] = x[cf2];
219 }
220
221 inline void reset_state(G4int i, G4double value)
222 {
223 G4double* const x = simulator->x;
224 G4double* const q = simulator->q;
225 G4double* const tq = simulator->tq;
226 G4double* const tx = simulator->tx;
227 const G4int idx = 4 * i;
228
229 x[idx] = value;
230
231 simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
232 if (simulator->lqu[i] < simulator->dQMin[i])
233 {
234 simulator->lqu[i] = simulator->dQMin[i];
235 }
236 q[idx] = value;
237 q[idx + 1] = q[idx + 2] = tq[i] = tx[i] = 0;
238 }
239
240 inline G4double evaluate_x_poly(G4int i, G4double dt, G4double* p)
241 {
242 return ((p[i + 3] * dt + p[i + 2]) * dt + p[i + 1]) * dt + p[i];
243 }
244
245 inline void advance_time_q(G4int i, G4double dt) // __attribute__((hot))
246 {
247 G4double* const p = simulator->q;
248 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
249 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
250 p[i1] = p[i1] + 2 * dt * p[i2];
251 }
252
253 inline void advance_time_x(G4int i, G4double dt) // __attribute__((hot))
254 {
255 G4double* const p = simulator->x;
256 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1, i3 = i2 + 1;
257 p[i0] = ((p[i3] * dt + p[i2]) * dt + p[i1]) * dt + p[i0];
258 p[i1] = (3 * p[i3] * dt + 2 * p[i2]) * dt + p[i1];
259 p[i2] = p[i2] + 3 * dt * p[i3];
260 }
261
262 G4double min_pos_root(G4double* coeff, G4int order);
263
264 inline G4double min_pos_root_2(G4double* coeff)
265 {
266 G4double mpr = Qss_misc::INF;
267
268 if (coeff[2] == 0 || (1000 * std::fabs(coeff[2])) < std::fabs(coeff[1]))
269 {
270 if (coeff[1] == 0) {
271 mpr = Qss_misc::INF;
272 } else {
273 mpr = -coeff[0] / coeff[1];
274 }
275
276 if (mpr < 0) { mpr = Qss_misc::INF; }
277 }
278 else
279 {
280 G4double disc;
281 disc = coeff[1] * coeff[1] - 4 * coeff[2] * coeff[0];
282 if (disc < 0) // no real roots
283 {
284 mpr = Qss_misc::INF;
285 }
286 else
287 {
288 G4double sd, r1;
289 G4double cf2_d2 = 2 * coeff[2];
290
291 sd = std::sqrt(disc);
292 r1 = (-coeff[1] + sd) / cf2_d2;
293 if (r1 > 0) {
294 mpr = r1;
295 } else {
296 mpr = Qss_misc::INF;
297 }
298 r1 = (-coeff[1] - sd) / cf2_d2;
299 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
300 }
301 }
302
303 return mpr;
304 } // __attribute__((hot))
305
306 inline G4double min_pos_root_3(G4double* coeff)
307 {
308 G4double mpr = Qss_misc::INF;
309 static const G4double sqrt3 = std::sqrt(3);
310
311 if ((coeff[3] == 0) || (1000 * std::fabs(coeff[3]) < std::fabs(coeff[2])))
312 {
313 mpr = min_pos_root_2(coeff);
314 }
315 else if (coeff[0] == 0)
316 {
317 if (coeff[1] == 0)
318 {
319 mpr = -coeff[2] / coeff[3];
320 }
321 else
322 {
323 coeff[0] = coeff[1];
324 coeff[1] = coeff[2];
325 coeff[2] = coeff[3];
326 mpr = min_pos_root_2(coeff);
327 }
328 }
329 else
330 {
331 G4double q, r, disc, q3;
332 G4double val = coeff[2] / 3 / coeff[3];
333 G4double cf32 = coeff[3] * coeff[3];
334 G4double cf22 = coeff[2] * coeff[2];
335 G4double denq = 9 * cf32;
336 G4double denr = 6 * coeff[3] * denq;
337 G4double rcomm = 9 * coeff[3] * coeff[2] * coeff[1] - 2 * cf22 * coeff[2];
338
339 q = (3 * coeff[3] * coeff[1] - cf22) / denq;
340 q3 = q * q * q;
341
342 r = (rcomm - 27 * cf32 * coeff[0]) / denr;
343 disc = q3 + r * r;
344 mpr = Qss_misc::INF;
345
346 if (disc >= 0)
347 {
348 G4double sd, sx, t, r1, rsd;
349 sd = std::sqrt(disc);
350 rsd = r + sd;
351 if (rsd > 0) {
352 sx = std::cbrt(rsd);
353 } else {
354 sx = -std::cbrt(std::fabs(rsd));
355 }
356
357 rsd = r - sd;
358 if (rsd > 0) {
359 t = std::cbrt(rsd);
360 } else {
361 t = -std::cbrt(std::fabs(rsd));
362 }
363
364 r1 = sx + t - val;
365
366 if (r1 > 0) { mpr = r1; }
367 }
368 else
369 {
370 // three real roots
371 G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1;
372 rho = std::sqrt(-q3);
373 th = std::acos(r / rho);
374 rho13 = std::cbrt(rho);
375 costh3 = std::cos(th / 3);
376 sinth3 = std::sin(th / 3);
377 spt = rho13 * 2 * costh3;
378 smti32 = -rho13 * sinth3 * sqrt3;
379 r1 = spt - val;
380 if (r1 > 0) { mpr = r1; }
381 r1 = -spt / 2 - val + smti32;
382 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
383 r1 = r1 - 2 * smti32;
384 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
385 }
386 }
387
388 return mpr;
389 } // __attribute__((hot))
390
391 inline G4double min_pos_root_2_alt(G4double* coeff, G4double cf0Alt)
392 {
393 G4double mpr = Qss_misc::INF;
394 G4double mpr2;
395
396 if (coeff[2] == 0 || (1000 * std::fabs(coeff[2])) < std::fabs(coeff[1]))
397 {
398 if (coeff[1] == 0)
399 {
400 mpr = Qss_misc::INF;
401 }
402 else
403 {
404 mpr = -coeff[0] / coeff[1];
405 mpr2 = -cf0Alt / coeff[1];
406 if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
407 }
408
409 if (mpr < 0) { mpr = Qss_misc::INF; }
410 }
411 else
412 {
413 G4double cf1_2 = coeff[1] * coeff[1];
414 G4double cf2_4 = 4 * coeff[2];
415 G4double disc1 = cf1_2 - cf2_4 * coeff[0];
416 G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
417 G4double cf2_d2 = 2 * coeff[2];
418
419 if (unlikely(disc1 < 0 && disc2 < 0))
420 {
421 mpr = Qss_misc::INF;
422 }
423 else if (disc2 < 0)
424 {
425 G4double sd, r1;
426 sd = std::sqrt(disc1);
427 r1 = (-coeff[1] + sd) / cf2_d2;
428 if (r1 > 0) {
429 mpr = r1;
430 } else {
431 mpr = Qss_misc::INF;
432 }
433 r1 = (-coeff[1] - sd) / cf2_d2;
434 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
435 }
436 else if (disc1 < 0)
437 {
438 G4double sd, r1;
439 sd = std::sqrt(disc2);
440 r1 = (-coeff[1] + sd) / cf2_d2;
441 if (r1 > 0) {
442 mpr = r1;
443 } else {
444 mpr = Qss_misc::INF;
445 }
446 r1 = (-coeff[1] - sd) / cf2_d2;
447 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
448 }
449 else
450 {
451 G4double sd1, r1, sd2, r2;
452 sd1 = std::sqrt(disc1);
453 sd2 = std::sqrt(disc2);
454 r1 = (-coeff[1] + sd1) / cf2_d2;
455 r2 = (-coeff[1] + sd2) / cf2_d2;
456
457 if (r1 > 0) {
458 mpr = r1;
459 } else {
460 mpr = Qss_misc::INF;
461 }
462 r1 = (-coeff[1] - sd1) / cf2_d2;
463 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
464
465 if (r2 > 0 && r2 < mpr) { mpr = r2; }
466 r2 = (-coeff[1] - sd2) / cf2_d2;
467 if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
468 }
469 }
470
471 return mpr;
472 } // __attribute__((hot))
473
474 inline G4double min_pos_root_3_alt(G4double* coeff, G4double cf0Alt)
475 {
476 G4double mpr = Qss_misc::INF;
477 static const G4double sqrt3 = std::sqrt(3);
478
479 if ((coeff[3] == 0) || (1000 * std::fabs(coeff[3]) < std::fabs(coeff[2])))
480 {
481 mpr = min_pos_root_2_alt(coeff, cf0Alt);
482 }
483 else if (coeff[0] == 0)
484 {
485 G4double mpr2;
486 coeff[0] = cf0Alt;
487 mpr = min_pos_root_3(coeff);
488
489 if (coeff[1] == 0)
490 {
491 mpr2 = -coeff[2] / coeff[3];
492 }
493 else
494 {
495 coeff[0] = coeff[1];
496 coeff[1] = coeff[2];
497 coeff[2] = coeff[3];
498 mpr2 = min_pos_root_2(coeff);
499 }
500
501 if (mpr2 > 0 && mpr2 < mpr) { mpr = mpr2; }
502 }
503 else if (cf0Alt == 0)
504 {
505 G4double mpr2;
506 mpr = min_pos_root_3(coeff);
507
508 if (coeff[1] == 0)
509 {
510 mpr2 = -coeff[2] / coeff[3];
511 }
512 else
513 {
514 coeff[0] = coeff[1];
515 coeff[1] = coeff[2];
516 coeff[2] = coeff[3];
517 mpr2 = min_pos_root_2(coeff);
518 }
519
520 if (mpr2 > 0 && mpr2 < mpr) { mpr = mpr2; }
521 }
522 else
523 {
524 G4double q, r, rAlt, disc, discAlt, q3;
525 G4double val = coeff[2] / 3 / coeff[3];
526 G4double cf32 = coeff[3] * coeff[3];
527 G4double cf22 = coeff[2] * coeff[2];
528 G4double denq = 9 * cf32;
529 G4double denr = 6 * coeff[3] * denq;
530 G4double rcomm = 9 * coeff[3] * coeff[2] * coeff[1] - 2 * cf22 * coeff[2];
531
532 q = (3 * coeff[3] * coeff[1] - cf22) / denq;
533 q3 = q * q * q;
534
535 r = (rcomm - 27 * cf32 * coeff[0]) / denr;
536 rAlt = (rcomm - 27 * cf32 * cf0Alt) / denr;
537
538 disc = q3 + r * r;
539 discAlt = q3 + rAlt * rAlt;
540 mpr = Qss_misc::INF;
541
542 if (disc >= 0)
543 {
544 G4double sd, sx, t, r1, rsd;
545 sd = std::sqrt(disc);
546 rsd = r + sd;
547 if (rsd > 0) {
548 sx = std::cbrt(rsd);
549 } else {
550 sx = -std::cbrt(std::fabs(rsd));
551 }
552
553 rsd = r - sd;
554 if (rsd > 0) {
555 t = std::cbrt(rsd);
556 } else {
557 t = -std::cbrt(std::fabs(rsd));
558 }
559
560 r1 = sx + t - val;
561
562 if (r1 > 0) { mpr = r1; }
563
564 if (discAlt >= 0)
565 {
566 G4double sdAlt, sAlt, tAlt, r1Alt, rsdAlt;
567 sdAlt = std::sqrt(discAlt);
568 rsdAlt = rAlt + sdAlt;
569 if (rsdAlt > 0) {
570 sAlt = std::cbrt(rsdAlt);
571 } else {
572 sAlt = -std::cbrt(std::fabs(rsdAlt));
573 }
574
575 rsdAlt = rAlt - sdAlt;
576 if (rsdAlt > 0) {
577 tAlt = std::cbrt(rsdAlt);
578 } else {
579 tAlt = -std::cbrt(std::fabs(rsdAlt));
580 }
581
582 r1Alt = sAlt + tAlt - val;
583
584 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
585 }
586 else
587 {
588 G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1Alt;
589
590 rho = std::sqrt(-q3);
591 th = std::acos(rAlt / rho);
592 rho13 = std::cbrt(rho);
593 costh3 = std::cos(th / 3);
594 sinth3 = std::sin(th / 3);
595 spt = rho13 * 2 * costh3;
596 smti32 = -rho13 * sinth3 * sqrt3;
597 r1Alt = spt - val;
598 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
599 r1Alt = -spt / 2 - val + smti32;
600 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
601 r1Alt = r1Alt - 2 * smti32;
602 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
603 }
604 }
605 else
606 {
607 G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1;
608
609 rho = std::sqrt(-q3);
610 th = std::acos(r / rho);
611 rho13 = std::cbrt(rho);
612 costh3 = std::cos(th / 3);
613 sinth3 = std::sin(th / 3);
614 spt = rho13 * 2 * costh3;
615 smti32 = -rho13 * sinth3 * sqrt3;
616 r1 = spt - val;
617 if (r1 > 0) { mpr = r1; }
618 r1 = -spt / 2 - val + smti32;
619 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
620 r1 = r1 - 2 * smti32;
621 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
622
623 if (discAlt >= 0)
624 {
625 G4double sdAlt, sAlt, tAlt, r1Alt, rsdAlt;
626 sdAlt = std::sqrt(discAlt);
627 rsdAlt = rAlt + sdAlt;
628 if (rsdAlt > 0) {
629 sAlt = std::cbrt(rsdAlt);
630 } else {
631 sAlt = -std::cbrt(std::fabs(rsdAlt));
632 }
633
634 rsdAlt = rAlt - sdAlt;
635 if (rsdAlt > 0) {
636 tAlt = std::cbrt(rsdAlt);
637 } else {
638 tAlt = -std::cbrt(std::fabs(rsdAlt));
639 }
640
641 r1Alt = sAlt + tAlt - val;
642
643 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
644 }
645 else
646 {
647 G4double thAlt, costh3Alt, sinth3Alt, sptAlt, smti32Alt, r1Alt;
648 thAlt = std::acos(rAlt / rho);
649 costh3Alt = std::cos(thAlt / 3);
650 sinth3Alt = std::sin(thAlt / 3);
651 sptAlt = rho13 * 2 * costh3Alt;
652 smti32Alt = -rho13 * sinth3Alt * sqrt3;
653 r1Alt = sptAlt - val;
654 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
655 r1Alt = -sptAlt / 2 - val + smti32Alt;
656 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
657 r1Alt = r1Alt - 2 * smti32Alt;
658 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
659 }
660 }
661 }
662
663 return mpr;
664 }
665
666 private:
667
668 QSS_simulator simulator;
669 };
670
671 #endif
672