Unstable Isotopes in NeutronIrradiation

[jackie-k-garcia]

We are using the Isotope class for RadWatch analysis, specifically for NAA. We found that there is a discrepancy between what NNDC says about Eu-151 versus what nucleardata.lu regarding stability; NNDC says Eu-151 is not stable. This is complicating our attempts to use this class for calculations for finding concentration of Eu-151 in our samples because the NeutronIrradiation class currently does not allow for unstable isotopes as the initial isotope.

[jccurtis]

Hi Jackie, can you post some example code for this (probably just a couple of lines...)? It would make it easier to test. We should add this to the pytest suite.

@markbandstra what do you think?

[bplimley]

The IsotopeQuantity class doesn't support activation of an unstable isotope, yet. It just requires someone to do the math for activation concurrent with decay.

As a temporary workaround, for an isotope like Eu-151 which has a very very long half-life, I think you can manually set the IsotopeQuantity object's half life, eu151.half_life = np.inf, before you do the activate call. I don't have a chance to test this right now though.

[alihanks]

Yes, the real issue was that an isotope that effectively has an infinite half-life (and does in at least one nuclear data table) does not according to nndc.

Good suggestion @bplimley, this works. I guess in the long run, before implementing the handling of activation concurrent with decay, it might be worth relaxing the requirement that the half-life be strictly infinite? ~10^25 years seems close enough to infinite to me. :)

@jccurtis I think the current behavior is 'right' in the sense that it does what you would expect and what it should, so no pytest suite additions seem necessary. Our issue was more with the strict requirement of an infinite half-life.

[markbandstra]

From looking at the IsotopeQuantity code, it looks like the main purpose of knowing whether an isotope is stable or unstable is to know whether to deal in units of mass (g) or activity (Bq or uCi).

I think it would be fine to set a large but arbitrary limit for what would be considered stable in this context. We should definitely add a test to the suite to track this issue. I might have time soon to add this fix to PR #62.

Incidentally, it was not known until 10 years ago that Eu-151 is actually unstable!

[markbandstra]

@jackie-k-garcia please pull the latest from the isotope-2 branch and try it out, hopefully it fixes your issue.

[bplimley]

@markbandstra thank you for the quick fix.

Should we leave this issue open until a general solution to activation of unstable isotopes can be implemented?

[markbandstra]

Just to make it clear, the solution I implemented allows the user to define their own limit for what is considered stable versus unstable. The default limit is 1E18 seconds but it can be changed when instantiating IsotopeQuantity with keyword stability. So if anyone encounters an isotope for which this fix does not work, you can change the limit to suit that particular case.

@bplimley It sounds like you mean reworking all of the math to not care whether an isotope is stable or unstable, and that will require some more attention from someone. We can keep this issue open to track progress.

[bplimley]

Yes, I like what you have currently but I was thinking forward to a general solution with reworked math. It is low priority because I don't have a particular use case, it just seems like a good idea.

[alihanks]

It looks like this issue was only addressed in cases where an IsotopeQuantity is used as initial.

if not initial.is_stable:
        raise NotImplementedError(
            'Activation not implemented for a radioactive initial isotope')`

This will still cause a failure if an Isotope is given as initial, where is_stable does not come from some stability threshold the user can specify but rather from:

def is_stable(self):
    """Return a boolean of whether the isotope is stable.
    Returns:
      a boolean.
    """

    df = self._wallet_card()
    data = df['T1/2 (txt)'].tolist()
    assert len(np.unique(data)) == 1
    return 'STABLE' in data

[jccurtis]

@alihanks has a resolution: add a kwarg for the minimum half life which is considered stable in NAA activate method. This will directly check the half life property to determine stability.

[jccurtis]

@alihanks Can this be resolved with the above method?

[jccurtis]

@alihanks is this resolved?

[jccurtis]

@alihanks is this ready for primetime?

[alihanks]

Yeah, this is ready.

[alihanks]

This is still an issue - the Isotope class does not have a good way of handling cases where isotopes are technically unstable but for all practical purposes are stable.

[markbandstra]

Got it, it looks like your fix never got merged in. Sorry! I'll handle that with a PR.