In This Package:

u238.py

Go to the documentation of this file.

#!/usr/bin/env python
'''
238U decay chain
'''

import decay,chain
from data import units


class U238Tests(object):

    _chain = chain.G4DataChain(92,238)

    def __init__(self):
        import random
        self.random = random.Random()
        return

    def isotope(self): return self.chain().head
    def chain(self): return U238Tests._chain

    def setAbundance(self,grams):
        self.isotope().abundance = grams * 6.02e23/238.0

    def meanRate(self):
        return self.isotope().lifetime/self.isotope().abundance

    def singleDecays(self):
        target = 10000
        count = 0
        while count < target:
            count += 1
            rad = self.isotope().decay()
            #print 'Resulting radiation:'
            while rad:
                time,trans = rad
                #print '\t@%.4E: %s'%rad
                rad = trans.final.decay()
                continue
            continue
        for iso in self.chain().groundStates:
            print '%s abundance=%.3f'%(iso,iso.abundance)
        return

    def radiateTransitions(self,dt_trans):
        '''Given list of (dt,transition), radiate the transitions.  If
        a transition results in a excited daughter, immediately radate
        the gamma.  A copy of dt_trans, with any gamma transitions
        appended, is returned.'''

        ret = []

        # radiate all transitions to increase daugher abundance
        for dt,tran in dt_trans:
            #print '%f %s'%(dt,tran)
            tran.radiate()

            ret.append((dt,tran))

            try: # if final state is excited, radiate gamma immediately
                t = tran.final.transitions[0]
                if type(t.radiation) == decay.GammaDecay:
                    ret.append((dt,t))
            except IndexError:
                pass
            continue
        return ret


    def decayChainInterval(self,interval):
        import decay
        dt_trans = []
        for s in self.chain().groundStates:
            if s.abundance == 0.0: continue
            dks = s.decayInterval(interval)
            #print 'Got %d %s decays in %f seconds'%(len(dks),s,interval/units.second)
            if not dks: continue
            dt_trans += dks
            continue
        return self.radiateTransitions(dt_trans)

    def trackAbundances(self,cycles,interval):
        abs = {}
        for s in self.chain().groundStates:
            abs[s] = []
        ntrans = []
        while cycles:
            cycles -= 1
            dks = self.decayChainInterval(interval)
            ntrans.append(len(dks))
            for s in self.chain().groundStates:
                abs[s].append(s.abundance)
                continue
            continue
        return abs,ntrans

    def plotAbundances(self,abs):
        import matplotlib.pyplot as plt
        plt.figure(1)

        nplots = sum([s.abundance > 0.0 and s.symbol != '238U' for s in abs.keys()])

        nstates = len(abs)
        n = 0
        for s,ab in abs.iteritems():
            if s.abundance == 0.0 or s.symbol == '238U': 
                print 'Skipping %s'%s
                continue
            n += 1
            end = len(ab)
            print 'Plotting %s with final abundance %f from 0 to %d'%(s,s.abundance,end)
            plt.subplot(nplots,1,n)
            plt.plot(range(0,end),ab,'r')
            plt.title(s.symbol)
            continue
        plt.show()
        return

    def setSecularEquilibrium(self):
        '''Sets 8 isotopes starting subchains in secular equilibrium
        with the 238U.  Returns a list with 238U and the 8 daughters'''
        daughters = ['234Th','234U','230Th','226Ra','222Rn','210Pb','210Bi','210Po']
        u238 = self.isotope()
        ret = [u238]
        for s in self.chain().groundStates:
            if not s.symbol in daughters: continue
            s.abundance = u238.abundance*s.halflife/u238.halflife
            ret.append(s)
        print '\n'.join(map(str,ret))
        return ret

    def generateDecays(self,states,ndecays):
        '''Generate ndecays starting from the given list of parent
        states.  Parent states are considered independent in the sense
        that a decay subchain is fortced to occur until either a
        stable state is reached or a state in the list of states is
        reached.  A list of the resulting radiations will be returned.
        '''
        totalrate = 0
        for s in states:
            totalrate += s.abundance/s.halflife
            continue

        decays = []

        while ndecays:
            ndecays -= 1

            # choose which state gets to do the honors
            r = self.random.uniform(0,totalrate)
            goal = 0
            for s in states:
                goal += s.abundance/s.halflife
                if r < goal: 
                    goal = s
                    break
                continue
        
            print 'Starting %s (%f/%f)'%(goal,r,totalrate)

            # decay until stable or daughter is in list of states
            while True:
                rad = goal.decay()
                decays.append(rad)
                trans = rad[1]
                trans.radiate()
                final = trans.final
                if final in states or final.symbol == '206Pb':
                    break
                goal = final
                continue
            continue
        return decays
    pass

def seceq():
    u238 = U238Tests()
    u238.setAbundance(grams=20.0e-6)
    abs,ntrans = u238.trackAbundances(10000,1000)
    u238.plotAbundances(abs)

def decays():
    u238 = U238Tests()
    u238.setAbundance(grams=20.0e-6)
    heads = u238.setSecularEquilibrium()
    print '\n'.join(map(str,heads))
    dt_trans = u238.generateDecays(heads,10)
    print '%d transitions:'%len(dt_trans)
    for dt,trans in dt_trans:
        print '@%8.3f %s'%(dt,trans)
    return u238

if '__main__' == __name__:
    #seceq()
    decays()


    

---