|
Lines 102-111
G4DensityEffectCalculator::G4DensityEffectCalculator(const G4Material* mat, G4in
Link Here
|
| 102 |
for(G4int i=0; i<nlev; ++i) { |
102 |
for(G4int i=0; i<nlev; ++i) { |
| 103 |
sum += sternf[i]; |
103 |
sum += sternf[i]; |
| 104 |
} |
104 |
} |
| 105 |
sum = (sum > 0.0) ? 1./sum : 0.0; |
105 |
sum += fConductivity; |
|
|
106 |
|
| 107 |
const G4double invsum = (sum > 0.0) ? 1./sum : 0.0; |
| 106 |
for(G4int i=0; i<nlev; ++i) { |
108 |
for(G4int i=0; i<nlev; ++i) { |
| 107 |
sternf[i] *= sum; |
109 |
sternf[i] *= invsum; |
| 108 |
} |
110 |
} |
|
|
111 |
fConductivity *= invsum; |
| 112 |
|
| 109 |
plasmaE = fMaterial->GetIonisation()->GetPlasmaEnergy()/CLHEP::eV; |
113 |
plasmaE = fMaterial->GetIonisation()->GetPlasmaEnergy()/CLHEP::eV; |
| 110 |
meanexcite = fMaterial->GetIonisation()->GetMeanExcitationEnergy()/CLHEP::eV; |
114 |
meanexcite = fMaterial->GetIonisation()->GetMeanExcitationEnergy()/CLHEP::eV; |
| 111 |
} |
115 |
} |
|
Lines 209-224
G4double G4DensityEffectCalculator::FermiDeltaCalculation(G4double x)
Link Here
|
| 209 |
|
213 |
|
| 210 |
// Calculate the Sternheimer adjusted energy levels and parameters l_i given |
214 |
// Calculate the Sternheimer adjusted energy levels and parameters l_i given |
| 211 |
// the Sternheimer parameter rho. |
215 |
// the Sternheimer parameter rho. |
| 212 |
sternrho /= plasmaE; |
|
|
| 213 |
for(G4int i=0; i<nlev; ++i) { |
216 |
for(G4int i=0; i<nlev; ++i) { |
| 214 |
sternEbar[i] = levE[i] * sternrho; |
217 |
sternEbar[i] = levE[i] * sternrho/plasmaE; |
| 215 |
sternl[i] = std::sqrt(gpow->powN(sternEbar[i], 2) + 2./3.*sternf[i]); |
218 |
sternl[i] = std::sqrt(gpow->powN(sternEbar[i], 2) + 2./3.*sternf[i]); |
| 216 |
} |
219 |
} |
| 217 |
|
220 |
|
| 218 |
// Make imphirical initial guess |
221 |
// The derivative of the function we are solving for is strictly |
| 219 |
const G4double sternL = Newton(sternrho, false); |
222 |
// negative for positive (physical) values, so if the value at |
| 220 |
if(sternL > -1.) { |
223 |
// zero is less than zero, it has no solution, and there is no |
| 221 |
return DeltaOnceSolved(sternL); |
224 |
// density effect in the Sternheimer "exact" treatment (which is |
|
|
225 |
// still an approximation). |
| 226 |
if(Ell(0) <= 0) return 0; |
| 227 |
|
| 228 |
// Attempt to find the root from 40 starting points evenly distributed |
| 229 |
// in log space. Trying a single starting point is not sufficient for |
| 230 |
// convergence in most cases. |
| 231 |
for(G4int startLi = -10; startLi < 30; startLi++){ |
| 232 |
const G4double sternL = Newton(gpow->powN(2, startLi), false); |
| 233 |
if(sternL != -1.) { |
| 234 |
return DeltaOnceSolved(sternL); |
| 235 |
} |
| 222 |
} |
236 |
} |
| 223 |
|
237 |
|
| 224 |
return -1.; // Signal the caller to use the Sternheimer approximation, |
238 |
return -1.; // Signal the caller to use the Sternheimer approximation, |
|
Lines 251-257
G4double G4DensityEffectCalculator::Newton(G4double start, G4bool first)
Link Here
|
| 251 |
const G4double del = value/dvalue; |
265 |
const G4double del = value/dvalue; |
| 252 |
lambda -= del; |
266 |
lambda -= del; |
| 253 |
|
267 |
|
| 254 |
const G4double eps = std::abs(del); |
268 |
const G4double eps = std::abs(del/lambda); |
| 255 |
if(eps <= 1.e-12) { |
269 |
if(eps <= 1.e-12) { |
| 256 |
++ngood; |
270 |
++ngood; |
| 257 |
if(ngood == 2) { |
271 |
if(ngood == 2) { |
|
Lines 263-269
G4double G4DensityEffectCalculator::Newton(G4double start, G4bool first)
Link Here
|
| 263 |
} else { |
277 |
} else { |
| 264 |
++nbad; |
278 |
++nbad; |
| 265 |
} |
279 |
} |
| 266 |
if(nbad > maxIter || eps > 1.) { break; } |
280 |
if(nbad > maxIter || std::isnan(value) || std::isinf(value)) { break; } |
| 267 |
} |
281 |
} |
| 268 |
if(fVerbose > 2) { |
282 |
if(fVerbose > 2) { |
| 269 |
G4cout << " Failed to converge last value= " << value |
283 |
G4cout << " Failed to converge last value= " << value |
|
Lines 318-323
G4double G4DensityEffectCalculator::DEll(G4double L)
Link Here
|
| 318 |
ans += sternf[i]/gpow->powN(y + L*L, 2); |
332 |
ans += sternf[i]/gpow->powN(y + L*L, 2); |
| 319 |
} |
333 |
} |
| 320 |
} |
334 |
} |
|
|
335 |
|
| 336 |
ans += fConductivity/gpow->powN(L*L, 2); |
| 337 |
|
| 321 |
ans *= (-2*L); // pulled out of the loop for efficiency |
338 |
ans *= (-2*L); // pulled out of the loop for efficiency |
| 322 |
return ans; |
339 |
return ans; |
| 323 |
} |
340 |
} |
|
Lines 332-343
G4double G4DensityEffectCalculator::Ell(G4double L)
Link Here
|
| 332 |
ans += sternf[i]/(gpow->powN(sternEbar[i], 2) + L*L); |
349 |
ans += sternf[i]/(gpow->powN(sternEbar[i], 2) + L*L); |
| 333 |
} |
350 |
} |
| 334 |
} |
351 |
} |
|
|
352 |
|
| 353 |
if(fConductivity > 0. && L != 0.) { |
| 354 |
ans += fConductivity/(L*L); |
| 355 |
} |
| 356 |
|
| 335 |
ans -= gpow->powZ(10, -2 * sternx); |
357 |
ans -= gpow->powZ(10, -2 * sternx); |
| 336 |
return ans; |
358 |
return ans; |
| 337 |
} |
359 |
} |
| 338 |
|
360 |
|
| 339 |
/** |
361 |
/** |
| 340 |
* Given the Sternheimer parameter l^2 (called 'sternL' here), and that |
362 |
* Given the Sternheimer parameter l (called 'sternL' here), and that |
| 341 |
* the l_i and adjusted energies have been found with SetupFermiDeltaCalc(), |
363 |
* the l_i and adjusted energies have been found with SetupFermiDeltaCalc(), |
| 342 |
* return the value of delta. Helper function for DoFermiDeltaCalc(). |
364 |
* return the value of delta. Helper function for DoFermiDeltaCalc(). |
| 343 |
*/ |
365 |
*/ |
|
Lines 350-355
G4double G4DensityEffectCalculator::DeltaOnceSolved(G4double sternL)
Link Here
|
| 350 |
+ gpow->powN(sternL, 2))/gpow->powN(sternl[i], 2)); |
372 |
+ gpow->powN(sternL, 2))/gpow->powN(sternl[i], 2)); |
| 351 |
} |
373 |
} |
| 352 |
} |
374 |
} |
|
|
375 |
|
| 376 |
// sternl for the conduction electrons is sqrt(fConductivity), with |
| 377 |
// no factor of 2./3 as with the other levels. |
| 378 |
if(fConductivity > 0) |
| 379 |
ans += fConductivity * G4Log((fConductivity |
| 380 |
+ gpow->powN(sternL, 2))/fConductivity); |
| 381 |
|
| 353 |
ans -= gpow->powN(sternL, 2)/(1 + gpow->powZ(10, 2 * sternx)); |
382 |
ans -= gpow->powN(sternL, 2)/(1 + gpow->powZ(10, 2 * sternx)); |
| 354 |
return ans; |
383 |
return ans; |
| 355 |
} |
384 |
} |
