Repository navigation
Conversation
| x, y = (b-a)*t, (c-a)*t | ||
| if y > 0.01: | ||
| # Use the largest separation in the denominator. | ||
| return (_two_rate_integral(t, a, b)-_two_rate_integral(t, b, c))/(c-a) |
There was a problem hiding this comment.
I verified that 883 matches the summation in the Bateman equation for N=3 Wikipedia, which matches the summation in column T.
With (a, b, c) sorted, this should limit cancellation errors from b-a and c-b in the denominators, so this looks good.
| homogeneous = y*homogeneous + x_power | ||
| factorial *= order+2 | ||
| coefficient += (-1)**order * homogeneous/factorial | ||
| return t*(t*exp(-a*t))*coefficient |
There was a problem hiding this comment.
I didn't verify these equations directly, but the tests show that they are sufficient when a and c are comparable.
| decay = exp(-min(a, b)*t) | ||
| if x < 1e-4: | ||
| return t*decay*(-expm1(-x)/x if x else 1.) | ||
| return decay*(-expm1(-x))/delta |
There was a problem hiding this comment.
Line 874 for x != 0 is equal to line 875 (the t cancels the 1/t in 1/x), so just need to check x==0.
Agreed that for the x=0 case the result should be λt exp(-at) = λt exp(-bt). This would have required an extra condition on the original line 841.
The expression on 875 is equivalent to what I had in 841, so I don't expect any improvement in accuracy.
Rewrite using the sorting technique from _three_rate_integral:
def _two_rate_integral(t: float, a: float, b: float) -> float:
"""Convolve exp(-a*t) and exp(-b*t), including coincident rates."""
if a == b:
return t * exp(-a*t)
if a > b:
a, b = b, a
return -exp(-a*t) * expm1((a-b)*t) / (b - a)There was a problem hiding this comment.
Restored the spreadsheet equations and their mapping. The two-rate helper now sorts its inputs; it divides before applying the decay and retains the limit for subnormal (b-a)*t. Tests cover both underflow cases and input orders.
| def _two_rate_integral(t: float, a: float, b: float) -> float: | ||
| """Convolve exp(-a*t) and exp(-b*t), including coincident rates.""" | ||
| if a == b: | ||
| return t*exp(-a*t) |
There was a problem hiding this comment.
Drop this if-block. The condition is captured by 901 below.
There was a problem hiding this comment.
Removed the redundant equality branch; coincident rates use the existing zero-product limit. The source and installed-wheel suites pass.
Sample.calculate_activation()can divide by zero when capture and removal rates coincide. For a 1 g cobalt sample with a flux of117540904144107.45and a one-hour exposure, the current code raisesZeroDivisionError; the next representable flux produces a Co-60 activity about 10.3% above the finite two-state limit.Evaluate the single- and double-capture exponential convolutions with stable divided differences, including repeated-rate limits. Add real Sample API regressions at equal and adjacent rates, plus independent high-precision checks of the three-rate expression. Update one existing Co-61 reference to its 100-digit value while keeping the
1e-12tolerance.Validation on Windows/Python 3.13: 93 configured source checks passed, sdist and wheel built, and 60 installed-wheel tests passed. The new API regressions fail on unchanged master (six failures, three zero-exposure passes). This change leaves the beta/post-irradiation chain calculation unchanged.