data.py

Go to the documentation of this file.

#!/usr/bin/env python

class units:

    second = 1.
    minute = 60*second
    hour = 60*minute
    day = 24*hour
    year = 365*day

    MeV = 1.
    eV = 1.0e-6*MeV
    keV = 1.0e-3*MeV
    GeV = 1.0e+3*MeV

    electron_mass_c2 = 0.51099906 * MeV
    pass


class Isotope(object):
    '''A nuclide'''

    def __init__(self,Z=0,A=0,halflife=0.0,
                 level=0.0,abundance=0):
        import elements
        self.name = '%s(%d)'%(elements.name(Z),A)
        self.symbol = '%d%s'%(A,elements.symbol(Z))
        self.A=A
        self.Z=Z
        self.halflife=halflife
        self.level=level
        self.abundance=abundance
        self.transitions=[]
        return

    def __getattr__(self,name):
        if name == 'N': return self.A - self.Z
        if name == 'lifetime':
            import math
            return self.halflife/math.log(2.0)
        raise AttributeError,'no such attribute ' + name

    def __str__(self):
        return '%s (%.3f)'%(self.symbol,self.level/units.MeV)

    def decayInterval(self,interval):
        '''Decay a number of isotopes drawn from the mean for the
        given interval.  Returns a list of tuples (dt,tranition)'''
        from numpy import random
        rate = self.abundance/self.lifetime
        mean = rate*interval
        ndecays = random.poisson(mean)
        #print "Got %d from %f %s"%(ndecays,mean,self)
        ret = []

        # reduce abundance by hand to avoid subtracting by small numbers
        self.abundance -= ndecays

        # generate decay products
        while ndecays:
            ndecays -= 1
            dt,trans = self.decay(False)
            dt = random.uniform(0,interval)
            ret.append((dt,trans))
            continue

        # sorts by time
        ret.sort()

        return ret

    def decay(self,reduceAbundance = True):
        'Decay this isotope, return [time,Transition].  Does NOT radiate() the Transition'
        # t(n) = t(n-1) * exp(-t/tL)
        # u = uniform()
        # dt = ((-1*math.log(u)) * lifetime)
        from numpy import random
        import math

        total = sum([t.fraction for t in self.transitions])
        u = random.uniform(0.0,total)
        total = 0.0
        trans = None
        #print 'Checking',len(self.transitions),'transitions with u',u
        for t in self.transitions:
            total += t.fraction            
            if u < total:
                trans  = t
                break
            continue
        if trans is None: return ()
        
        u = random.uniform(0,1)
        dt = ((-1*math.log(u)) * self.lifetime/self.abundance)
        if reduceAbundance:
            self.abundance -= 1
        #trans.radiate()
        return  (dt,trans)

    pass
    

        

class Transition(object):
    '''Transition from an initial to final nuclide states radioactive
emission'''

    def __init__(self,initial=None,final=None,radiation=None,
                 lifetime=0.0,fraction=1.0):
        self.initial=initial
        self.final=final
        self.radiation=radiation
        self.fraction=fraction
        self.count = 0
        self.initial.transitions.append(self)
        #print self
        return

    def __str__(self):
        return 'Transition[%7.3f%%]: %s --> %s + %s'%(self.fraction*100,self.initial,self.final,self.radiation)

    def radiate(self):
        self.count += 1
        self.final.abundance += 1
        return

    pass

class PrettyPrint(object):

    _tab = '  '

    def __init__(self,iso,depth=0):
        self.iso = iso
        self.depth = depth
        return

    def tab(self): return PrettyPrint._tab*self.depth

    def __str__(self):
        import decay
        ret = [self.tab()+str(self.iso)]

        iso2process = []
        for t in self.iso.transitions: # assume only be non-GammaDecay

            # Add this transition
            #print self.tab()+str(t)
            ret.append(self.tab()+str(t))

            # if daughter not stable
            if t.final.transitions: 
                ft = t.final.transitions[0] 
                # and if daugther decays by gamma
                if type(ft.radiation) == decay.GammaDecay: 
                    # save string rep and ground state for later processing
                    ret.append(self.tab()+PrettyPrint._tab+str(ft))
                    if not ft.final in iso2process:
                        iso2process.append(ft.final)
                    continue

            # otherwise, save for recursion
            if not t.final in iso2process: 
                iso2process.append(t.final)
            continue
        for iso in iso2process:
            pp = PrettyPrint(iso,self.depth+1)
            ret.append(str(pp))

        return '\n'.join(ret)
                       


if __name__ == '__main__':
    plot_64Cu()