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