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