""" Collection of routines for isotopes, most prominently the Ion class. """ # contains all the isotope routines from __future__ import division from __future__ import print_function from __future__ import with_statement import string import re import os import copy import numpy as np from matplotlib.cbook import is_numlike from collections import Iterable import time import textwrap import os.path import time2human import datetime from logged import Logged from utils import cachedmethod, project, OutFile, stuple from version2human import version2human # Would b better if this was a class so we could write # elements.H elements = ( 'nt', 'h', 'he', 'li', 'be', 'b', 'c', 'n', 'o', 'f', 'ne', 'na', 'mg', 'al', 'si', 'p', 's', 'cl', 'ar', 'k', 'ca', 'sc', 'ti', 'v', 'cr', 'mn', 'fe', 'co', 'ni', 'cu', 'zn', 'ga', 'ge', 'as', 'se', 'br', 'kr', 'rb', 'sr', 'y', 'zr', 'nb', 'mo', 'tc', 'ru', 'rh', 'pd', 'ag', 'cd', 'in', 'sn', 'sb', 'te', 'i', 'xe', 'cs', 'ba', 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu', 'gd', 'tb', 'dy', 'ho', 'er', 'tm', 'yb', 'lu', 'hf', 'ta', 'w', 're', 'os', 'ir', 'pt', 'au', 'hg', 'tl', 'pb', 'bi', 'po', 'at', 'rn', 'fr', 'ra', 'ac', 'th', 'pa', 'u', 'np', 'pu', 'am', 'cm', 'bk', 'cf', 'es', 'fm', 'md', 'no', 'lr', 'rf', 'db', 'sg', 'bh', 'hs', 'mt', 'ds', 'rg', 'cn', 'ut', 'uq', 'up', 'uh', 'us', 'uo') # should be Uut Uuq Uup Uuh Uus Uuo but I will not ue 3 letters # as this breaks too many things Elements = [element.capitalize() for element in elements] Elements[0] = 'n' Elements = tuple(Elements) z2Name = ( "Neutron", "Hydrogen", "Helium", "Lithium", "Beryllium", "Boron", "Carbon", "Nitrogen", "Oxygen", "Fluorine", "Neon", "Sodium", "Magnesium", "Aluminium", "Silicon", "Phosphorus", "Sulfur", "Chlorine", "Argon", "Potassium", "Calcium", "Scandium", "Titanium", "Vanadium", "Chromium", "Manganese", "Iron", "Cobalt", "Nickel", "Copper", "Zinc", "Gallium", "Germanium", "Arsenic", "Selenium", "Bromine", "Krypton", "Rubidium", "Strontium", "Yttrium", "Zirconium", "Niobium", "Molybdenum", "Technetium", "Ruthenium", "Rhodium", "Palladium", "Silver", "Cadmium", "Indium", "Tin", "Antimony", "Tellurium", "Iodine", "Xenon", "Caesium", "Barium", "Lanthanum", "Cerium", "Praseodymium", "Neodymium", "Promethium", "Samarium", "Europium", "Gadolinium", "Terbium", "Dysprosium", "Holmium", "Erbium", "Thulium", "Ytterbium", "Lutetium", "Hafnium", "Tantalum", "Tungsten", "Rhenium", "Osmium", "Iridium", "Platinum", "Gold", "Mercury", "Thallium", "Lead", "Bismuth", "Polonium", "Astatine", "Radon", "Francium", "Radium", "Actinium", "Thorium", "Protactinium", "Uranium", "Neptunium", "Plutonium", "Americium", "Curium", "Berkelium", "Californium", "Einsteinium", "Fermium", "Mendelevium", "Nobelium", "Lawrencium", "Rutherfordium", "Dubnium", "Seaborgium", "Bohrium", "Hassium", "Meitnerium", "Darmstadtium", "Roentgenium", "Copernicium", "(Ununtrium)", "(Ununquadium)", "(Ununpentium)", "(Ununhexium)", "(Ununseptium)", "(Ununoctium)") z2name = tuple((name.lower() for name in z2Name)) # el2z provides easy convesion from element symbol to charge number el2z = {el:i for i,el in enumerate(elements)} el2z.update({el:i for i,el in enumerate(Elements[1:],start=1)}) def el2name(sym): return z2name[el2z[sym]] class Ion(object): """ This is a generic class for 'ions'. It is designed to handel elements, isotioes, isomers, and isobars. Isomers have an energy level appended by 'e' or 'E'; if the level is '0' then it should not be printed. It provides also an "index" field (idx) that can be used for sorting, but this may be obsolte for external use as it provides a __cmp__ operator. __add__ and __sub__ should help with networks. R([[I],O]) allows to formulate reactions will return ground state nuclei Use (solely) A = ... for 'isobars' Z = ... for 'elements' N = ... for 'isotones' To specify isotopes you need to supply at least 2 of A, Z, N g = ( 0, 0, 0, 0) e- = ( 1,-1, 0, 0) e+ = ( 1,+1, 0, 0) void = (-1, 0, 0, 0) TODO: ... maybe always put E1 as 'm' ... maybe not ... ... destinguish isotopes and isomers - Al26g is not Al26 """ # in the future we want to use these # should flags be first or last? # maybe last? AMAX = 512 ZMAX = 128 EMAX = 256 FMAX = 128 # flags - 7 bit # leave the sign bit for int64 EMUL = 1 AMUL = EMAX ZMUL = AMUL * AMAX FMUL = ZMUL * ZMAX BMUL = FMUL * FMAX # what lies Beyond ... (extra information for subclasses) # if these are used, the index no longer fits into int64 # bit 0-2 F_BOSON = 0 F_ISOMER = 1 F_ISOTOPE = 2 F_ELEMENT = 3 F_ISOBAR = 4 F_ISOTONE = 5 F_HADRON = 6 # currently not used F_LEPTON = 7 F_GROUP_MASK = 7 # bit 3 F_REACTION = 8 # bit 4 F_KEPION = 16 Excited = 'E' excited = 'e' isomer = 'm' ground = 'g' excite = (Excited, excited, isomer) VOID = (0, -1, 0, 0, 0) VOID_IDX = -1 * FMUL VOID_STRING = '-' GAMMA = (0, F_BOSON, 0, 0, +1) special = {0: 'nt', 1 * EMUL: 'g', (F_LEPTON * FMUL) + (+1 * ZMUL) : 'e+', (F_LEPTON * FMUL) + (-1 * ZMUL) : 'e-', VOID_IDX : VOID_STRING} def __init__(self, name = None, Z = None, A = None, E = None, F = None, B = None, idx = None, N = None, element = True): """ Initialize Isotope. """ self.F = self.F_ISOMER if F == None else F self.Z = 0 if Z == None else Z self.A = 0 if A == None else A self.E = 0 if E == None else E self.B = 0 if B == None else B if N != None: if Z != None: self.A = Z + N assert (A == self.A) or (A == None), "Conflict." elif A != None: self.Z = A - N self.A = N else: self.F = self.F_ISOTONE self.A = N assert (self.F == F) or (F == None), "Conflict." assert self.E == 0, "Not sure what this would be." assert self.A >= 0, "Not sure what this would be." if Z != None: if (N == None) and (A == None): self.F = self.F_ELEMENT assert (self.F == F) or (F == None), "Conflict." assert self.E == 0, "Not sure what this would be." if A != None: if (Z == None) and (N == None): self.F = self.F_ISOBAR assert (self.F == F) or (F == None), "Conflict." assert self.E == 0, "Not sure what this would be." assert self.A > 0, "Not sure what this would be." if (E == None) or (E == -1): if self.F & self.F_GROUP_MASK == self.F_ISOMER: self.F = self.F_ISOTOPE self.E == 0 s = name if isinstance(s, Ion): self.B, self.F, self.Z, self.A, self.E = s.tuple() # conversion to Ion class # or could add keyword if not issubclass(type(self),type(s)): self.B = 0 self.F = 0 s = None if (type(s) == type(0)) and (idx == None): if 0 < s < len(Elements): self.Z = s self.A = Z if Z is not None else 0 self.E = A if A is not None else 0 self.F = self.F_ISOMER if A == None: self.F = self.F_ISOTOPE if Z == None: self.F = self.F_ELEMENT else: idx = s s = None if type(s) == type(()): if len(s) == 1: self.Z = s[0] self.F = self.F_ELEMENT if len(s) == 2: self.Z = s[0] self.A = s[1] self.F = self.F_ISOTOPE if len(s) == 3: self.Z = s[0] self.A = s[1] self.E = s[2] self.F = self.F_ISOMER if len(s) == 4: self.F = s[0] self.Z = s[1] self.A = s[2] self.E = s[3] if len(s) == 5: self.B = s[0] self.F = s[1] self.Z = s[2] self.A = s[3] self.E = s[4] s = None if s != None: self.B, self.F, self.Z, self.A, self.E = self.ion2bfzae(s, element = element) if idx != None: self.B, self.F, self.Z, self.A, self.E = self.idx2bfzae(idx) self.idx = self.bfzae2idx(F = self.F, Z = self.Z, A = self.A, E = self.E, B = self.B) self.N = self.get_N() super(Ion, self).__init__() # super(self.__class__,self).__init__() # py3: super().__init__() @classmethod def min_iso(cls, Z): if issubclass(type(Z),cls): Z = Z.Z try: Z = int(Z) except: Z = el2z[Z] return cls(Z = Z, N = 0) @classmethod def max_iso(cls, Z): if issubclass(type(Z),cls): Z = Z.Z try: Z = int(Z) except: Z = el2z[Z] return cls(Z = Z, A = cls.AMAX -1) @classmethod def element_slice(cls, Z): return slice(cls.min_iso(Z), cls.max_iso(Z)) @classmethod def bfzae2idx(cls, F = F_ISOTOPE, Z = 0, A = 0, E = 0, B = 0): """ Compute index from F, Z, A, E. """ assert ((F == -1) and (A == 0) and (Z == 0) and (E == 0) and (B == 0)) \ or ((F >= 0) and (A >= 0) and (E >= 0) and (B >= 0)) \ or (F >= cls.F_REACTION), "Wrong numbers." return (F * cls.FMUL + E * cls.EMUL + Z * cls.ZMUL + A * cls.AMUL + B * cls.BMUL) @classmethod def idx2bfzae(cls, idx): """ Return things that are difined within basic Ion class. """ aidx = abs(idx) f = (aidx // cls.FMUL) % cls.FMAX z = (aidx // cls.ZMUL) % cls.ZMAX a = (aidx // cls.AMUL) % cls.AMAX e = (aidx // cls.EMUL) % cls.EMAX b = (aidx // cls.BMUL) if z > 0: z *= cmp(idx,0) elif idx == cls.VOID_IDX: b, f, z, a, e = cls.VOID else: e *= cmp(idx,0) return b, f, z, a, e @classmethod def ion2idx(cls,s=''): return cls.bfzae2idx(cls.ion2bfzae(s)) def update_idx(self): """ Update idx assuming B,F,A,Z,E are all current """ self.idx = (self.F * self.FMUL + self.E * self.EMUL + self.Z * self.ZMUL + self.A * self.AMUL + self.B * self.BMUL ) def get_N(self): N = self.A - self.Z if self.F in (self.F_BOSON, self.F_ELEMENT, self.F_ISOBAR, self.F_HADRON, self.F_LEPTON): N = 0 return N def tuple(self): return (self.B, self.F, self.Z, self.A, self.E ) def dict(self): return { 'F': self.F, 'Z': self.Z, 'A': self.A, 'E': self.E, 'B' : self.B } def index(self): return self.idx def index_AZE(self): return ( self.E * self.EMUL + self.Z * self.ZMUL + self.A * self.AMUL ) def get_idx(self): return self.idx def get_Z(self): return self.Z def get_A(self): return self.A def get_E(self): return self.E def get_F(self): return self.F def get_B(self): return self.B def Ye(self): if self.Z < 0 or self.A == 0: return -1 return self.Z / self.A def mu(self): if self.Z < 0: return -1 return self.A / (self.Z + 1) def mue(self): if self.Z <= 0: return -1 return self.A / self.Z def Name(self, align = 1, width = 0, upcase = True): """ Return capitalized ion name; same as name otherwise. """ return self.name(align = align, width = width, upcase = upcase) def name(self, align = 1, width = 0, upcase = False): """ Return ion name. align [+1]: -1: left 0: center +1: right (default) width [0]: size of field upcase [False]: capitalize first letter """ s = self._name(upcase = upcase) if width > 0: if align == -1: s = s.ljust(width) elif align == 1: s = s.rjust(width) else: s = s.center(width) return s def element(self): return Ion(Z = self.Z) def isobar(self): return Ion(A = self.A) def isotone(self): return Ion(N = self.N) def isotope(self): return Ion(A = self.A, Z = self.Z) def element_symbol(self, upcase = True): s = '' if self.Z >= 0: if upcase: s = Elements[self.Z] else: s = elements[self.Z] return s def LaTeX(self, math = False): if self.is_element(): el = Elements[self.Z] if math: return r"\mathrm{"+el+"}" else: return el a = str(self.A) el = Elements[self.Z] if math: return r"^{"+a+"}\mathrm{"+el+"}" else: return r'$^{'+a+'}$'+el return s @classmethod def isomer_name(cls, E): if E == 0: return cls.ground elif E == 1: return "{:s}".format(cls.isomer) else: return "{:s}{:d}".format(cls.isomer, E) def mpl(self): if self.is_element(): el = Elements[self.Z] return r"$\mathsf{{{:}}}$".format(el) if self.is_isobar(): return r"$\mathsf{{{:d}}}$".format(self.A) if self.is_isotone(): return r"$\mathsf{{{:d}}}$".format(self.N) if self.is_isotope(): el = Elements[self.Z] a = str(self.A) return r"$\mathsf{{^{{{:}\!}}{:}}}$".format(a,el) if self.is_isomer(): el = Elements[self.Z] a = str(self.A) m = self.isomer_name(self.E) return r"$\mathsf{{^{{{:}\!}}{:}^{{{:}\!}}}}$".format(a,el,m) return '*' def LaTeX_Table(self, math = False, nozero = False, protons = False, manual = False, width = 11, align = 0 ): a = str(self.A) el= Element[self.Z] if manual: s = '$^{'+a+'}$&'+el else: s = "\I{"+a+"}{"+el+"}" if nozero and idx == VOID_IDX: s = '--' if align == -1: s = s.ljust(width) elif align == 1: s = s.rjust(width) else: s = s.center(width) if math and not manual: s = '$'+s+'$' return s def __str__(self): return self._name(upcase = True) def _name(self, upcase = True): """ Convert normal isotopes and isomers to string. For isomers use 'g', 'm', 'm2', m3', ... """ s = '' if self.Z >= 0: if upcase: s = Elements[self.Z] else: s = elements[self.Z] if self.F & self.F_GROUP_MASK == self.F_ISOBAR: s = 'A:' if self.F & self.F_GROUP_MASK == self.F_ISOTONE: s = 'N:' if self.A != 0 or (self.F & self.F_GROUP_MASK == self.F_ISOTONE): s += "{:d}".format(self.A) if self.A == 1 and self.Z == 0 and self.F == 0: s='n' if self.F & self.F_GROUP_MASK == self.F_ISOMER: s += self.isomenr_name(self.E) s = self.special.get(self.idx, s) return s def __cmp__(self,x): if not isinstance(x,Ion): x = Ion(x) return cmp(self.idx, x.idx) @classmethod def ion2bfzae(cls,s='', element = True): sz, sa, se = cls._decompose(s, element = element) valid = True f = cls.F_ISOMER b = 0 try: a = int(sa) except ValueError: a = 0 if sa != '': valid = False try: z = Elements.index(sz) except ValueError: try: z = elements.index(sz) except ValueError: try: z = z2Name.index(sz.strip('-')) except ValueError: try: z = z2name.index(sz.strip('-')) except ValueError: z=0 if sz != '': valid = False if sz in ('n','N') and a > 1: z = 7 if sz in ('p','P') and a > 1: z = 15 if sz == 'e-' and sa == '' and se == '': z = -1 e = 0 f = cls.F_LEPTON valid = True if (sz == 'e+') and sa == '' and se == '': z = +1 e = 0 f = cls.F_LEPTON valid = True try: e = int(se) except ValueError: e = 0 f = cls.F_ISOTOPE if se != '': valid = False if f == cls.F_ISOTOPE and z > 0 and a == 0: f = cls.F_ELEMENT if not valid: b, f, z, a, e = cls.VOID return b, f, z, a, e @classmethod def _decompose(cls,s='',element = True): """ This routine TRIES to decompose string in 1) element symbol 2) mass number 3) excited state number All are returned as strings. The optional parameter 'element' enforces that 'p' is taken as *element* potassium (not H1) 'n' is taken as *element* nitrogen (not neutron) """ name=string.strip(s) n=len(name) el = '' a = '' e = '' # get numbers n=re.split("\D| ",s) while '' in n: n.remove('') # get strings c=re.split("\d| ",s) while '' in c: c.remove('') if len(n) > 0: a = n[0] if len(n) > 1: e = n[1] if len(n) > 2: raise ValueError("Can't understand isotope") if len(c) > 0: el = c[0] if len(el) > 0: if el[-1] in cls.excite and len(c) == 1 and len(n) == 2: c.append(el[-1]) el = el[:-1] if len(c) == 2: if c[1] == 'g': e = '0' if c[1] == 'm' and len(n) == 1: e = '1' if s == 'a' and a == '': el = 'He' a = '4' # this is a possible conflict with potassium if (not element) and s in ('p',): el = 'H' a = '1' if el in ('d','D'): el = 'H' a = '2' if el in ('t','T'): el = 'H' a = '3' if (element) and s == 'nt': el = 'nt' if s in ('nt','nt1'): el = 'nt' a = '1' if s in ('pn','pn1'): el = 'H' a = '1' if (not element) and s in ('n','n1','1n'): el = 'nt' a = '1' # if (not element) and s in ('p','p1','1p'): # el = 'H' # a = '1' if (element) and s == 'n': el = 'N' if (element) and s == 'p': el = 'P' if s == 'g': el = '' a = '' e = '1' if (s.lower() in ('e-','b-','bd','pc')): el = 'e-' if ((s.lower() in ('e+','b+','ec')) or ((not elements) and (s.lower() == 'pd'))): el = 'e+' el = el.strip() a = a.strip() e = e.strip() return el,a,e def dexcite(self): self.E = 0 self.idx = self.bfzae2idx(F = self.F, Z = self.Z, A = self.A, E = self.E, B = self.B) def dexcited(self): return Ion(F = self.F, A = self.A, Z = self.Z, B = self.B) def __add__(self, x): if not isinstance(x,Ion): x=self.__class__(x) return self.__class__(Z = self.Z + x.Z, A = self.A + x.A) __radd__ = __add__ def __sub__(self,x): if not isinstance(x,Ion): x=self.__class__(x) return self.__class__(Z = self.Z - x.Z, A = self.A - x.A) def __rsub__(self,x): if not isinstance(x,Ion): x=self.__class__(x) return self.__class__(Z = x.Z - self.Z, A = x.A - self.A) def __mul__(self,x): return self.__class__(Z = self.Z * x, A = self.A * x, F = self.F, E = self.E) __rmul__ = __mul__ def __rmul__(self,x): return self.__mul__(x) def __neg__(self): return self.__class__(Z = - self.Z, A = - self.A, E = self.E, F = self.F, B = self.B) def __nonzero__(self): return self.idx != self.VOID_IDX def __hash__(self): return self.idx def __len__(self): return 1 def __repr__(self): return self.__class__.__name__ + "('"+self.Name()+"')" def R(self,I=None,O=None): """ Reaction rate function. If only one parameter is given it is interpreted as 'outgoing' If no parameter is provided, just a deep copy of self is returned """ if O == None: if I == None: return self + '' return self - I return (self + I) - O def __getstate__(self): # need to add version number # should save index instead # print('xxx_save') return self.tuple() def __setstate__(self,x): # need to add check of version number # print('xxx_load') self.__init__(x) # this does not appear to be called for a new-style class ever. # def __getinitargs__(self): # print 'xxx_old' # return (self.tuple(),) # the following could work as well, but does appear less flexible # def __getnewargs__(self): # print 'xxx_new' # return (self.tuple(),) def element_name(self): return z2name[self.Z] def is_lepton(self): return self.F & self.F_GROUP_MASK == self.F_LEPTON def is_isotope(self): return self.F & self.F_GROUP_MASK == self.F_ISOTOPE def is_isobar(self): return self.F & self.F_GROUP_MASK == self.F_ISOBAR def is_isotone(self): return self.F & self.F_GROUP_MASK == self.F_ISOTONE def is_element(self): return self.F & self.F_GROUP_MASK == self.F_ELEMENT def is_isomer(self): return self.F & self.F_GROUP_MASK == self.F_ISOMER def is_hadron(self): return self.F & self.F_GROUP_MASK == self.F_HADRON ufunc_A = np.frompyfunc(lambda x: x.A, 1, 1) ufunc_Z = np.frompyfunc(lambda x: x.Z, 1, 1) ufunc_N = np.frompyfunc(lambda x: x.N, 1, 1) ufunc_E = np.frompyfunc(lambda x: x.E, 1, 1) ufunc_F = np.frompyfunc(lambda x: x.F, 1, 1) ufunc_B = np.frompyfunc(lambda x: x.B, 1, 1) ufunc_is_lepton = np.frompyfunc(lambda x: x.is_lepton() , 1, 1) ufunc_is_isotope = np.frompyfunc(lambda x: x.is_isotope(), 1, 1) ufunc_is_isobar = np.frompyfunc(lambda x: x.is_isobar() , 1, 1) ufunc_is_isotone = np.frompyfunc(lambda x: x.is_isotone(), 1, 1) ufunc_is_element = np.frompyfunc(lambda x: x.is_element(), 1, 1) ufunc_is_isomer = np.frompyfunc(lambda x: x.is_isomer() , 1, 1) ufunc_is_hadron = np.frompyfunc(lambda x: x.is_hadron() , 1, 1) def type(self): if self.is_lepton(): return 'lepton' if self.is_isotope(): return 'isotope' if self.is_isobar(): return 'isobar' if self.is_isotone(): return 'isotone' if self.is_element(): return 'element' if self.is_isomer(): return 'isomer' if self.is_hadron(): return 'hadron' return 'unkown' GAMMA = Ion(Ion.GAMMA) VOID = Ion(Ion.VOID) NEUTRON = Ion(Z=0,A=1) PROTON = Ion(Z=1,A=1) ALPHA = Ion(Z=2,A=4) # some convenience functions def Element(Z): return Ion(Z = Z) def Isotone(N): return Ion(N = N) def Isobar(A): return Ion(A = A) def Isotope(A, Z): return Ion(A = A, Z = Z) def Isomer(A, Z, E = 0): return Ion(A = A, Z = Z, E = E) def isoarr(arr): return [Ion(a) for a in arr] class KepIon(Ion): """ Special Ions for KEPLER APPOX/NSE/QSE networks. TODO: need to overwrite add etc. to make sure nothing bad happens ...? """ ionlist={ 'nt1' : ( 1, Ion.F_KEPION, 0, 1, 0), 'h1' : ( 2, Ion.F_KEPION, 1, 1, 0), 'pn1' : ( 3, Ion.F_KEPION, 1, 1, 0), 'he3' : ( 4, Ion.F_KEPION, 2, 3, 0), 'he4' : ( 5, Ion.F_KEPION, 2, 4, 0), 'c12' : ( 6, Ion.F_KEPION, 6,12, 0), 'n14' : ( 7, Ion.F_KEPION, 7,14, 0), 'o16' : ( 8, Ion.F_KEPION, 8,16, 0), 'ne20' : ( 9, Ion.F_KEPION,10,20, 0), 'mg24' : (10, Ion.F_KEPION,12,24, 0), 'si28' : (11, Ion.F_KEPION,14,28, 0), 's32' : (12, Ion.F_KEPION,16,32, 0), 'ar36' : (13, Ion.F_KEPION,18,36, 0), 'ca40' : (14, Ion.F_KEPION,20,40, 0), 'ti44' : (15, Ion.F_KEPION,22,44, 0), 'cr48' : (16, Ion.F_KEPION,24,48, 0), 'fe52' : (17, Ion.F_KEPION,26,52, 0), 'fe54' : (18, Ion.F_KEPION,26,54, 0), 'ni56' : (19, Ion.F_KEPION,28,56, 0), 'ye' : (30, Ion.F_KEPION, 0, 0, 0), 'yq' : (31, Ion.F_KEPION, 0, 0, 0), 'eb0' : (32, Ion.F_KEPION, 0, 0, 0), 'yf' : (33, Ion.F_KEPION, 0, 0, 0), 'fe56' : (34, Ion.F_KEPION,26,56, 0), "'fe'" : (35, Ion.F_KEPION,26,56, 0)} @classmethod def approx_ion_names(cls): """ return names of APPROX19 isotopes """ return [name for name,ion in sorted(cls.ionlist.items(), key = lambda x: x[1][0]) if ion[0] <= 20] def __init__(self, *args, **kwargs): if len(args) == 1 and len(kwargs) == 0: ion = args[0] if isinstance(ion, Ion) and ion != VOID: for n,t in self.ionlist.iteritems(): if (t[2] == ion.Z and t[3] == ion.A and t[4] == ion.E and ((ion.F == Ion.F_KEPION and ion.B == t[0]) or ion.F != Ion.F_KEPION)): name = n super(KepIon, self).__init__(name) return super(KepIon, self).__init__(self.VOID) return super(KepIon, self).__init__(*args, **kwargs) @classmethod def ion2bfzae(cls,s='',element = None): s = s.strip() if s.startswith('ion'): s=s[3:] return cls.ionlist.get(s.lower(), cls.VOID) @classmethod def ionn2bfzae(cls, ionn): for ion in cls.ionlist.itervalues(): if ion[0] == ionn: return ion return cls.VOID @classmethod def from_ionn(cls, ionn): ion = cls.ionn2bfzae(ionn) if ion != cls.VOID: return KepIon(B = ion[0], F = ion[1], Z = ion[2], A = ion[3], E = ion[4]) return KepIon(cls.VOID) @classmethod def nuclei(cls): # should be replaced by a set return np.array([KepIon(n) for n,t in cls.ionlist.iteritems() if t[3] > 0]) def Ion(self): return Ion(A = self.A, Z = self.Z) ufunc_Ion = np.frompyfunc(lambda x: x.Ion(), 1, 1) @classmethod def from_ion(cls, ion): if not isinstance(ion, Ion): raise TypeError() return KepIon(ion.name()) def _name(self, upcase = True): t = self.tuple() for k,v in self.ionlist.items(): if v == t: return k return self.VOID_STRING def ionn(self): return self.B class DecIon(Ion): """ Class for decays. Will store A, Z, E in ion class portion Will set flag F_REACTION Will store extra info in "B" field, for now just the 12 recognized decay modes Provides list of ions, a dictionary with numbers - number of ions/particles that come out. The general reaction class should have these be negative for things that go in and positive for things that come out whereas the A, Z, E fields will store IN - OUT net values. LIMITATION: Currently this class only handles a pre-defined set of decays. PLAN: Cxtension for sf - complicated decays - maybe a separate derived class? TODO: For general class provide a generalized *index*. probbaly multiply the FZAE indices (plus BMUL/2 to take care of negative values in Z) with according numbers of BMUL? Add in flag for IN/OUT? It does not have to be nice, just unique... """ # list of decay names accepted reaction_list = { '' : 0, 's' : 0, 'g' : 1, 'b-' : 2, 'b+' : 3, 'ec' : 4, 'n' : 5, 'p' : 6, 'a' : 7, '2a' : 7, 'a2' : 8, '2p' : 8, 'p2' : 9, 'n2' : 9, '2n' : 10, 'np' : 10, 'pn' : 11, 'bn' : 12, 'b2n' : 13 } # OK, we want names to be unique on output reaction_names = { 0 : 's' , 1 : 'g' , 2 : 'b-' , 3 : 'b+' , 4 : 'ec' , 5 : 'n' , 6 : 'p' , 7 : 'a' , 8 : '2a' , 9 : '2p' , 10 : '2n' , 11 : 'pn' , 12 : 'bn' , 13 : 'b2n' } # list of particles that come out... particle_list = { 0 : {}, 1 : {Ion('g') : +1}, 2 : {Ion('e-') : 1}, 3 : {Ion('e+') : 1}, 4 : {Ion('e-') : -1}, 5 : {Ion('nt1') : 1}, 6 : {Ion('h1') : 1}, 7 : {Ion('he4') : 1}, 8 : {Ion('he4') : 2}, 9 : {Ion('h1') : 2}, 10 : {Ion('nt1') : 2}, 11 : {Ion('nt1') : 1, Ion('h1') : 1}, 12 : {Ion('b-') : 1, Ion('nt1') : 1}, 13 : {Ion('b-') : 1, Ion('nt1') : 2}} def _add(self, x, sign1 = +1, sign2 = +1): if not isinstance(x,Ion): x = Ion(x) if self.A == 0 and self.Z == 0: return Ion(Z = sign2 * x.Z, A = sign2 * x.A, E = max(sign2 * x.E + sign1 * self.E, 0) ) else: return Ion(Z = sign1 * self.Z + sign2 * x.Z, A = sign1 * self.A + sign2 * x.A) def __add__(self, x): self._add(x) __radd__ = __add__ def __sub__(self, x): self._add(x, +1, -1) def __rsub__(self, x): self._add(x, -1, +1) def __mul__(self, x): raise NotImplementedError("Can't multiply decays yet.") # later we could return a general reaction here def __init__(self, s): """ Set up decay reaction. Currently only a string is allowed for initialization. The main purpose is the interface to the decay.dat file and the Decay class setup in maxradio.py. """ if isinstance(s,str): self.B, self.F, self.Z, self.A, self.E = self.ion2bfzae(s) elif isinstance(s,type(self)): self.B, self.F, self.Z, self.A, self.E = s.tuple() else: raise AttributeError('wrong type') assert 0 <= self.B < len(self.reaction_names), "Unknown decay/reaction." self._particles = self.particle_list[self.B] self.update_idx() # for deepcopy: def __getstate__(self): return self.B def __setstate__(self, x): self.__init__(self.reaction_names[x]) @classmethod def ion2bfzae(cls,s=''): s = s.strip() if s.startswith('(') and s.endswith(')'): s=s[1:-1].strip() b = cls.reaction_list.get(s,-1) z, a, e = cls.particles2zae(cls.particle_list.get(b,{})) return b, cls.F_REACTION, z, a, e @staticmethod def particles2zae(particles): z, a, e = 0, 0, 0 for i,k in particles.items(): a -= i.A * k z -= i.Z * k e -= i.E * k return z, a, e def _name(self,low=None): """ Return pre-defined name from list. """ return self.reaction_names.get(self.B,self.VOID_STRING) def particles(self): """ Returns all paricles in the reaction/decay. """ return self._particles def hadrons(self): """ Returns just the hardrons in the reaction/decay. This is useful for networks like decay where photons and leptons usually are not accounted for. """ h = dict() for i,j in self._particles.items(): if i.A > 0: h[i] = j return h def leptons(self): """ Returns just the leptons in the reaction/decay. """ h = dict() for i,j in self._particles.items(): if i.F & self.F_LEPTON != 0: h[i] = j return h def photons(self): """ Returns just the photons (if any) in the reaction/decay. """ h = dict() for i,j in self._particles.items(): if (i.E != 0 and i.A == 0 and i.Z == 0 and i.F == 0): h[i] = j return h def isstable(self): """ Return if 'reaction'/'decay' is 'stable'. """ return self.B == 0 # not sure this will become sub- or super-class of DecIon class ReactIon(DecIon): """ Class for pre-defined set of reactions. I suppose we don't need separate 'In' and 'Out' list - just positive and negative values in the dictionary of particles. We could support 'In' and 'Out', though. We should add pre-defined strings for the common set of p,n,a,g,b+/b-. Will store A, Z, E in ion class portion Will set flag F_REACTION Will store extra info in "B" field, for now just the recognized reactions Provides list of ions, a dictionary with numbers - number of ions/particles that come out or come in. The general reaction class should have these be negative for things that go in and positive for things that come out whereas the A, Z, E fields will store IN - OUT net values. LIMITATION: Currently this class only handles a pre-defined set of reactions PLAN: Cxtension for sf - complicated reactions - maybe a separate derived class? TODO: For general class provide a generalized *index*. probbaly multiply the FZAE indices (plus BMUL/2 to take care of negative values in Z) with according numbers of BMUL? Add in flag for IN/OUT? It does not have to be nice, just unique... """ # list of reactions names accepted reaction_list = { '' : 0, 'b-' : 1, 'b+' : 2, 'ec' : 3, 'g,n' : 4, 'g,p' : 5, 'g,a' : 6, 'n,g' : 7, 'n,p' : 8, 'n,a' : 9, 'p,g' : 10, 'p,n' : 11, 'p,a' : 12, 'a,g' : 13, 'a,n' : 14, 'a,p' : 15, 'g' : 16, 'g*' : 17 } # OK, we want names to be unique on output reaction_names = { 0 : '' , 1 : 'b-' , 2 : 'b+' , 3 : 'ec' , 4 : 'g,n', 5 : 'g,p', 6 : 'g,a', 7 : 'n,g', 8 : 'n,p', 9 : 'n,a', 10 : 'p,g', 11 : 'p,n', 12 : 'p,a', 13 : 'a,g', 14 : 'a,n', 15 : 'a,p', 16 : 'g' , 17 : 'g*' } # list of particles that come out... particle_list = { 0 : {}, 1 : { Ion('e-') : 1}, 2 : { Ion('e+') : 1}, 3 : {Ion('e-') : -1 }, 4 : {Ion('g') : -1, Ion('nt1') : 1}, 5 : {Ion('g') : -1, Ion('h1') : 1}, 6 : {Ion('g') : -1, Ion('he4') : 1}, 7 : {Ion('nt1') : -1, Ion('g') : 1}, 8 : {Ion('nt1') : -1, Ion('h1') : 1}, 9 : {Ion('nt1') : -1, Ion('he4') : 1}, 10 : {Ion('h1') : -1, Ion('g') : 1}, 11 : {Ion('h1') : -1, Ion('nt1') : 1}, 12 : {Ion('h1') : -1, Ion('he4') : 1}, 13 : {Ion('he4') : -1, Ion('g') : 1}, 14 : {Ion('he4') : -1, Ion('nt1') : 1}, 15 : {Ion('he4') : -1, Ion('h1') : 1}, 16 : {Ion('g') : -1 }, 17 : { Ion('g') : 1}, } def _name(self,low=None): """ Return pre-defined name from list. """ return "({})".format(self.reaction_names.get(self.B,self.VOID_STRING)) def __mul__(self, x): raise NotImplementedError("Can't multiply tabulated reactions yet.") # later we could return a general reaction here class Reaction(ReactIon): """ This should become general reaction class. """ pass class IonSet(object): """ Provide sorted list of isotopes. TODO: efficient additions of arrays of ions. I suppose in principle we do want :set: properties but it also is supposed to be sorted by index at all times if _sort is set. Maybe there need be 2 kinds to allow interface with maxradio: maxradio also allows to have elements, etc., mixed. TODO - Maybe not KepIon and Ion mixed. """ def __init__(self, *args, **kwargs): self._ions = np.array((), dtype = np.object) self._sort = kwargs.get('sort', True) self._type = kwargs.get('type', None) assert self._type is None or not issubclass(self._type, Ion), 'type needs to be subclass of Ion' self.add(args) def add(self, ions = None): """ Add ion(s) to list. Ions need to be same "type" """ if len(ions) == 1: ions = ions[0] if (ions is not None): if np.isscalar(ions): self._add_one(ions) elif len(ions) == 1: while isinstance(ions, Iterable): ions = ions[0] self._add_one(ions) else: for ion in ions: self._add_one(ion, update = False) self._update() def _add_one(self, ion, update = True): """ add one ion and check compatibility """ if self._type == None: self._type = type(ion) if not issubclass(self._type, Ion): self._type = Ion if isinstance(ion, Ion) and type(ion) != self._type: raise TypeError("All ions need compatible type: {:s}, {:s}".format( str(self._type),str(type(ion)))) if not isinstance(ion, self._type): ion = self._type(ion) self._ions = np.append(self._ions, ion) if update: self._update() def _update(self): if self._sort: self._ions.sort() self._ions = np.unique(self._ions) def __str__(self): return "[{}]".format(','.join((str(i) for i in self._ions))) def __repr__(self): return "{}({})".format(self.__class__.__name__,str(self)) def __getitem__(self, index): return self._ions[index] def __len__(self): """ Return number of isotopes. """ return len(self._ions) def __iter__(self): i=0 while i < self.__len__(): yield self._ions[i] i+=1 class Abu(object): """ Ion plus single abundance. Provide some general functionallity. """ def __init__(self, ion, abu = 0): self.ion = ion self.abu = np.float64(abu) def __str__(self): return "{!s}:{:>12.5f}".format(self.ion,self.abu) def X(self): return self.abu def Y(self): return self.abu / self.ion.A class KepAbuDist(object): """ KepIon and distribution information, e.g., structure. """ # should be initialized with KepIon and abundace vector, # mass and radius would be good. pass class KepAbu(object): """ Singe KepIon with abundance. """ pass class KepAbuSet(object): """ Single Zone Kepler abudance. """ def __init__(self, iso = None, abu = None, comment = None, mixture = None, silent = False, **kwargs): """ Initialize KEPLER abundance set. Currently allow call with * list of isotopes only (abu set to 0) * list of isotopes and list of abundances * keyword iso=abu * dictionary {'iso':abu, ...} TODO - should always be sorted? TODO - should this be subclass of AbuSet? would need to define isotope class to use (KepIon) """ self.comment = stuple(comment) self.mixture = mixture if isinstance(iso, dict): self.iso, self.abu = self._ion_abu_from_dict(**iso) elif iso == None: assert abu is None, "Need isotope name" self.iso = np.array([],dtype=np.object) self.abu = np.array([],dtype=np.float64) else: self.iso = np.array([KepIon(i) for i in np.atleast_1d(iso)], dtype=np.object) if abu is not None: self.abu = np.array(np.atleast_1d(abu), dtype=np.float64) assert len(self.abu) == len(self.iso), "Need equal number of elements." else: self.abu = np.zeros_like(self.iso, dtype=np.float64) if len(kwargs) > 0: self._append(*self._ion_abu_from_dict(**kwargs)) @staticmethod def _ion_abu_from_dict(**kwargs): #todo: add sorting? return ( np.array([KepIon(i) for i in kwargs.keys()], dtype=np.object), np.array(kwargs.values(), dtype=np.float64)) def _append(self, iso, abu): # sorting? self.iso = np.append(self.iso, iso) self.abu = np.append(self.abu, abu) def __str__(self): return ("ion(" + ", ".join(['{:s}: {:f}'\ .format(iso.Name(),abu) for iso, abu in zip(self.iso,self.abu)]) + ")") __repr__ = __str__ def iteritems(self): """ Return pairs of (isotope, value). """ if self.iso.size == 0: raise StopIteration for iso, abu in zip(self.iso, self.abu): yield (iso, abu) def metallicity(self): """ Return 'metallicity' Z of composition. """ z = sum([abu for abu,iso in zip(self.abu,self.iso) if iso.Z >= 3]) return z def mu(self): """ Return mean molecular weight of composition. """ xmu = sum([abu/iso.A*(iso.Z+1) for abu,iso in zip(self.abu,self.iso)]) return 1./xmu def normalize(self, total = None): """ Normalize abundances to one. If sum == 0 just return. Note: Clone of AbuSet.normalize() """ abusum = self.abu.sum() if abusum == 0.: return self.abu /= abusum if total != None: self.abu *= total @classmethod def from_abu_to_approx(cls, abuset, comment = None): """ Map AbuSet to APPROX and return KepAbuSet """ from maxradio import MapBurn assert isinstance(abuset, AbuSet) approx_map = MapBurn(abuset.iso) # return self.approx_map(self.abu) return cls(abu = approx_map(abuset.abu), iso = approx_map.decions, comment = stuple(abuset.comment, comment), mixture = abuset.mixture) # maybe this routine is not needed ... ? def write_approx(self, outfile, comment = None, mixture = None, overwrite = False, silent = False): """ Write APPROX abundances generator to file. If outfile is file use this. If outfile is a filename open outfile for writing. """ self.setup_logger(silent = silent) # this code seems to repeat with OutFile(outfile) as f: # # TODO - test utils.OutFile context manager # if not isinstance(outfile, file): # filename = os.path.expanduser(os.path.expandvars(outfile)) # assert overwrite or not os.path.exists(filename) # f = open(filename,'w') # else: # f = outfile if not mixture: mixture = self.mixture assert mixture is not None for c in stuple(comment, self.comment): f.write('c {:s}'.format(c)) for iso, abu in zip(self.iso, self,abu): f.write('m {:s} {:25.17E} {:s}'.format( mixture, abu, iso.name())) # if not isinstance(outfile, file) # f.close() self.close_logger(timing='BURN generator written to "{:s}" in'.format(f.name)) class KepAbuData(object): """ Kepler abudance set. Helpful to convert ppn data with different networks as a function of zone into usable data. """ def __init__(self, ppn, netnum, ionn, wind = None, xm = None, zm = None): """ Initialize woth KEPLER network information. """ self.ppn = ppn self.netnum = netnum self.ionn = ionn self.wind = wind self.xm = xm self.zm = zm self.ions = np.ndarray(ionn.shape, dtype = np.object) nflat = self.ionn.reshape(-1) iflat = self.ions.reshape(-1) for i in xrange(len(nflat)): iflat[i] = KepIon.from_ionn(nflat[i]) # @cachedmethod def abu(self, ion, molfrac = False, missing = np.nan): """ Return ion or kepion abundace, depending on type. Sort of the cheap version of maxradio for ppn. """ if isinstance(ion, KepIon): yfac = max(1, ion.A) if molfrac: yfac = 1 value = np.ndarray(self.ppn.shape[1], dtype = np.float64) value.fill(missing) netlist = np.argwhere(ion.ionn() == self.ionn) for inet in xrange(netlist.shape[0]): zonelist = np.argwhere((netlist[inet,1] + 1) == self.netnum[1:-1]) if len(zonelist) > 0: value[zonelist + 1] = self.ppn[netlist[inet,0], zonelist + 1] * yfac # WIND data is in g not per gram # if self.wind is not None: # value[-1] = self.wind[ion.ionn-1] return value else: kepion = KepIon(ion) if kepion == VOID: value = np.ndarray(self.ppn.shape[1], dtype = np.float64) value.fill(missing) return value if kepion.B == 2: # pn1 seems to be defined in all three netwoks, but h1 not value = self.abu(kepion, molfrac = molfrac, missing = 0.) value += self.abu(KepIon.from_ionn(3), molfrac = molfrac, missing = missing) return value return self.abu(kepion, molfrac = molfrac, missing = missing) @cachedmethod def iron(self, molfrac = False, missing = np.nan): """ Return ion mass fraction (A > 46). """ value = np.zeros(self.ppn.shape[1], dtype = np.float64) for netnum in xrange(self.ionn.shape[1]): ion_index = np.where(Ion.ufunc_A(self.ions[:,netnum]) > 46)[0] zone_index = np.where(self.netnum[1:-1] == (netnum + 1))[0] + 1 if molfrac: A = np.ones(len(ion_index)) else: A = np.array([ion.A for ion in self.ions[ion_index,netnum]]) if len(zone_index) > 0: value[zone_index] = np.tensordot( self.ppn[np.ix_(ion_index,zone_index)], A, axes=(0,0)) # add wind values! Only reasoable for mass. value[[0,-1]] = missing return value @cachedmethod def metallicity(self, molfrac = False, missing = np.nan): """ Return mass fraction sum of A >= 5. """ value = np.zeros(self.ppn.shape[1], dtype = np.float64) for netnum in xrange(self.ionn.shape[1]): ion_index = np.where(Ion.ufunc_Z(self.ions[:,netnum]) > 2)[0] zone_index = np.where(self.netnum[1:-1] == (netnum + 1))[0] + 1 if molfrac: A = np.ones(len(ion_index)) else: A = np.array([ion.A for ion in self.ions[ion_index,netnum]]) if len(zone_index) > 0: value[zone_index] = np.tensordot( self.ppn[np.ix_(ion_index,zone_index)], A, axes=(0,0)) value[[0,-1]] = missing return value class AbuDist(object): """ Ion and distribution information, e.g., structure. """ # should be initialized with Ion and abundace vector, # mass and radius would be good. pass class AbuData(object): """ n-D abu data plus extra info. Probably base class (if not identical) to starfit data. """ def __init__(self, ions=None, fields=None, formats=None, parameters=None): super(self.__class__,self).__init__(self) #py3: super().__init__(self) class AbuSet(Logged): """ Prototype Class for single abundace set like solar abundance Should there be a subclass or a field for normalized/partial abundances? - Normalization is difficult if addition occurs piecewise ... How about support for different formats, like [] ... ? *** TODO *** - derive using common functions with IonSet - replace internal variables with "_"-names, most prominently: iso --> _ions abu --> _abu """ def __init__(self, iso = None, abu = None, comment = None, mixture = None, dat_file = None, bg_file = None, silent = False, sorted = True, **kwargs): """ Initialize abundance set. Currently allow call with * list of isotopes only (abu set to 0) * list of isotopes and list of abundances * keyword iso=abu * dictionary {'iso':abu, ...} Initialize from data file: * dat_file (KEPLER abundance data file, like solabu.dat) * bg_file (KEPLER BURN generator data file) TODO: sorting implementation """ self.comment = stuple(comment) self.mixture = mixture self.is_sorted = sorted if dat_file is not None: self._from_dat(dat_file, silent = silent) return if bg_file is not None: self._from_bg(bg_file, silent = silent) return if issubclass(dict, type(iso)): self.iso, self.abu = self._ion_abu_from_dict(**iso) elif iso == None: assert abu is None, "Need isotope name" self.iso = np.array([],dtype=np.object) self.abu = np.array([],dtype=np.float64) else: self.iso = np.array([Ion(i) for i in np.atleast_1d(iso)], dtype=np.object) if abu is not None: self.abu = np.array(np.atleast_1d(abu), dtype=np.float64) assert len(self.abu) == len(self.iso), "Need equal number of elements." else: self.abu = np.zeros_like(self.iso, dtype=np.float64) if len(kwargs) > 0: self._append(*self._ion_abu_from_dict(**kwargs)) if self.is_sorted: self.sort() @staticmethod def _ion_abu_from_dict(**kwargs): #todo: add sorting? #todo: check uniqueness return ( np.array([Ion(i) for i in kwargs.keys()], dtype=np.object), np.array(kwargs.values(), dtype=np.float64)) def _append(self, iso, abu): #todo: add sorting? #todo: check uniqueness self.iso = np.append(self.iso, iso) self.abu = np.append(self.abu, abu) def __str__(self): return ("abu(" + ", ".join(['{:s}: {:8G}'\ .format(iso.Name(),abu) for iso, abu in zip(self.iso,self.abu)]) + ")") __repr__ = __str__ def __add__(self, x): """ Note: duplicates are added NOTE - does return AbuSet class! """ assert isinstance(x, type(self)) # if empty we are done if x.iso.size == 0: return copy.deepcopy(self) # check for uniqueness in x assert np.unique(x.iso).size == x.iso.size, "Operand contains duplicate elements." abu_new = self.abu iso_new = self.iso size = self.iso.size # distinuish sorted and non-sorted cases if self.is_sorted: rev = np.arange(size) iso = iso_new else: i = self.iso.argsort() iso = self.iso[i] rev = i.argsort() add_list = list() for j,(i,a) in enumerate(zip(x.iso, x.abu)): idx = iso.searchsorted(i) if idx < size: if iso[idx] == i: abu_new[rev[idx]] += a else: add_list.append(j) else: add_list.append(j) if len(add_list) > 0: iso_new = np.append(iso_new, x.iso[add_list]) abu_new = np.append(abu_new, x.abu[add_list]) # here we can only return things with the same constructor # a = self.__class__(iso_new, abu_new, sorted = self.is_sorted) a = AbuSet(iso_new, abu_new, sorted = self.is_sorted or x.is_sorted) return a def __mul__(self, factor): # the following should have worked; there seems a bug in python 2.7.x # x = copy.deepcopy(self) # x.abu *= factor # <<< Instead we do x = copy.copy(self) x.iso = x.iso.copy() x.abu = x.abu * factor # >>> return x __rmul__ = __mul__ def normalize(self, total = None): """ Normalize abundances to one. If sum == 0 just return. """ abusum = self.abu.sum() if abusum == 0.: return self.abu /= abusum if total != None: self.abu *= total def normalized(self, total = None): """ Return normalized copy of abu. """ x = copy.copy(self) x.normalize(total = total) return x def _from_bg(self, file_name, silent = False): """ Generate abundace set from BURN gen file. TODO - return tuple with mutiple mixtures if file has several Actually, this should be a module function, not part of the constructor. TODO - add option to show comment """ self.setup_logger(silent = silent) self.comment = stuple(self.comment) xre = re.compile('[-+a-zA-Z0-9.]+') self.iso = np.array([],dtype=np.object) self.abu = np.array([],dtype=np.float64) with open(file_name,'r') as f: self.logger_file_info(f) for line in f: if line.startswith('c '): self.comment += (line[2:].rstrip(),) elif line.startswith('m '): xdata = xre.findall(line) xnum = len(xdata) if xnum < 2: continue mixture = xdata[1] if self.mixture is None: self.mixture = mixture self.logger.info('Loading mixture "{:s}".'.format(self.mixture)) if self.mixture == mixture: if xnum < 3: continue assert xnum % 2 == 0 xion = xdata[3::2] xabu = xdata[2::2] for i,a in zip(xion,xabu): self._append(Ion(i),np.double(a.replace('D','E').replace('d','e'))) self.close_logger(timing='Data loaded in') def _from_dat(self, file_name, silent = False): """ Load abundace set from "dat" file. TODO - add option to show comment """ self.setup_logger(silent = silent) self.comment = stuple(self.comment) xre = re.compile('[-+a-zA-Z0-9.]+') self.iso = np.array([],dtype=np.object) self.abu = np.array([],dtype=np.float64) with open(file_name,'r') as f: self.logger_file_info(f) self.comment += ('', 'Generated from file "{:s}".'.format(file_name), 'Original file comments follow:', '') for line in f: if not line.startswith((';','#')): xdata = xre.findall(line) xnum = len(xdata) if xnum == 0: continue if xnum == 2: xion,xabu = tuple(xdata) else: print(line) raise IOError('bad format') self._append(Ion(xion),np.double(xabu)) else: self.comment += (line[2:].rstrip(),) message = "{:3d} isotopes loaded in".format(len(self.iso)) self.close_logger(timing = message) def _new_order(self): # reset things that depend on order # this should be called every time things are added or removed. # try: # del self.approx_map # except: # pass pass def ionset(self): return IonSet(self.ions) def sort(self): """ Sort ions. """ if self.iso.size == 0: return sort = self.iso.argsort() self.iso = self.iso[sort] self.abu = self.abu[sort] self._new_order() # def map_approx(self, approx_comment = []): # """ # Map to APPROX and return KepAbuSet # """ # import maxradio # if not hasattr(self,'approx_map'): # self.approx_map = maxradio.MapBurn(self.iso) # # return self.approx_map(self.abu) # return KepAbuSet(abu = self.approx_map(self.abu), # iso = self.approx_map.decions, # comment = stuple(self.comment, approx_comment), # mixture = self.mixture) def write_bg(self, outfile, net = 1, mixture = None, zmax = 83, overwrite = False, silent = False): """ Write out BURN generator. If outfile is file use this. If outfile is a filename open outfile for writing. We need to assert gaps around be8, b9 (though it is not clear why KEPLER needs these gaps) We do need to assert, however, because ADAPNET preserves gaps for this reason, not to introduce gaps in new network. Which also makes the work of ADAPNET easier. Eventually KEPLER should check for these gaps and issue errors if present. Maybe it already does? """ self.setup_logger(silent = silent) # minnet from adapnet general purpose network # we could load this, but for now, we just hard-code. # note the gaps for Be8 and B9 that appear to be requires # by KEPLER - though it is not clear why # maybe because these istopes are not in bdat # TODO - ask RDH to add to future versions of bdat. minnet = [[[ 1, 1]], [[ 1, 3]], [[ 3, 4]], [[ 6, 7]], [[ 7, 7], [ 9, 9]], [[ 8, 8], [10, 11]], [[11,13]], [[13,15]], [[14,18]], [[17,19]], [[19,22]], [[21,23]], [[23,26]], [[25,27]], [[27,30]], [[30,31]], [[31,36]]] version = 10000 default_mixture = 'x' card_cmt = 'c' card_mix = 'm' card_grd = 'gg' card_ntw = 'netw' zmax = min(zmax+1, len(elements)) netw = np.ndarray((zmax,2,2), dtype=np.int64) netw_count = np.zeros(zmax, dtype=np.int64) for iz,z in enumerate(minnet): for ia,a in enumerate(z): netw[iz,ia,:] = a netw_count[iz] = len(z) niso = 0 iso = np.zeros_like(self.iso) abu = np.zeros_like(self.abu) for i,a in zip(self.iso, self.abu): Z = i.Z if Z < zmax: A = i.A n = netw_count[Z] if n == 0: netw_count[Z] = 1 netw[Z,0,:] = A else: netw[Z,0,0] = min(A, netw[Z,0,0]) netw[Z,n-1,1] = max(A, netw[Z,n-1,1]) add = True if n == 2: if (A > netw[Z,0,1]) and (A < netw[Z,1,0]): self.logger.error('{:s} inside required gap.'.format(i.name())) add = False if add: iso[niso] = i abu[niso] = a niso += 1 iso = iso[:niso] abu = abu[:niso] if mixture is None: try: mixture = self.mixture except: pass if mixture is None: mixture = default_mixture if not isinstance(outfile, file): filename = os.path.expanduser(os.path.expandvars(outfile)) assert overwrite or not os.path.exists(filename),'file exisits' f = open(filename,'w') else: f = outfile f.write(card_cmt + ' COMPUTER-GENERATED BURN GENERATOR FILE\n') f.write(card_cmt + ' VERSION {:s}'.format(version2human(version))+'\n') f.write(card_cmt + ' ' + time.asctime(time.gmtime())+' UTC\n') f.write(card_cmt + '\n') for c in self.comment: f.write("{:s} {:s}\n".format(card_cmt,c)) f.write(card_cmt + '\n') f.write(card_cmt + ' define network (Z_max = {:d})\n'.format(zmax-1)) for i in xrange(zmax): nc = netw_count[i] if nc > 0: c = " ".join(["{:3d} {:3d}".format(*netw[i,j,:]) for j in xrange(nc)]) f.write("{:s} {:d} {:2s} {:s}\n".format(card_ntw,net,elements[i],c)) f.write(card_cmt + '\n') f.write(card_cmt + ' define composition '+mixture+'\n') for i,a in zip(iso, abu): f.write("{:s} {:s} {:14.8E} {:s}\n".format(card_mix,mixture,a,i.name())) f.write(card_cmt + '\n') f.write(card_cmt + ' specify grid composition (homogeneous star)\n') f.write(card_cmt + ' NOTE: call only after all "g" cards in main generator\n') f.write("{:s} {:d} {:s}\n".format(card_grd,net,mixture)) if not isinstance(outfile, file): f.close() self.close_logger(timing='BURN generator written to "{:s}" in'.format(f.name)) def write_dat(self, outfile, overwrite = False, silent = False): """ Write out dat file. If file is file use this. If file is file name open file for write. """ self.setup_logger(silent = silent) version = 10000 card_cmt = ';' if not isinstance(outfile, file): filename = os.path.expanduser(os.path.expandvars(outfile)) assert overwrite or not os.path.exists(filename),' file exisits: '+filename f = open(filename,'w') else: f = outfile f.write(card_cmt + ' COMPUTER-GENERATED ABUNDANCE DATA FILE\n') f.write(card_cmt + ' VERSION {:s}'.format(version2human(version))+'\n') f.write(card_cmt + ' ' + time.asctime(time.gmtime())+' UTC\n') f.write(card_cmt + '\n') for c in self.comment: f.write("{:s} {:s}\n".format(card_cmt,c)) f.write(card_cmt + '\n') f.write(card_cmt + '----------------------------------------------------------------------\n') for i,a in zip(self.iso, self.abu): f.write("{:6s} {:13.7E}\n".format(i.name(),a)) if not isinstance(outfile, file): f.close() self.close_logger(timing = '"{:s}" written in'.format(f.name)) # some property routines def A(self): # return np.array([x.A() for x in self.iso],dtype=np.float64) return np.array([x.A for x in self.iso]) def Z(self): # return np.array([x.Z() for x in self.iso],dtype=np.float64) return np.array([x.Z for x in self.iso]) def N(self): # return np.array([x.Z() for x in self.iso],dtype=np.float64) return np.array([x.N for x in self.iso]) def E(self): # return np.array([x.Z() for x in self.iso],dtype=np.float64) return np.array([x.E for x in self.iso]) def element(self): return np.array([Ion(Z=x.Z) for x in self.iso]) def element_name(self): return np.array([x.element_name() for x in self.iso]) def element_symbol(self, upcase = True): return np.array([x.element_symbol(upcase = upcase) for x in self.iso]) def idx(self): return np.array([x.idx for x in self.iso]) def XYZ(self): """ Return 'astronomical' X, Y, Z of composition. """ x = sum(self.X()[self.Z() == 1]) y = sum(self.X()[self.Z() == 2]) z = sum(self.X()[self.Z() >= 3]) return x,y,z def metallicity(self): """ Return 'metallicity' Z of composition. """ z = sum(self.X()[self.Z() >= 3]) return z def Ye(self): """ Return electron to baryon ratio. """ Ye = np.sum(self.Z()*self.Y()) return Ye def mue(self): """ Return mean molecular weight per electron of composition. """ return 1./self.Ye() def eta(self): """ Return neutron excess of mixture. """ return 1.-2.*self.Ye() def mu(self): """ Return mean molecular weight of composition. """ xmu = np.sum((self.Z()+1)*self.Y()) return 1./xmu # here to add some general routines that return data for any Ion type # (isotope, isotone, isobar, element) def _X(self, selection = None, data = 'X'): """ Return data for selected (isotope, isotone, isobar, element). If nothing is specified, all isotopes will be returned. Default is mass fraction 'X'. """ # setup uis_func = [Ion.ufunc_is_isotope, Ion.ufunc_is_element, Ion.ufunc_is_isobar, Ion.ufunc_is_isotone] # we would only have to switch out this one definition to get them all... if data == 'X': val_func = [self.isotopes_X, self.elements_X, self.isobars_X, self.isotones_X] elif data == 'Y': val_func = [self.isotopes_Y, self.elements_Y, self.isobars_Y, self.isotones_Y] elif data == 'log_eps': val_func = [self.isotopes_log_eps, self.elements_log_eps, self.isobars_log_eps, self.isotones_log_eps] else: raise ValueError("Invalid Data request: '"+data+"'.") # selection setup if selection is None: return val_func[0]() if np.shape(selection) == (): selection = np.array([selection]) if issubclass(type(selection[0]),str): selection = np.array([Ion(ion) for ion in selection]) selections = np.ndarray([len(uis_func),len(selection)]) # pure cases # here we test the most likely cases first ... for i,(uis,val) in enumerate(zip(uis_func,val_func)): selections[i,:] = uis(selection) if np.all(selections[i,:]): return val(selection) # now mixed cases x = np.zeros(len(selection), dtype=np.float64) for i,val in enumerate(val_func): if np.any(selections[i,:]): index, = np.nonzero(selections[i,:]) x[index] = val(selection[index]) return x def _Y(self): return self.abu / self.A() def X(self, selection = None): """ Return mass fraction for selected (isotope, isotone, isobar, element). If nothing is specified, all isotopes will be returned. """ return self._X(selection, data = 'X') def Y(self, selection = None): """ Return mole/g for selected (isotope, isotone, isobar, element). If nothing is specified, all isotopes will be returned. """ return self._X(selection, data = 'Y') def log_eps(self, selection = None): """ Return log(eps) for selected (isotope, isotone, isobar, element). If nothing is specified, all isotopes will be returned. """ return self._X(selection, data = 'log_eps') # add geophysical relative to si28? # add functions to compute [] , \delta, # isotope routines def _isotopes_selection(self, x, selection = None): if selection != None: # we should really use ion sets if np.shape(selection) == (): selection = np.array([selection]) if issubclass(type(selection[0]),str): selection = np.array([Ion(ion) for ion in selection]) if issubclass(type(selection[0]),Ion): selection = np.array([ion.idx for ion in selection]) else: raise AttributeError("Selection argument invalid.") x = x[np.in1d(self.idx(), selection)] if len(x) < len(selection): raise KeyError("Selection contained invalid entries.") return x def isotopes_X(self, selection = None): """ Return isotope mass fractions. """ return self._isotopes_selection(self.abu, selection = selection) def isotopes_Y(self, selection = None): """ Return isotope number fractions (mol per gram). """ return self._isotopes_selection(self.abu / self.A(), selection = selection) def isotopes_log_eps(self, selection = None): """ Return isotope log(eps). log10(# of atoms relative to H) + 12 """ h = self.Z() == 1 A = self.A() x = sum(self.abu[h] / A[h]) y = self._isotopes_selection(self.abu / A, selection = selection) return np.log10(y / x) + 12. # element routines def elements(self): """ Return list of all elements in abundance set. """ return np.array([Ion(Z = x) for x in self.elements_Z()], dtype=np.object) def elements_Z(self): """ Return charge number of all elements. """ return np.unique(self.Z()) def _elements_selection(self, x, selection = None): if selection != None: # we should really use ion sets if np.shape(selection) == (): selection = np.array([selection]) if issubclass(type(selection[0]),str): selection = np.array([Ion(ion).Z for ion in selection]) if issubclass(type(selection[0]),Ion): selection = np.array([ion.Z for ion in selection]) x = x[np.in1d(self.elements_Z(),selection)] if len(x) < len(selection): raise KeyError("Selection contained invalid entries.") return x def elements_X(self, selection = None): """ Return elemental mass fractions. """ return self._elements_selection(project(self.abu,self.Z()), selection = selection) def elements_Y(self, selection = None): """ Return elemental number fractions (mol per gram). """ return self._elements_selection(project(self._Y(),self.Z()), selection = selection) def elements_log_eps(self, selection = None): """ Return elemental log(eps). log10(# of atoms relative to H) + 12 """ h = self.Z() == 1 A = self.A() x = sum(self.abu[h] / A[h]) y = self._elements_selection(project(self._Y(),self.Z()), selection = selection) return np.log10(y / x) + 12. def elements_name(self): """ Return element name of each isotope. """ return np.array([Elements[x] for x in self.elements_Z()]) def elements_set(self, selection = None): """ Return set of projected elements. """ abu = self._elements_selection(project(self.abu,self.Z()), selection = selection) iso = self._elements_selection(np.array([ Ion(Z = Z) for Z in self.elements_Z()], dtype = object), selection = selection) return AbuSet(abu = abu, iso = iso, comment = self.comment) # we may want to add something signifying elements ... # isobar routines def isobars(self): """ Return list of all elements in abundance set. """ return np.array([Ion(A = x) for x in self.isobars_A()], dtype=np.object) def isobars_A(self): """ Return mass number of all isobars. """ return np.unique(self.A()) def _isobars_selection(self, x, selection = None): if selection != None: # we should really use ion sets if np.shape(selection) == (): selection = np.array([selection]) if issubclass(type(selection[0]),Ion): selection = np.array([ion.A for ion in selection]) x = x[np.in1d(self.isobars_A(),selection)] if len(x) < len(selection): raise KeyError("Selection contained invalid entries.") return x def isobars_X(self, selection=None): """ Return isobar mass fractions. """ return self._isobars_selection(project(self.abu,self.A()), selection = selection) def isobars_Y(self, selection=None): """ Return isobar number fractions (mol per gram). """ return self._isobars_selection(project(self._Y(),self.A()), selection = selection) def isobars_log_eps(self, selection = None): """ Return isobar log(eps). log10(# of atoms relative to H) + 12 """ h = self.Z() == 1 A = self.A() x = sum(self.abu[h] / A[h]) y = self._isobars_selection(project(self._Y(),A), selection = selection) return np.log10(y / x) + 12. def isobars_set(self, selection = None): """ Return set of projected isobars. """ abu = self._isobars_selection(project(self.abu,self.A()), selection = selection) iso = self._isobars_selection(np.array([ Ion(A = A) for A in self.isobars_A()], dtype = object), selection = selection) return AbuSet(abu = abu, iso = iso, comment = self.comment) # we may want to add something signifying isobars ... # isotone routines def isotones(self): """ Return list of all isotones in abundance set. """ return np.array([Ion(N = x) for x in self.isotones_N()], dtype=np.object) def isotones_N(self): """ Return neutron number of all isotones. """ return np.unique(self.N()) def _isotones_selection(self, x, selection = None): if selection != None: # we should really use ion sets if np.shape(selection) == (): selection = np.array([selection]) if issubclass(type(selection[0]),Ion): selection = np.array([ion.N for ion in selection]) x = x[np.in1d(self.isotones_N(),selection)] if len(x) < len(selection): raise KeyError("Selection contained invalid entries.") return x def isotones_X(self, selection=None): """ Return isotone mass fractions. """ return self._isotones_selection(project(self.abu,self.N()), selection = selection) def isotones_Y(self, selection=None): """ Return isotone number fractions (mol per gram). """ return self._isotones_selection(project(self._Y(),self.N()), selection = selection) def isotones_log_eps(self, selection = None): """ Return isotone log(eps). log10(# of atoms relative to H) + 12 """ h = self.Z() == 1 A = self.A() x = sum(self.abu[h] / A[h]) y = self._isotones_selection(project(self._Y(),self.N()), selection = selection) return np.log10(y / x) + 12. def isotones_set(self, selection = None): """ Return set of projected isotones. """ abu = self._isotones_selection(project(self.abu, self.N()), selection = selection) iso = self._isotones_selection(np.array([ Ion(N = N) for N in self.isotones_N()], dtype = object), selection = selection) return AbuSet(abu = abu, iso = iso, comment = self.comment) # we may want to add something signifying isotones ... # general access interfaces def __getitem__(self, index): try: # this does not work for numpy # TDOD add "where" or "searchsorted" resolution if isinstance(index, str): index = Ion(index) if isinstance(index, Ion): index = np.where(self.iso == index)[0][0] return self.abu[index] except: raise AttributeError('Isotope not found.') def __setitem__(self, index, item): # TODO add isotope if not in list? # maybe add parameter allow_new try: # this does not work for numpy if isinstance(index, str): index = Ion(index) if isinstance(index, Ion): index = np.where(self.iso == index) if not is_numlike(item): raise AttributeError('Abundance needs to be numlike.') self.abu[index] = item except: raise AttributeError('Isotope not found.') def __len__(self): """ Return number of isotopes. """ return len(self.iso) def __getattr__(self, attr): if self.__dict__.has_key('iso'): x = np.where(self.iso == attr) if len(x) > 0: x=x[0] if len(x) == 1: # return np.array(self.abu[x[0]],dtype=np.float64,ndmin=1) return self.abu[x[0]] # raise AttributeError('Isotope/Attribute not found.') # if hasattr(super(AbuSet, self), '__getattr__'): super(AbuSet, self).__getattr__(attr) #py3: super().__getattr__(attr) def __setattr__(self,attr,value): if not self.__dict__.has_key(attr): if self.__dict__.has_key('iso'): x=np.where(self.iso == attr)[0] if len(x) == 1: i = x[0] fac = 1 - value self.abu[i] = 0. self._normalize(fac) self.abu[i] = value return # it does suck that python 2.x requires class super(AbuSet, self).__setattr__(attr,value) #py3: super().__setattr__(attr,value) def __iter__(self): i=0 while i < self.__len__(): yield (self.iso[i],self.abu[i]) i+=1 def __call__(self, *args, **kwargs): """ TODO: 1) add new isotopes to abundance pattern 2) renormailzie 3) check value is in Range 0 1.: jj, = np.argwhere(scaled.iso == Ion('He4')) # make sure we have same set of istopes bbn = (sun * 0.) + bbn for j in np.argwhere(scaled.abu < sun.abu).flat: scaled.abu[jj] += scaled.abu[j] scaled.abu[j] = sun.abu[j] * np.exp( (scale - 1.)*(1. - bbn.abu[j]/sun.abu[j])) scaled.abu[jj] -= scaled.abu[j] scaled.normalize() self.iso = scaled.iso self.abu = scaled.abu self.is_sorted = scaled.is_sorted self.comment = ( "Version {:6s} - {:s}".format(self.version,time.asctime(time.gmtime())+' UTC'), "Scaled solar abundances: {:g} solar".format(scale), "Sun: {:s} - {:s}".format(solar,solar_file), "BBN: {:s} - {:s}".format(zero,zero_file), "X = {:8G}, Y = {:8G}, Z = {:8G}".format(*self.XYZ())) self.logger.info("X = {:8G}, Y = {:8G}, Z = {:8G}".format(*self.XYZ())) self.close_logger(timing='Abundance set generated in') class Mass(Logged): """ Object to hold ion masses. Not clear at this point where the informastion will come from. Currently load audi 2003 For the rest: for now just use A TODO: use N*m_n + Z*m_p - Q/c**2 (from bdat) TODO: Option to just use A """ default_data = 'Au03' def __init__(self, data = default_data, silent = False): """ TODO - implement data """ self.setup_logger(silent = silent) path = os.getenv('KEPLER_DATA') if not path: path=os.path.join(os.path.expanduser('~'),'kepler','local_data') self.logger.warning('using default path ' + path) file_name = os.path.join(path,'masses_audi_2003.dat') self.comment = () self.iso = np.array([],dtype=np.object) self.mass = np.array([],dtype=np.float64) xre = re.compile('[-+a-zA-Z0-9.]+') with open(file_name,'r') as f: self.logger_file_info(f) for line in f: if not line.startswith((';','#')): xdata = xre.findall(line) xnum = len(xdata) if xnum == 0: continue if xnum == 2: xion,xabu = tuple(xdata) else: print(line) raise IOError('bad format') self._append(Ion(xion),np.double(xabu)) else: self.comment += (line[2:],) message = "{:3d} masses loaded in".format(len(self.iso)) self.close_logger(timing = message) def _append(self, iso, abu): self.iso = np.append(self.iso, iso) self.mass = np.append(self.mass, abu) def __call__(self, ion): """ Return mass of ion in amu (or approximate) TODO: accept list """ try: i, = np.argwhere(self.iso == ion) return self.mass[i[0]] except: pass return np.double(Ion(ion).A) def __str__(self): return ("mass(" + ", ".join(['{:s}: {:f}'\ .format(iso.Name(),mass) for iso, mass in zip(self.iso,self.mass)]) + ")") __repr__ = __str__ class Asplund2009Data(AbuSet): """ Routine to load Asplund 2009 solar abundances Here is what I used to call: x = isotope.Asplund2009Data() x.write_dat('~/kepler/local_data/solas12.dat') x.write_bg('~/kepler/local_data/solas12g') """ def __init__(self, file_name = '~/Plots/solar/Asplund2009-isotopes_protosun.dat', comment = None, silent = False): """ Load abundace set from Aspund "dat" file. TODO - add option to show comment """ self.setup_logger(silent = silent) comment = stuple(comment) xre = re.compile('[-+a-zA-Z0-9.]+') iso = np.array([],dtype=np.object) abu = np.array([],dtype=np.float64) file_name = os.path.expanduser(file_name) with open(file_name,'r') as f: self.logger_file_info(f) comment += ('', 'Generated from file "{:s}".'.format(file_name), 'Original file comments follow:', '') for line in f: if not line.startswith((';','#')): xdata = xre.findall(line) xnum = len(xdata) if xnum == 0: continue if xnum == 5: xiso = Ion(A = int(xdata[2]), Z = int(xdata[1])) xabu = np.double(xdata[4]) else: print(line) raise IOError('bad format') iso = np.append(iso, xiso) abu = np.append(abu, xabu) else: comment += (line[2:].rstrip(),) m = Mass() # well, this could require tests... abu = np.array([a*m(i) for i,a in zip(iso,abu)]) abu = abu/abu.sum() super(self.__class__, self).__init__( iso = iso, abu = abu, comment = comment) #py3: super().__init__() message = "{:3d} isotopes loaded in".format(iso.size) self.close_logger(timing = message)