@@ -503,7 +503,11 @@ void PhotonInteraction::compton_doppler(
503503 }
504504
505505 // Sample value on bounded cdf
506- c = prn (seed) * c_max;
506+ if (pz_max > 0.0 ) {
507+ c = prn (seed) * c_max;
508+ } else {
509+ c = 1.0 + (c_max - 1.0 ) * prn (seed);
510+ }
507511
508512 // Determine pz corresponding to sampled cdf value
509513 tensor::View<const double > cdf_shell = profile_cdf_.slice (shell);
@@ -528,29 +532,24 @@ void PhotonInteraction::compton_doppler(
528532 double momentum_sq = std::pow ((pz / FINE_STRUCTURE), 2 );
529533 double f = 1.0 + alpha * (1.0 - mu);
530534 double a = momentum_sq - f * f;
531- double b = 2.0 * E * (f - momentum_sq * mu);
532- c = E * E * (momentum_sq - 1.0 );
533-
534- double quad = std::sqrt (b * b - 4.0 * a * c);
535- double E_out1 = -(b + quad) / (2.0 * a);
536- double E_out2 = -(b - quad) / (2.0 * a);
537-
538- // Determine solution to quadratic equation that is positive
539- if (E_out1 > 0.0 ) {
540- if (E_out2 > 0.0 ) {
541- // If both are positive, pick one at random
542- *E_out = prn (seed) < 0.5 ? E_out1 : E_out2;
543- } else {
544- *E_out = E_out1;
545- }
546- } else {
547- if (E_out2 > 0.0 ) {
548- *E_out = E_out2;
549- } else {
550- // No positive solution -- resample
551- continue ;
552- }
535+ double b = 2.0 * (f - momentum_sq * mu);
536+ c = (momentum_sq - 1.0 );
537+ double quad = b * b - 4.0 * a * c;
538+ if (quad < 0 ) {
539+ write_message (" Negative quad {}" , quad);
540+ continue ;
553541 }
542+ quad = std::sqrt (quad);
543+ double E_out1 = -(b + quad) / (2.0 * a) * E;
544+ double E_out2 = -(b - quad) / (2.0 * a) * E;
545+ // Determine solution to quadratic equation that is positive and minimal
546+ if ((E_out1 < 0.0 ) && (E_out2 < 0.0 ))
547+ continue ;
548+ if ((E_out1 > 0.0 ) && (E_out2 > 0.0 ))
549+ *E_out = std::min (E_out1, E_out2);
550+ else
551+ *E_out = std::max (E_out1, E_out2);
552+
554553 if (prn (seed) * E <= *E_out)
555554 break ;
556555 }
0 commit comments