| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4AnalyticalPolSolver class implementation 26 // G4AnalyticalPolSolver class implementation
27 // 27 //
28 // Author: V.Grichine, 24.04.97 28 // Author: V.Grichine, 24.04.97
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "globals.hh" 31 #include "globals.hh"
32 #include <complex> 32 #include <complex>
33 33
34 #include "G4AnalyticalPolSolver.hh" 34 #include "G4AnalyticalPolSolver.hh"
35 35
36 ////////////////////////////////////////////// 36 //////////////////////////////////////////////////////////////////////////////
37 37
38 G4AnalyticalPolSolver::G4AnalyticalPolSolver() 38 G4AnalyticalPolSolver::G4AnalyticalPolSolver() { ; }
39 39
40 ////////////////////////////////////////////// 40 //////////////////////////////////////////////////////////////////////////////
41 41
42 G4AnalyticalPolSolver::~G4AnalyticalPolSolver( 42 G4AnalyticalPolSolver::~G4AnalyticalPolSolver() { ; }
43 43
44 ////////////////////////////////////////////// 44 //////////////////////////////////////////////////////////////////////////////
45 // 45 //
46 // Array r[3][5] p[5] 46 // Array r[3][5] p[5]
47 // Roots of poly p[0] x^2 + p[1] x+p[2]=0 47 // Roots of poly p[0] x^2 + p[1] x+p[2]=0
48 // 48 //
49 // x = r[1][k] + i r[2][k]; k = 1, 2 49 // x = r[1][k] + i r[2][k]; k = 1, 2
50 50
51 G4int G4AnalyticalPolSolver::QuadRoots(G4doubl 51 G4int G4AnalyticalPolSolver::QuadRoots(G4double p[5], G4double r[3][5])
52 { 52 {
53 G4double b, c, d2, d; 53 G4double b, c, d2, d;
54 54
55 b = -p[1] / p[0] / 2.; 55 b = -p[1] / p[0] / 2.;
56 c = p[2] / p[0]; 56 c = p[2] / p[0];
57 d2 = b * b - c; 57 d2 = b * b - c;
58 58
59 if(d2 >= 0.) 59 if(d2 >= 0.)
60 { 60 {
61 d = std::sqrt(d2); 61 d = std::sqrt(d2);
62 r[1][1] = b - d; 62 r[1][1] = b - d;
63 r[1][2] = b + d; 63 r[1][2] = b + d;
64 r[2][1] = 0.; 64 r[2][1] = 0.;
65 r[2][2] = 0.; 65 r[2][2] = 0.;
66 } 66 }
67 else 67 else
68 { 68 {
69 d = std::sqrt(-d2); 69 d = std::sqrt(-d2);
70 r[2][1] = d; 70 r[2][1] = d;
71 r[2][2] = -d; 71 r[2][2] = -d;
72 r[1][1] = b; 72 r[1][1] = b;
73 r[1][2] = b; 73 r[1][2] = b;
74 } 74 }
75 75
76 return 2; 76 return 2;
77 } 77 }
78 78
79 ////////////////////////////////////////////// 79 //////////////////////////////////////////////////////////////////////////////
80 // 80 //
81 // Array r[3][5] p[5] 81 // Array r[3][5] p[5]
82 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0 82 // Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
83 // x=r[1][k] + i r[2][k] k=1,...,3 83 // x=r[1][k] + i r[2][k] k=1,...,3
84 // Assumes 0<arctan(x)<pi/2 for x>0 84 // Assumes 0<arctan(x)<pi/2 for x>0
85 85
86 G4int G4AnalyticalPolSolver::CubicRoots(G4doub 86 G4int G4AnalyticalPolSolver::CubicRoots(G4double p[5], G4double r[3][5])
87 { 87 {
88 G4double x, t, b, c, d; 88 G4double x, t, b, c, d;
89 G4int k; 89 G4int k;
90 90
91 if(p[0] != 1.) 91 if(p[0] != 1.)
92 { 92 {
93 for(k = 1; k < 4; k++) 93 for(k = 1; k < 4; k++)
94 { 94 {
95 p[k] = p[k] / p[0]; 95 p[k] = p[k] / p[0];
96 } 96 }
97 p[0] = 1.; 97 p[0] = 1.;
98 } 98 }
99 x = p[1] / 3.0; 99 x = p[1] / 3.0;
100 t = x * p[1]; 100 t = x * p[1];
101 b = 0.5 * (x * (t / 1.5 - p[2]) + p[3]); 101 b = 0.5 * (x * (t / 1.5 - p[2]) + p[3]);
102 t = (t - p[2]) / 3.0; 102 t = (t - p[2]) / 3.0;
103 c = t * t * t; 103 c = t * t * t;
104 d = b * b - c; 104 d = b * b - c;
105 105
106 if(d >= 0.) 106 if(d >= 0.)
107 { 107 {
108 d = std::pow((std::sqrt(d) + std::fabs(b)) 108 d = std::pow((std::sqrt(d) + std::fabs(b)), 1.0 / 3.0);
109 109
110 if(d != 0.) 110 if(d != 0.)
111 { 111 {
112 if(b > 0.) 112 if(b > 0.)
113 { 113 {
114 b = -d; 114 b = -d;
115 } 115 }
116 else 116 else
117 { 117 {
118 b = d; 118 b = d;
119 } 119 }
120 c = t / b; 120 c = t / b;
121 } 121 }
122 d = std::sqrt(0.75) * (b - c); 122 d = std::sqrt(0.75) * (b - c);
123 r[2][2] = d; 123 r[2][2] = d;
124 b = b + c; 124 b = b + c;
125 c = -0.5 * b - x; 125 c = -0.5 * b - x;
126 r[1][2] = c; 126 r[1][2] = c;
127 127
128 if((b > 0. && x <= 0.) || (b < 0. && x > 0 128 if((b > 0. && x <= 0.) || (b < 0. && x > 0.))
129 { 129 {
130 r[1][1] = c; 130 r[1][1] = c;
131 r[2][1] = -d; 131 r[2][1] = -d;
132 r[1][3] = b - x; 132 r[1][3] = b - x;
133 r[2][3] = 0; 133 r[2][3] = 0;
134 } 134 }
135 else 135 else
136 { 136 {
137 r[1][1] = b - x; 137 r[1][1] = b - x;
138 r[2][1] = 0.; 138 r[2][1] = 0.;
139 r[1][3] = c; 139 r[1][3] = c;
140 r[2][3] = -d; 140 r[2][3] = -d;
141 } 141 }
142 } // end of 2 equal or complex roots 142 } // end of 2 equal or complex roots
143 else 143 else
144 { 144 {
145 if(b == 0.) 145 if(b == 0.)
146 { 146 {
147 d = std::atan(1.0) / 1.5; 147 d = std::atan(1.0) / 1.5;
148 } 148 }
149 else 149 else
150 { 150 {
151 d = std::atan(std::sqrt(-d) / std::fabs( 151 d = std::atan(std::sqrt(-d) / std::fabs(b)) / 3.0;
152 } 152 }
153 153
154 if(b < 0.) 154 if(b < 0.)
155 { 155 {
156 b = std::sqrt(t) * 2.0; 156 b = std::sqrt(t) * 2.0;
157 } 157 }
158 else 158 else
159 { 159 {
160 b = -2.0 * std::sqrt(t); 160 b = -2.0 * std::sqrt(t);
161 } 161 }
162 162
163 c = std::cos(d) * b; 163 c = std::cos(d) * b;
164 t = -std::sqrt(0.75) * std::sin(d) * b - 0 164 t = -std::sqrt(0.75) * std::sin(d) * b - 0.5 * c;
165 d = -t - c - x; 165 d = -t - c - x;
166 c = c - x; 166 c = c - x;
167 t = t - x; 167 t = t - x;
168 168
169 if(std::fabs(c) > std::fabs(t)) 169 if(std::fabs(c) > std::fabs(t))
170 { 170 {
171 r[1][3] = c; 171 r[1][3] = c;
172 } 172 }
173 else 173 else
174 { 174 {
175 r[1][3] = t; 175 r[1][3] = t;
176 t = c; 176 t = c;
177 } 177 }
178 if(std::fabs(d) > std::fabs(t)) 178 if(std::fabs(d) > std::fabs(t))
179 { 179 {
180 r[1][2] = d; 180 r[1][2] = d;
181 } 181 }
182 else 182 else
183 { 183 {
184 r[1][2] = t; 184 r[1][2] = t;
185 t = d; 185 t = d;
186 } 186 }
187 r[1][1] = t; 187 r[1][1] = t;
188 188
189 for(k = 1; k < 4; k++) 189 for(k = 1; k < 4; k++)
190 { 190 {
191 r[2][k] = 0.; 191 r[2][k] = 0.;
192 } 192 }
193 } 193 }
194 return 0; 194 return 0;
195 } 195 }
196 196
197 ////////////////////////////////////////////// 197 //////////////////////////////////////////////////////////////////////////////
198 // 198 //
199 // Array r[3][5] p[5] 199 // Array r[3][5] p[5]
200 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0 200 // Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
201 // x=r[1][k] + i r[2][k] k=1,...,4 201 // x=r[1][k] + i r[2][k] k=1,...,4
202 202
203 G4int G4AnalyticalPolSolver::BiquadRoots(G4dou 203 G4int G4AnalyticalPolSolver::BiquadRoots(G4double p[5], G4double r[3][5])
204 { 204 {
205 G4double a, b, c, d, e; 205 G4double a, b, c, d, e;
206 G4int i, k, j; 206 G4int i, k, j;
207 207
208 if(p[0] != 1.0) 208 if(p[0] != 1.0)
209 { 209 {
210 for(k = 1; k < 5; k++) 210 for(k = 1; k < 5; k++)
211 { 211 {
212 p[k] = p[k] / p[0]; 212 p[k] = p[k] / p[0];
213 } 213 }
214 p[0] = 1.; 214 p[0] = 1.;
215 } 215 }
216 e = 0.25 * p[1]; 216 e = 0.25 * p[1];
217 b = 2 * e; 217 b = 2 * e;
218 c = b * b; 218 c = b * b;
219 d = 0.75 * c; 219 d = 0.75 * c;
220 b = p[3] + b * (c - p[2]); 220 b = p[3] + b * (c - p[2]);
221 a = p[2] - d; 221 a = p[2] - d;
222 c = p[4] + e * (e * a - p[3]); 222 c = p[4] + e * (e * a - p[3]);
223 a = a - d; 223 a = a - d;
224 224
225 p[1] = 0.5 * a; 225 p[1] = 0.5 * a;
226 p[2] = (p[1] * p[1] - c) * 0.25; 226 p[2] = (p[1] * p[1] - c) * 0.25;
227 p[3] = b * b / (-64.0); 227 p[3] = b * b / (-64.0);
228 228
229 if(p[3] < 0.) 229 if(p[3] < 0.)
230 { 230 {
231 CubicRoots(p, r); 231 CubicRoots(p, r);
232 232
233 for(k = 1; k < 4; k++) 233 for(k = 1; k < 4; k++)
234 { 234 {
235 if(r[2][k] == 0. && r[1][k] > 0) 235 if(r[2][k] == 0. && r[1][k] > 0)
236 { 236 {
237 d = r[1][k] * 4; 237 d = r[1][k] * 4;
238 a = a + d; 238 a = a + d;
239 239
240 if(a >= 0. && b >= 0.) 240 if(a >= 0. && b >= 0.)
241 { 241 {
242 p[1] = std::sqrt(d); 242 p[1] = std::sqrt(d);
243 } 243 }
244 else if(a <= 0. && b <= 0.) 244 else if(a <= 0. && b <= 0.)
245 { 245 {
246 p[1] = std::sqrt(d); 246 p[1] = std::sqrt(d);
247 } 247 }
248 else 248 else
249 { 249 {
250 p[1] = -std::sqrt(d); 250 p[1] = -std::sqrt(d);
251 } 251 }
252 252
253 b = 0.5 * (a + b / p[1]); 253 b = 0.5 * (a + b / p[1]);
254 254
255 p[2] = c / b; 255 p[2] = c / b;
256 QuadRoots(p, r); 256 QuadRoots(p, r);
257 257
258 for(i = 1; i < 3; i++) 258 for(i = 1; i < 3; i++)
259 { 259 {
260 for(j = 1; j < 3; j++) 260 for(j = 1; j < 3; j++)
261 { 261 {
262 r[j][i + 2] = r[j][i]; 262 r[j][i + 2] = r[j][i];
263 } 263 }
264 } 264 }
265 p[1] = -p[1]; 265 p[1] = -p[1];
266 p[2] = b; 266 p[2] = b;
267 QuadRoots(p, r); 267 QuadRoots(p, r);
268 268
269 for(i = 1; i < 5; i++) 269 for(i = 1; i < 5; i++)
270 { 270 {
271 r[1][i] = r[1][i] - e; 271 r[1][i] = r[1][i] - e;
272 } 272 }
273 273
274 return 4; 274 return 4;
275 } 275 }
276 } 276 }
277 } 277 }
278 if(p[2] < 0.) 278 if(p[2] < 0.)
279 { 279 {
280 b = std::sqrt(c); 280 b = std::sqrt(c);
281 d = b + b - a; 281 d = b + b - a;
282 p[1] = 0.; 282 p[1] = 0.;
283 283
284 if(d > 0.) 284 if(d > 0.)
285 { 285 {
286 p[1] = std::sqrt(d); 286 p[1] = std::sqrt(d);
287 } 287 }
288 } 288 }
289 else 289 else
290 { 290 {
291 if(p[1] > 0.) 291 if(p[1] > 0.)
292 { 292 {
293 b = std::sqrt(p[2]) * 2.0 + p[1]; 293 b = std::sqrt(p[2]) * 2.0 + p[1];
294 } 294 }
295 else 295 else
296 { 296 {
297 b = -std::sqrt(p[2]) * 2.0 + p[1]; 297 b = -std::sqrt(p[2]) * 2.0 + p[1];
298 } 298 }
299 299
300 if(b != 0.) 300 if(b != 0.)
301 { 301 {
302 p[1] = 0; 302 p[1] = 0;
303 } 303 }
304 else 304 else
305 { 305 {
306 for(k = 1; k < 5; k++) 306 for(k = 1; k < 5; k++)
307 { 307 {
308 r[1][k] = -e; 308 r[1][k] = -e;
309 r[2][k] = 0; 309 r[2][k] = 0;
310 } 310 }
311 return 0; 311 return 0;
312 } 312 }
313 } 313 }
314 314
315 p[2] = c / b; 315 p[2] = c / b;
316 QuadRoots(p, r); 316 QuadRoots(p, r);
317 317
318 for(k = 1; k < 3; k++) 318 for(k = 1; k < 3; k++)
319 { 319 {
320 for(j = 1; j < 3; j++) 320 for(j = 1; j < 3; j++)
321 { 321 {
322 r[j][k + 2] = r[j][k]; 322 r[j][k + 2] = r[j][k];
323 } 323 }
324 } 324 }
325 p[1] = -p[1]; 325 p[1] = -p[1];
326 p[2] = b; 326 p[2] = b;
327 QuadRoots(p, r); 327 QuadRoots(p, r);
328 328
329 for(k = 1; k < 5; k++) 329 for(k = 1; k < 5; k++)
330 { 330 {
331 r[1][k] = r[1][k] - e; 331 r[1][k] = r[1][k] - e;
332 } 332 }
333 333
334 return 4; 334 return 4;
335 } 335 }
336 336
337 ////////////////////////////////////////////// 337 //////////////////////////////////////////////////////////////////////////////
338 338
339 G4int G4AnalyticalPolSolver::QuarticRoots(G4do 339 G4int G4AnalyticalPolSolver::QuarticRoots(G4double p[5], G4double r[3][5])
340 { 340 {
341 G4double a0, a1, a2, a3, y1; 341 G4double a0, a1, a2, a3, y1;
342 G4double R2, D2, E2, D, E, R = 0.; 342 G4double R2, D2, E2, D, E, R = 0.;
343 G4double a, b, c, d, ds; 343 G4double a, b, c, d, ds;
344 344
345 G4double reRoot[4]; 345 G4double reRoot[4];
346 G4int k; 346 G4int k;
347 347
348 for(k = 0; k < 4; k++) 348 for(k = 0; k < 4; k++)
349 { 349 {
350 reRoot[k] = DBL_MAX; 350 reRoot[k] = DBL_MAX;
351 } 351 }
352 352
353 if(p[0] != 1.0) 353 if(p[0] != 1.0)
354 { 354 {
355 for(k = 1; k < 5; k++) 355 for(k = 1; k < 5; k++)
356 { 356 {
357 p[k] = p[k] / p[0]; 357 p[k] = p[k] / p[0];
358 } 358 }
359 p[0] = 1.; 359 p[0] = 1.;
360 } 360 }
361 a3 = p[1]; 361 a3 = p[1];
362 a2 = p[2]; 362 a2 = p[2];
363 a1 = p[3]; 363 a1 = p[3];
364 a0 = p[4]; 364 a0 = p[4];
365 365
366 // resolvent cubic equation cofs: 366 // resolvent cubic equation cofs:
367 367
368 p[1] = -a2; 368 p[1] = -a2;
369 p[2] = a1 * a3 - 4 * a0; 369 p[2] = a1 * a3 - 4 * a0;
370 p[3] = 4 * a2 * a0 - a1 * a1 - a3 * a3 * a0; 370 p[3] = 4 * a2 * a0 - a1 * a1 - a3 * a3 * a0;
371 371
372 CubicRoots(p, r); 372 CubicRoots(p, r);
373 373
374 for(k = 1; k < 4; k++) 374 for(k = 1; k < 4; k++)
375 { 375 {
376 if(r[2][k] == 0.) // find a real root 376 if(r[2][k] == 0.) // find a real root
377 { 377 {
378 reRoot[k] = r[1][k]; 378 reRoot[k] = r[1][k];
379 } 379 }
380 else 380 else
381 { 381 {
382 reRoot[k] = DBL_MAX; // kInfinity; 382 reRoot[k] = DBL_MAX; // kInfinity;
383 } 383 }
384 } 384 }
385 y1 = DBL_MAX; // kInfinity; 385 y1 = DBL_MAX; // kInfinity;
386 for(k = 1; k < 4; k++) 386 for(k = 1; k < 4; k++)
387 { 387 {
388 if(reRoot[k] < y1) 388 if(reRoot[k] < y1)
389 { 389 {
390 y1 = reRoot[k]; 390 y1 = reRoot[k];
391 } 391 }
392 } 392 }
393 393
394 R2 = 0.25 * a3 * a3 - a2 + y1; 394 R2 = 0.25 * a3 * a3 - a2 + y1;
395 b = 0.25 * (4 * a3 * a2 - 8 * a1 - a3 * a3 395 b = 0.25 * (4 * a3 * a2 - 8 * a1 - a3 * a3 * a3);
396 c = 0.75 * a3 * a3 - 2 * a2; 396 c = 0.75 * a3 * a3 - 2 * a2;
397 a = c - R2; 397 a = c - R2;
398 d = 4 * y1 * y1 - 16 * a0; 398 d = 4 * y1 * y1 - 16 * a0;
399 399
400 if(R2 > 0.) 400 if(R2 > 0.)
401 { 401 {
402 R = std::sqrt(R2); 402 R = std::sqrt(R2);
403 D2 = a + b / R; 403 D2 = a + b / R;
404 E2 = a - b / R; 404 E2 = a - b / R;
405 405
406 if(D2 >= 0.) 406 if(D2 >= 0.)
407 { 407 {
408 D = std::sqrt(D2); 408 D = std::sqrt(D2);
409 r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D 409 r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D;
410 r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D 410 r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D;
411 r[2][1] = 0.; 411 r[2][1] = 0.;
412 r[2][2] = 0.; 412 r[2][2] = 0.;
413 } 413 }
414 else 414 else
415 { 415 {
416 D = std::sqrt(-D2); 416 D = std::sqrt(-D2);
417 r[1][1] = -0.25 * a3 + 0.5 * R; 417 r[1][1] = -0.25 * a3 + 0.5 * R;
418 r[1][2] = -0.25 * a3 + 0.5 * R; 418 r[1][2] = -0.25 * a3 + 0.5 * R;
419 r[2][1] = 0.5 * D; 419 r[2][1] = 0.5 * D;
420 r[2][2] = -0.5 * D; 420 r[2][2] = -0.5 * D;
421 } 421 }
422 if(E2 >= 0.) 422 if(E2 >= 0.)
423 { 423 {
424 E = std::sqrt(E2); 424 E = std::sqrt(E2);
425 r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E 425 r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E;
426 r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E 426 r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E;
427 r[2][3] = 0.; 427 r[2][3] = 0.;
428 r[2][4] = 0.; 428 r[2][4] = 0.;
429 } 429 }
430 else 430 else
431 { 431 {
432 E = std::sqrt(-E2); 432 E = std::sqrt(-E2);
433 r[1][3] = -0.25 * a3 - 0.5 * R; 433 r[1][3] = -0.25 * a3 - 0.5 * R;
434 r[1][4] = -0.25 * a3 - 0.5 * R; 434 r[1][4] = -0.25 * a3 - 0.5 * R;
435 r[2][3] = 0.5 * E; 435 r[2][3] = 0.5 * E;
436 r[2][4] = -0.5 * E; 436 r[2][4] = -0.5 * E;
437 } 437 }
438 } 438 }
439 else if(R2 < 0.) 439 else if(R2 < 0.)
440 { 440 {
441 R = std::sqrt(-R2); 441 R = std::sqrt(-R2);
442 G4complex CD2(a, -b / R); 442 G4complex CD2(a, -b / R);
443 G4complex CD = std::sqrt(CD2); 443 G4complex CD = std::sqrt(CD2);
444 444
445 r[1][1] = -0.25 * a3 + 0.5 * real(CD); 445 r[1][1] = -0.25 * a3 + 0.5 * real(CD);
446 r[1][2] = -0.25 * a3 - 0.5 * real(CD); 446 r[1][2] = -0.25 * a3 - 0.5 * real(CD);
447 r[2][1] = 0.5 * R + 0.5 * imag(CD); 447 r[2][1] = 0.5 * R + 0.5 * imag(CD);
448 r[2][2] = 0.5 * R - 0.5 * imag(CD); 448 r[2][2] = 0.5 * R - 0.5 * imag(CD);
449 G4complex CE2(a, b / R); 449 G4complex CE2(a, b / R);
450 G4complex CE = std::sqrt(CE2); 450 G4complex CE = std::sqrt(CE2);
451 451
452 r[1][3] = -0.25 * a3 + 0.5 * real(CE); 452 r[1][3] = -0.25 * a3 + 0.5 * real(CE);
453 r[1][4] = -0.25 * a3 - 0.5 * real(CE); 453 r[1][4] = -0.25 * a3 - 0.5 * real(CE);
454 r[2][3] = -0.5 * R + 0.5 * imag(CE); 454 r[2][3] = -0.5 * R + 0.5 * imag(CE);
455 r[2][4] = -0.5 * R - 0.5 * imag(CE); 455 r[2][4] = -0.5 * R - 0.5 * imag(CE);
456 } 456 }
457 else // R2=0 case 457 else // R2=0 case
458 { 458 {
459 if(d >= 0.) 459 if(d >= 0.)
460 { 460 {
461 D2 = c + std::sqrt(d); 461 D2 = c + std::sqrt(d);
462 E2 = c - std::sqrt(d); 462 E2 = c - std::sqrt(d);
463 463
464 if(D2 >= 0.) 464 if(D2 >= 0.)
465 { 465 {
466 D = std::sqrt(D2); 466 D = std::sqrt(D2);
467 r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * 467 r[1][1] = -0.25 * a3 + 0.5 * R + 0.5 * D;
468 r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * 468 r[1][2] = -0.25 * a3 + 0.5 * R - 0.5 * D;
469 r[2][1] = 0.; 469 r[2][1] = 0.;
470 r[2][2] = 0.; 470 r[2][2] = 0.;
471 } 471 }
472 else 472 else
473 { 473 {
474 D = std::sqrt(-D2); 474 D = std::sqrt(-D2);
475 r[1][1] = -0.25 * a3 + 0.5 * R; 475 r[1][1] = -0.25 * a3 + 0.5 * R;
476 r[1][2] = -0.25 * a3 + 0.5 * R; 476 r[1][2] = -0.25 * a3 + 0.5 * R;
477 r[2][1] = 0.5 * D; 477 r[2][1] = 0.5 * D;
478 r[2][2] = -0.5 * D; 478 r[2][2] = -0.5 * D;
479 } 479 }
480 if(E2 >= 0.) 480 if(E2 >= 0.)
481 { 481 {
482 E = std::sqrt(E2); 482 E = std::sqrt(E2);
483 r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * 483 r[1][3] = -0.25 * a3 - 0.5 * R + 0.5 * E;
484 r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * 484 r[1][4] = -0.25 * a3 - 0.5 * R - 0.5 * E;
485 r[2][3] = 0.; 485 r[2][3] = 0.;
486 r[2][4] = 0.; 486 r[2][4] = 0.;
487 } 487 }
488 else 488 else
489 { 489 {
490 E = std::sqrt(-E2); 490 E = std::sqrt(-E2);
491 r[1][3] = -0.25 * a3 - 0.5 * R; 491 r[1][3] = -0.25 * a3 - 0.5 * R;
492 r[1][4] = -0.25 * a3 - 0.5 * R; 492 r[1][4] = -0.25 * a3 - 0.5 * R;
493 r[2][3] = 0.5 * E; 493 r[2][3] = 0.5 * E;
494 r[2][4] = -0.5 * E; 494 r[2][4] = -0.5 * E;
495 } 495 }
496 } 496 }
497 else 497 else
498 { 498 {
499 ds = std::sqrt(-d); 499 ds = std::sqrt(-d);
500 G4complex CD2(c, ds); 500 G4complex CD2(c, ds);
501 G4complex CD = std::sqrt(CD2); 501 G4complex CD = std::sqrt(CD2);
502 502
503 r[1][1] = -0.25 * a3 + 0.5 * real(CD); 503 r[1][1] = -0.25 * a3 + 0.5 * real(CD);
504 r[1][2] = -0.25 * a3 - 0.5 * real(CD); 504 r[1][2] = -0.25 * a3 - 0.5 * real(CD);
505 r[2][1] = 0.5 * R + 0.5 * imag(CD); 505 r[2][1] = 0.5 * R + 0.5 * imag(CD);
506 r[2][2] = 0.5 * R - 0.5 * imag(CD); 506 r[2][2] = 0.5 * R - 0.5 * imag(CD);
507 507
508 G4complex CE2(c, -ds); 508 G4complex CE2(c, -ds);
509 G4complex CE = std::sqrt(CE2); 509 G4complex CE = std::sqrt(CE2);
510 510
511 r[1][3] = -0.25 * a3 + 0.5 * real(CE); 511 r[1][3] = -0.25 * a3 + 0.5 * real(CE);
512 r[1][4] = -0.25 * a3 - 0.5 * real(CE); 512 r[1][4] = -0.25 * a3 - 0.5 * real(CE);
513 r[2][3] = -0.5 * R + 0.5 * imag(CE); 513 r[2][3] = -0.5 * R + 0.5 * imag(CE);
514 r[2][4] = -0.5 * R - 0.5 * imag(CE); 514 r[2][4] = -0.5 * R - 0.5 * imag(CE);
515 } 515 }
516 } 516 }
517 return 4; 517 return 4;
518 } 518 }
519 519
520 // 520 //
521 // 521 //
522 ////////////////////////////////////////////// 522 //////////////////////////////////////////////////////////////////////////////
523 523