''' utils.py Utility class for holding extra methods from mesa.py, nuh5p.py ''' import numpy as np import scipy as sc import ascii_table as att from scipy import optimize import matplotlib.pyplot as pl import os class data_fitting(): ''' Wrapper for the scipy method optimize.leastsq Parameters ---------- func : string or function, optional If func is a string, then it must be one of the two strings 'linear' or 'powerlaw'. If func is not a string it must be a custom function. The default is "linear". coef : tuple or list A guess for the list of coeffiecients. For 'powerlaw' coef must have three enries. For 'linear' coef must have two enries. If you provide your own function, then provide as many coef entries as your function needs. Examples -------- Typically you have some data y(x) and you want to fit some function to this data. in order to test this create some sample data y(x) >>> x=arange(0,100) >>> y=(((random_sample(100)+0.5)*50.)-50)+x # this is the data to fit This class provides two simple fit functions: linear and power law. Feel free to add more, or provide your custom function as an argument. For example, in order to provide an exponential fit function, first define the function. >>> def ff(coef,x): >>> return coef[0]*sc.exp(coef[1]*x)+coef[2] and then call the data_fitting instance. >>> f=utils.data_fitting(ff,coef=(1,0,0)) Once you have initialized this class, the instance provides two methods; fit and plot to check the fit. >>> f.fit(x,y) >>> f.plot() The fit coefficients are stored in self.fcoef. (FH) ''' def __init__(self, func='linear', coef=(1, 1)): ''' Parameters ---------- func : string or function, optional If func is a string, then it must be one of the two strings 'linear' or 'powerlaw'. If func is not a string it must be a custom function. The default is "linear". coef : list, optional A guess for the list of coeffiecients. For 'powerlaw' coef must have three enries. For 'linear' coef must have two enries. If you provide your own function, then provide as many coef entries as your function needs. The default is (1, 1). ''' if func is 'linear': print "Information: 'linear' fit needs coef list with 2 entries" print " -> will use default: coef = "+str(coef) if len(coef) is not 2: print "Warning: you want a linear fit but you have not" print " provided a guess for coef with the" print " right length (2)." print " -> I will continue and assume coef=(1,1)" coef = (1,1) def ff(coef,x): return coef[0]*x + coef[1] self.func_name = func elif func is 'powerlaw': print "Information: 'powerlaw' fit needs coef list with 3 entries" print " -> will use default: coef = "+str(coef) if len(coef) is not 3: print "Warning: you want a power law fit but you have" print " not provided a guess for coef with the" print " right length (3)." print " -> I will continue and assume coef=(1,1,1)" coef = (1,1,1) def ff(coef,x): return coef[0]*x**coef[1] + coef[2] self.func_name = func else: print "Information: You provide a fit function yourself. I trust" print " you have provided a matching guess for the " print " coefficient list!" ff = func self.func_name = func.__name__ # we want to determine the coefficients that # fit the power law to the data # this is done by finding the minimum to a # residual function: # func(params) = ydata - f(xdata, params) # therefore we define a residual function def fres(coef,y,ff,x): return y-ff(coef,x) self.residual = fres self.coef = coef self.func = ff def fit(self, x, y, dcoef='none'): ''' performs the fit x, y : list Matching data arrays that define a numerical function y(x), this is the data to be fitted. dcoef : list or string You can provide a different guess for the coefficients, or provide the string 'none' to use the inital guess. The default is 'none'. Returns ------- ierr Values between 1 and 4 signal success. Notes ----- self.fcoef, contains the fitted coefficients. ''' self.x = x self.y = y if dcoef is not 'none': coef = dcoef else: coef = self.coef fcoef=optimize.leastsq(self.residual,coef,args=(y,self.func,x)) self.fcoef = fcoef[0].tolist() return fcoef[1] def plot(self, ifig=1, data_label='data', fit_label='fit', data_shape='o', fit_shape='-'): ''' plot the data and the fitted function. Parameters ---------- ifig : integer Figure window number. The default is 1. data_label : string Legend for data. The default is 'data'. fit_label : string Legend for fit. If fit_lable is 'fit', then substitute fit function type self.func_name. The default is 'fit'. data_shape : character Shape for data. The default is 'o'. fit_shape : character Shape for fit. The default is '-'. ''' if len(self.coef) is not len(self.fcoef): print "Warning: the fitted coefficient list is not same" print " length as guessed list - still I will try ..." pl.figure(ifig) pl.plot(self.x,self.y,data_shape,label=data_label) if fit_label is 'fit': fit_label=self.func_name pl.plot(self.x,self.func(self.fcoef,self.x),fit_shape,label=fit_label) pl.legend() class constants(): mass_sun=1.9891e+33 mass_sun_unit='g' one_year=31558149.984 avogadro=6.02214179e23 avogadro_unit='mol^-1' class Utils(): ''' This private class contains utilities that are used by methods, mostly in the ppn and mppnp classes. Users whould normally not use these methods directly. Things go here when it can be imagined that they may be used not in immediate conjunction with plotting. Otherwise they would go into the superclass data_plot. ''' #elements_names is marked for deletion (FH, Oct2011) and #should be replaced with self.stable_names from self._stable_names() elements_names = ['Neutron','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'] stable_el = [['Neutron','999'],['H',1, 2],['He', 3, 4],['Li', 6, 7],['Be', 9], ['B', 10, 11],['C', 12, 13],['N', 14, 15],['O', 16, 17, 18],['F', 19], ['Ne', 20, 21, 22],['Na', 23],['Mg', 24, 25, 26],['Al', 27],['Si', 28, 29, 30], ['P', 31],['S', 32, 33, 34, 36],['Cl', 35, 37],['Ar', 36, 38, 40],['K', 39, 40, 41], ['Ca', 40, 42, 43, 44, 46, 48],['Sc', 45],['Ti', 46, 47, 48, 49, 50],['V', 50, 51], ['Cr', 50, 52, 53, 54],['Mn', 55],['Fe', 54, 56, 57, 58],['Co', 59], ['Ni', 58, 60, 61, 62, 64],['Cu', 63, 65],['Zn', 64, 66, 67, 68, 70],['Ga', 69, 71], ['Ge', 70, 72, 73, 74, 76],['As', 75],['Se', 74, 76, 77, 78, 80, 82],['Br', 79, 81], ['Kr', 78, 80, 82, 83, 84, 86],['Rb', 85, 87],['Sr', 84, 86, 87, 88],['Y', 89], ['Zr', 90, 91, 92, 94, 96],['Nb', 93],['Mo', 92, 94, 95, 96, 97, 98, 100], ['Tc',999],['Ru', 96, 98, 99, 100, 101, 102, 104],['Rh', 103], ['Pd', 102, 104, 105, 106, 108, 110],['Ag', 107, 109], ['Cd', 106, 108, 110, 111, 112, 113, 114, 116],['In', 113, 115], ['Sn', 112, 114, 115, 116, 117, 118, 119, 120, 122, 124],['Sb', 121, 123], ['Te', 120, 122, 123, 124, 125, 126, 128, 130],['I', 127], ['Xe', 124, 126, 128, 129, 130, 131, 132, 134, 136],['Cs', 133], ['Ba', 130, 132, 134, 135, 136, 137, 138],['La', 138, 139],['Ce', 136, 138, 140, 142], ['Pr', 141],['Nd', 142, 143, 144, 145, 146, 148, 150],['Pm',999], ['Sm', 144, 147, 148, 149, 150, 152, 154],['Eu', 151, 153], ['Gd', 152, 154, 155, 156, 157, 158, 160],['Tb', 159], ['Dy', 156, 158, 160, 161, 162, 163, 164],['Ho', 165], ['Er', 162, 164, 166, 167, 168, 170],['Tm', 169],['Yb', 168, 170, 171, 172, 173, 174, 176], ['Lu', 175, 176],['Hf', 174, 176, 177, 178, 179, 180],['Ta', 180, 181], ['W', 180, 182, 183, 184, 186],['Re', 185, 187],['Os', 184, 186, 187, 188, 189, 190, 192], ['Ir', 191, 193],['Pt', 190, 192, 194, 195, 196, 198],['Au', 197], ['Hg', 196, 198, 199, 200, 201, 202, 204],['Tl', 203, 205],['Pb', 204, 206, 207, 208], ['Bi', 209],['Th', 232],['U',235,238],['Po',999]] def _stable_names(self): ''' This private method extracts the element names from stable_el. Note that stable_names is a misnomer as stable_el also contains unstable element names with a number 999 for the *stable* mass numbers. (?!??) ''' stable_names=[] for i in range(len(self.stable_el)): stable_names.append(self.stable_el[i][0]) self.stable_names=stable_names def _process_abundance_vector(self, a, z, isomers, yps): ''' This private method takes as input one vector definition and processes it, including sorting by charge number and mass number. It returns the processed input variables plus an element and isotope vector and a list of isomers. ''' tmp=[] isom=[] for i in xrange(len(a)): if z[i]!=0 and isomers[i]==1: #if its not 'NEUt and not an isomer' tmp.append([self.stable_names[int(z[i])]+'-'+str(int(a[i])),yps[i],z[i],a[i]]) elif isomers[i]!=1: #if it is an isomer if yps[i]==0: isom.append([self.stable_names[int(z[i])]+'-'+str(int(a[i]))+'-'+str(int(isomers[i]-1)),1e-99]) else: isom.append([self.stable_names[int(z[i])]+'-'+str(int(a[i]))+'-'+str(int(isomers[i]-1)),yps[i]]) tmp.sort(self.compar) tmp.sort(self.comparator) abunds=[] isotope_to_plot=[] z_iso_to_plot=[] a_iso_to_plot=[] el_iso_to_plot=[] for i in xrange(len(tmp)): isotope_to_plot.append(tmp[i][0]) abunds.append(tmp[i][1]) z_iso_to_plot.append(int(tmp[i][2])) a_iso_to_plot.append(int(tmp[i][3])) el_iso_to_plot.append(self.stable_names[int(tmp[i][2])]) return a_iso_to_plot,z_iso_to_plot,abunds,isotope_to_plot,el_iso_to_plot,isom def compar(self, x, y): ''' simple comparator method ''' indX=0 indY=0 a= int(x[0].split('-')[1]) b= int(y[0].split('-')[1]) if a>b: return 1 if a==b: return 0 if aindY: return 1 if indX==indY: return 0 if indX beta- decay # else beta > 0 => beta+ decay # now we need an index array on the scale of the abundance # distribution to be plotted that points to itself for stable species, # and to the stable element to which it decays in case of an unstable # species decay_index_pointer=np.zeros(len(isotope_to_plot), dtype='int')-1 idp_to_stables_in_isostoplot=[] for i in range(len(isotope_to_plot)): element_name=isotope_to_plot[i].split('-')[0] try: stable_a=stable_a_db[np.where(el_db==element_name)][0] # 4th column for that element in isotopedatabase.txt except IndexError: print "Can't find element "+element_name+" in isotopedatabase.txt" if a_iso_to_plot[i] <= 209 and stable_a <=209: # Bi209 is last stable element stable_mass_numbers=self.stable_el[self.stable_names.index(element_name)][1:] iso_db_index_range_el=np.where(el_db==element_name) beta_for_this_species=\ beta[iso_db_index_range_el][np.where(a_db[iso_db_index_range_el]==a_iso_to_plot[i])] if beta_for_this_species == 0: # if there are no stable species for an element (Tc,Pm) the cutoff specifies beta_for_this_species = -1 # the lowest mass beta- isotope if a_iso_to_plot[i] in stable_mass_numbers: # print isotope_to_plot[i]+" is stable" decay_index_pointer[i]=i idp_to_stables_in_isostoplot.append(i) elif a_iso_to_plot[i]==8: # Be8 -> He4 decay_index_pointer[i]=isotope_to_plot.index('He-4') else: # beta decay found_decay_target=False i_search=-1*beta_for_this_species while not found_decay_target: try: try_target_el=self.stable_names[charge_from_element_name[element_name]+i_search] except TypeError: print "Maybe information about species "+isotope_to_plot[i]+" is not available in isotopedatabase.txt" decay_index_pointer[i]=-1 break # print try_target_el try: stable_mass_numbers=self.stable_el[self.stable_names.index(try_target_el)][1:] except ValueError: print "Can not find decay target for "+isotope_to_plot[i] if a_iso_to_plot[i] in stable_mass_numbers: ind_range=np.where(np.array(el_iso_to_plot)==try_target_el)[0] if a_iso_to_plot[i] in np.array(a_iso_to_plot)[ind_range]: this_ind=\ ind_range[np.where(np.array(a_iso_to_plot)[ind_range]==a_iso_to_plot[i])[0]] # print isotope_to_plot[i]+" is unstable and decays to "+isotope_to_plot[this_ind] decay_index_pointer[i]=this_ind else: print "It seems unstable species "+isotope_to_plot[i]+" wants to decay to " \ +try_target_el+"-"+str(a_iso_to_plot[i])+", however this species is not in this run." \ +" This points to an inconsistency in the network build. Here we will ignore the abundance of " \ +isotope_to_plot[i]+'.' decay_index_pointer[i]=-1 found_decay_target=True else: i_search += -1*beta_for_this_species if self.debug: print "Decay rules:" for i in range(len(isotope_to_plot)): if decay_index_pointer[i]>= 0: print isotope_to_plot[i]+" -> "+isotope_to_plot[decay_index_pointer[i]] ind_tmp=idp_to_stables_in_isostoplot #ind_tmp=utils.strictly_monotonic(decay_index_pointer) # this would do the same, but the method above is more straight forward self.decay_idp=decay_index_pointer self.idp_to_stables_in_isostoplot=ind_tmp def is_stable(self,species): ''' This routine accepts input formatted like 'He-3' and checks with stable_el list if occurs in there. If it does, the routine returns True, otherwise False. Notes ----- this method is designed to work with an se instance from nugridse.py. In order to make it work with ppn.py some additional work is required. FH, April 20, 2013. ''' element_name_of_iso = species.split('-')[0] try: a_of_iso = int(species.split('-')[1]) except ValueError: # if the species name contains in addition to the # mass number some letters, e.g. for isomere, then # we assume it is unstable. This is not correct but # related to the fact that in nugridse.py we do not # identify species properly by the three numbers A, Z # and isomeric_state. We should do that!!!!!! a_of_iso = 999 idp_of_element_in_stable_names = self.stable_names.index(element_name_of_iso) if a_of_iso in self.stable_el[idp_of_element_in_stable_names][1:]: return True else: return False class iniabu(Utils): ''' This class in the utils package reads an abundance distribution file of the type iniab.dat. It then provides you with methods to change some abundances, modify, normalise and eventually write out the final distribution in a format that can be used as an initial abundance file for ppn. This class also contains a method to write initial abundance files for a MESA run, for a given MESA netowrk. ''' # clean variables that we will use in this class filename = '' def __init__(self,filename): ''' Init method will read file of type iniab.dat, as they are for example found in the frames/mppnp/USEPP directory. Notes ----- An instance of this class will have the following data arrays z charge number. a mass number. abu abundance. names name of species. habu a hash array of abundances, referenced by species name. hindex hash index returning index of species from name. Examples -------- E.g. if x is an instance then x.names[4] gives you the name of species 4, and x.habu['c 12'] gives you the abundance of C12, and x.hindex['c 12'] returns 4. Note, that you have to use the species names as they are provided in the iniabu.dat file. Example - generate modified input file ppn calculations: >>> import utils >>> p=utils.iniabu('iniab1.0E-02.ppn_asplund05') >>> sp={} >>> sp['h 1']=0.2 >>> sp['c 12']=0.5 >>> sp['o 16']=0.2 >>> p.set_and_normalize(sp) >>> p.write('p_ini.dat','header for this example') p.write_mesa allows you to write this NuGrid initial abundance file into a MESA readable initial abundance file. ''' f0=open(filename) sol=f0.readlines() f0.close # Now read in the whole file and create a hashed array: names=[] z=[] yps=np.zeros(len(sol)) mass_number=np.zeros(len(sol)) for i in range(len(sol)): z.append(int(sol[i][1:3])) names.extend([sol[i].split(" ")[0][4:]]) yps[i]=float(sol[i].split(" ")[1]) try: mass_number[i]=int(names[i][2:5]) except ValueError: print "WARNING:" print "This initial abundance file uses an element name that does" print "not contain the mass number in the 3rd to 5th position." print "It is assumed that this is the proton and we will change" print "the name to 'h 1' to be consistent with the notation used" print "in iniab.dat files" names[i]='h 1' mass_number[i]=int(names[i][2:5]) # now zip them together: hash_abu={} hash_index={} for a,b in zip(names,yps): hash_abu[a] = b for i in range(len(names)): hash_index[names[i]] = i self.z=z self.abu=yps self.a=mass_number self.names=names self.habu=hash_abu self.hindex=hash_index def write(self, outfile='initial_abundance.dat', header_string='initial abundances for a PPN run'): ''' Write initial abundance file (intended for use with ppn) Parameters ---------- outfile : string Name of output file. The default is 'initial_abundance.dat'. header_string : string A string with header line. The default is 'initial abundances for a PPN run'. ''' dcols=['Z', 'species','mass fraction'] data=[self.z,self.names,self.abu] hd=[header_string] att.write(outfile,hd,dcols,data) def write_mesa(self, mesa_isos_file='isos.txt', add_excess_iso='fe56', outfile='xa_iniabu.dat', header_string='initial abundances for a MESA run', header_char='!'): ''' Write initial abundance file, returns written abundances and mesa names. Parameters ---------- mesa_isos_file : string, optional List with isos copied from mesa network definition file in mesa/data/net_data/nets. The default is 'isos.txt'. add_excess_iso : string, optional Add 1.-sum(isos in mesa net) to this isotope. The defualt is 'fe56'. outfile : string, optional name of output file. The default file is 'xa_iniabu.dat'. header_string : string, optional Srting with header line. The default is 'initial abundances for a MESA run'. header_char : character, optional The default is '!'. Examples -------- >>> import utils >>> !ls ~/PPN/forum.astro.keele.ac.uk/frames/mppnp/USEEPP/ # find ppn initial abundance file >>> !cat ~/mesa/data/net_data/nets/agb.net # find isos needed in mesa net >>> !cat > isos.txt # paste needed isos into file >>> help(utils.iniabu) # check documentation of method >>> x=utils.iniabu('path_to_here/forum.astro.keele.ac.uk/frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn') >>> x.write_mesa? >>> mnames,mabus = x.write_mesa(add_excess_iso='ne22', ... header_string='mppnp/USEEPP/iniab2.0E-02GN93.ppn for mesa/agb.net', ... outfile='xa_2.0E-02GN93.mesa') ''' f=open('isos.txt') a=f.readlines() isos=[] for i in range(len(a)): isos.append(a[i].strip().rstrip(',')) mesa_names=[] abus=[] for i in range(len(self.z)): b=self.names[i].split() a='' a=a.join(b) if a in isos: mesa_names.append(a) abus.append(self.abu[i]) # mesa_names.append(elements_names[int(x.z[i])].lower()+str(int(x.a[i]))) for i in range(len(isos)): if isos[i] not in mesa_names: mesa_names.append(isos[i]) abus.append(0.0) excess=1.-np.sum(np.array(abus)) abus=np.array(abus) abus[mesa_names.index(add_excess_iso)] += excess dcols=['',''] data=[mesa_names,abus] hd=[header_string] att.write(outfile,hd,dcols,data,header_char=header_char) return mesa_names,abus def set_and_normalize(self,species_hash): ''' species_hash is a hash array in which you provide abundances referenced by species names that you want to set to some particular value; all other species are then normalised so that the total sum is 1. Examples -------- You can set up the argument array for this method for example in the following way. >>> sp={} >>> sp['he 4']=0.2 >>> sp['h 1']=0.5 ''' sum_before = 1. for i in range(len(species_hash)): sum_before -= self.abu[self.hindex[species_hash.keys()[i]]] print "sum_before = "+str(sum_before) normalization_factor=(1.-sum(species_hash.values()))/sum_before print "normalizing the rest witih factor "+str(normalization_factor) self.abu *= normalization_factor for i in range(len(species_hash)): self.abu[self.hindex[species_hash.keys()[i]]]=species_hash.values()[i] for name in self.habu: self.habu[name]=self.abu[self.hindex[name]] def isoratio_init(self,isos): ''' This file returns the isotopic ratio of two isotopes specified as iso1 and iso2. The isotopes are given as, e.g., ['Fe',56,'Fe',58] or ['Fe-56','Fe-58'] (for compatibility) -> list. ''' if len(isos) == 2: dumb = [] dumb = isos[0].split('-') dumb.append(isos[1].split('-')[0]) dumb.append(isos[1].split('-')[1]) isos = dumb ssratio = self.habu[isos[0].ljust(2).lower() + str(int(isos[1])).rjust(3)] / self.habu[isos[2].ljust(2).lower() + str(int(isos[3])).rjust(3)] return ssratio def iso_abundance(self,isos): ''' This routine returns the abundance of a specific isotope. Isotope given as, e.g., 'Si-28' or as list ['Si-28','Si-29','Si-30'] ''' if type(isos) == list: dumb = [] for it in range(len(isos)): dumb.append(isos[it].split('-')) ssratio = [] isos = dumb for it in range(len(isos)): ssratio.append(self.habu[isos[it][0].ljust(2).lower() + str(int(isos[it][1])).rjust(3)]) else: isos = isos.split('-') ssratio = self.habu[isos[0].ljust(2).lower() + str(int(isos[1])).rjust(3)] return ssratio def trajectory_SgConst(Sg=0.1, delta_logt_dex=-0.01): ''' setup trajectories for constant radiation entropy. S_gamma/R where the radiation constant R = N_A*k (Dave Arnett, Supernova book, p. 212) This relates rho and T but the time scale for this is independent. Parameters ---------- Sg : float S_gamma/R, values between 0.1 and 10. reflect conditions in massive stars. The default is 0.1. delta_logt_dex : float Sets interval between time steps in dex of logtimerev. The default is -0.01. ''' # reverse logarithmic time logtimerev=np.arange(5.,-6.,delta_logt_dex) logrho=np.linspace(0,8.5,len(logtimerev)) logT = (1./3.)*(logrho + 21.9161 + np.log10(Sg)) #rho_6=10**logrho/(0.1213*1.e6) #T9=rho_6**(1./3.) #logT_T3=np.log10(T9*1.e9) pl.close(3);pl.figure(3);pl.plot(logrho,logT,label='$S/\mathrm{N_Ak}='+str(Sg)+'$') pl.legend(loc=2);pl.xlabel('$\log \\rho$'); pl.ylabel('$\log T$') pl.close(5);pl.figure(5);pl.plot(logtimerev, logrho) pl.xlabel('$\log (t_\mathrm{final}-t)$'); pl.ylabel('$\log \\rho$') pl.xlim(8,-6) pl.close(6);pl.figure(6);pl.plot(logtimerev) pl.ylabel('$\log (t_\mathrm{final}-t)$'); pl.xlabel('cycle') # [t] logtimerev yrs # [rho] cgs # [T] K T9=10**logT/1.e9 data=[logtimerev,T9,logrho] att.writeTraj(filename='trajectory.input', data=data, ageunit=2, tunit=1, rhounit=1, idNum=1) # data: A list of 1D data vectors with time, T and rho # ageunit: If 1 ageunit = SEC, If 0 ageunit = YRS. If 2 agunit = logtimerev in yrs. Default is 0 # logtimerev is log of time until end def xlimrev(self): ''' reverse xrange''' xmax,xmin=pyl.xlim() pyl.xlim(xmin,xmax) def close_wins(win_min,win_max): ''' close all windows in a certain window number range win_min/max minumum and maximum window number to close ''' for i in range(win_min,win_max+1): close(i) def xlimrev(): ''' reverse xrange''' xmax,xmin=pl.xlim() pl.xlim(xmin,xmax) def species_list(what_list): ''' provide default lists of elements to plot. what_list : string String name of species lists provided. If what_list is "CNONe", then C, N, O and some other light elements. If what_list is "s-process", then s-process indicators. ''' if what_list is "CNONe": list_to_print = ['H-1','He-4','C-12','N-14','O-16','Ne-20'] elif what_list is "sprocess": list_to_print = ['Fe-56','Ge-70','Zn-70','Se-76','Kr-80','Kr-82','Kr-86','Sr-88','Ba-138','Pb-208'] elif what_list is "burn_stages": list_to_print = ['H-1','He-4','C-12','O-16','Ne-20','Si-28'] elif what_list is "list_marco_1": list_to_print = ['C-12','O-16','Ne-20','Ne-22','Na-23','Fe-54','Fe-56','Zn-70','Ge-70','Se-76','Kr-80','Kr-82','Sr-88','Y-89','Zr-96','Te-124','Xe-130','Xe-134','Ba-138'] return list_to_print def linestyle(i,a=5,b=3): ''' provide one out of 25 unique combinations of style, color and mark use in combination with markevery=a+mod(i,b) to add spaced points, here a would be the base spacing that would depend on the data density, modulated with the number of lines to be plotted (b) Parameters ---------- i : integer Number of linestyle combination - there are many.... a : integer Spacing of marks. The default is 5. b : integer Modulation in case of plotting many nearby lines. The default is 3. Examples -------- >>> plot(x,sin(x,linestyle(7)[0], markevery=linestyle(7)[1]) (c) 2014 FH ''' lines=['-','--','-.',':'] points=['v','^','<','>','1','2','3','4','s','p','*','h','H','+','x','D','d','o'] colors=['b','g','r','c','m','k'] ls_string = colors[sc.mod(i,6)]+lines[sc.mod(i,4)]+points[sc.mod(i,18)] mark_i = a+sc.mod(i,b) return ls_string,mark_i def symbol_list(what_list): ''' provide default symbol lists Parameters ---------- what_list : string String name of symbol lists provided; "list1", "list2", "lines1" or "lines2". ''' if what_list is "list1": symbol=['ro','bo','ko','go','mo'\ ,'r-','b-','k-','g-','m-','r--','b--','k--'\ ,'g--','r1'] #symbol=['r+','ro','r-'] elif what_list is "list2": symbol=['r-','b--','g-.','k:','md','.','o','v','^','<','>','1','2',\ '3','4','s','p','*','h','H','+'] elif what_list is "lines1": symbol=['b--','k--','r--','c--','m--','g--','b-','k-','r-','c-','m-','g-','b.','b-.','k-.','r-.','c-.','m-.','g-.','b:','k:','r:','c:','m:','g:'] elif what_list is "lines2": symbol=['g:','r-.','k-','b--','k-.','b+','r:','b-','c--','m--','g--','r-','c-','m-','g-','k-.','c-.','m-.','g-.','k:','r:','c:','m:','b-.','b:'] return symbol def make_list(default_symbol_list, len_list_to_print): ''' provide the list of symbols to use according for the list of species/arrays to plot. Parameters ---------- default_symbol_list : list Symbols that the user choose to use. len_list_to_print : integer len of list of species/arrays to print. ''' symbol_used = [] for i in range(len_list_to_print): symbol_used.append(default_symbol_list[sc.mod(i,len(default_symbol_list))]) return symbol_used def strictly_monotonic(bb): ''' bb is an index array which may have numerous double or triple occurrences of indices, such as for example the decay_index_pointer. This method removes all entries <= -, then all dublicates and finally returns a sorted list of indices. ''' cc=bb[np.where(bb>=0)] cc.sort() dc=cc[1:]-cc[:-1] # subsequent equal entries have 0 in db dc=np.insert(dc,0,1) # the first element is always unique (the second occurence is the dublicate) dc_mask=np.ma.masked_equal(dc,0) return np.ma.array(cc,mask=dc_mask.mask).compressed() def solar(filename_solar, solar_factor): ''' read solar abundances from filename_solar. Parameters ---------- filename_solar : string The file name. solar_factor : float The correction factor to apply, in case filename_solar is not solar, but some file used to get initial abundances at metallicity lower than solar. However, notice that this is really rude, since alpha-enahncements and things like that are not properly considered. Only H and He4 are not multiplied. So, for publications PLEASE use proper filename_solar at...solar, and use solar_factor = 1. Marco ''' f0=open(filename_solar) sol=f0.readlines() f0.close sol[0].split(" ") # Now read in the whole file and create a hashed array: global names_sol names_sol=[] global z_sol z_sol=[] yps=np.zeros(len(sol)) mass_number=np.zeros(len(sol)) for i in range(len(sol)): z_sol.append(int(sol[i][1:3])) names_sol.extend([sol[i].split(" ")[0][4:]]) yps[i]=float(sol[i].split(" ")[1]) * solar_factor try: mass_number[i]=int(names_sol[i][2:5]) except ValueError: print "WARNING:" print "This initial abundance file uses an element name that does" print "not contain the mass number in the 3rd to 5th position." print "It is assumed that this is the proton and we will change" print "the name to 'h 1' to be consistent with the notation used in" print "iniab.dat files" names_sol[i]='h 1' mass_number[i]=int(names_sol[i][2:5]) if mass_number[i] == 1 or mass_number[i] == 4: yps[i] = yps[i]/solar_factor # convert 'h 1' in prot, not needed any more?? #names_sol[0] = 'prot ' # now zip them together: global solar_abundance solar_abundance={} for a,b in zip(names_sol,yps): solar_abundance[a] = b z_bismuth = 83 global solar_elem_abund solar_elem_abund = np.zeros(z_bismuth) for i in range(z_bismuth): dummy = 0. for j in range(len(solar_abundance)): if z_sol[j] == i+1: dummy = dummy + float(solar_abundance[names_sol[j]]) solar_elem_abund[i] = dummy def convert_specie_naming_from_h5_to_ppn(isotope_names): ''' read isotopes names from h5 files, and convert them according to standard scheme used inside ppn and mppnp. Also Z and A are recalculated, for these species. Isomers are excluded for now, since there were recent changes in isomers name. As soon as the isomers names are settled, than Z and A provided here will be obsolete, and can be changed by usual Z and A. ''' spe_rude1 = [] spe_rude2 = [] spe_rude3 = [] for i in range(len(isotope_names)): spe_rude1.append(isotope_names[i].split('-')[0]) spe_rude2.append(isotope_names[i].split('-')[1]) # spe_rude1 is elem name and spe_rude2 is mass number. #print spe_rude1,spe_rude2 k = 0 for i in range(len(spe_rude1)): try: if int(spe_rude2[i]) < 10: spe_rude3.append(str(spe_rude1[i][0:2])+str(' ')+str(spe_rude2[i][0:3])) elif int(spe_rude2[i]) >= 10 and int(spe_rude2[i]) < 100 : spe_rude3.append(str(spe_rude1[i][0:2])+str(' ')+str(spe_rude2[i][0:3])) elif int(spe_rude2[i]) >= 100 : spe_rude3.append(str(spe_rude1[i][0:2])+str(spe_rude2[i][0:3])) except ValueError: k = k+1 None global spe spe = [] global n_array n_array = [] for i in range(len(spe_rude3)): if len(str(spe_rude1[i])) == 1: spe.append(str(spe_rude3[i][0:1])+str(' ')+str(spe_rude3[i][1:4])) else: spe.append(spe_rude3[i]) n_array.append(i) if spe[0]=='Ne 1': spe[0] = 'N 1' # spe_rude is the isotope name, in agreement with what we use in ppn, etc. # need to do this to can use other functions without changing them drastically. # here I skip isomers... global amass_int amass_int=np.zeros(len(spe_rude2)) for i in range(len(spe_rude2)-k): amass_int[i]=int(spe_rude2[i]) #print amass_int # here I have to create an array for the atomic numbers. # I need to this when I calculate and plot element abundances global znum_int znum_int=np.zeros(len(spe)) for i in range(len(spe)): if str(spe[i][0:2]) == 'H ': znum_int[i] = 1 elif str(spe[i][0:2]) == 'He': znum_int[i] = 2 elif str(spe[i][0:2]) == 'Li': znum_int[i] = 3 elif str(spe[i][0:2]) == 'Be': znum_int[i] = 4 elif str(spe[i][0:2]) == 'B ': znum_int[i] = 5 elif str(spe[i][0:2]) == 'C ': znum_int[i] = 6 elif str(spe[i][0:2]) == 'N ': znum_int[i] = 7 elif str(spe[i][0:2]) == 'O ': znum_int[i] = 8 elif str(spe[i][0:2]) == 'F ': znum_int[i] = 9 elif str(spe[i][0:2]) == 'Ne': znum_int[i] = 10 elif str(spe[i][0:2]) == 'Na': znum_int[i] = 11 elif str(spe[i][0:2]) == 'Mg': znum_int[i] = 12 elif str(spe[i][0:2]) == 'Al': znum_int[i] = 13 elif str(spe[i][0:2]) == 'Si': znum_int[i] = 14 elif str(spe[i][0:2]) == 'P ': znum_int[i] = 15 elif str(spe[i][0:2]) == 'S ': znum_int[i] = 16 elif str(spe[i][0:2]) == 'Cl': znum_int[i] = 17 elif str(spe[i][0:2]) == 'Ar': znum_int[i] = 18 elif str(spe[i][0:2]) == 'K ': znum_int[i] = 19 elif str(spe[i][0:2]) == 'Ca': znum_int[i] = 20 elif str(spe[i][0:2]) == 'Sc': znum_int[i] = 21 elif str(spe[i][0:2]) == 'Ti': znum_int[i] = 22 elif str(spe[i][0:2]) == 'V ': znum_int[i] = 23 elif str(spe[i][0:2]) == 'Cr': znum_int[i] = 24 elif str(spe[i][0:2]) == 'Mn': znum_int[i] = 25 elif str(spe[i][0:2]) == 'Fe': znum_int[i] = 26 elif str(spe[i][0:2]) == 'Co': znum_int[i] = 27 elif str(spe[i][0:2]) == 'Ni': znum_int[i] = 28 elif str(spe[i][0:2]) == 'Cu': znum_int[i] = 29 elif str(spe[i][0:2]) == 'Zn': znum_int[i] = 30 elif str(spe[i][0:2]) == 'Ga': znum_int[i] = 31 elif str(spe[i][0:2]) == 'Ge': znum_int[i] = 32 elif str(spe[i][0:2]) == 'As': znum_int[i] = 33 elif str(spe[i][0:2]) == 'Se': znum_int[i] = 34 elif str(spe[i][0:2]) == 'Br': znum_int[i] = 35 elif str(spe[i][0:2]) == 'Kr': znum_int[i] = 36 elif str(spe[i][0:2]) == 'Rb': znum_int[i] = 37 elif str(spe[i][0:2]) == 'Sr': znum_int[i] = 38 elif str(spe[i][0:2]) == 'Y ': znum_int[i] = 39 elif str(spe[i][0:2]) == 'Zr': znum_int[i] = 40 elif str(spe[i][0:2]) == 'Nb': znum_int[i] = 41 elif str(spe[i][0:2]) == 'Mo': znum_int[i] = 42 elif str(spe[i][0:2]) == 'Tc': znum_int[i] = 43 elif str(spe[i][0:2]) == 'Ru': znum_int[i] = 44 elif str(spe[i][0:2]) == 'Rh': znum_int[i] = 45 elif str(spe[i][0:2]) == 'Pd': znum_int[i] = 46 elif str(spe[i][0:2]) == 'Ag': znum_int[i] = 47 elif str(spe[i][0:2]) == 'Cd': znum_int[i] = 48 elif str(spe[i][0:2]) == 'In': znum_int[i] = 49 elif str(spe[i][0:2]) == 'Sn': znum_int[i] = 50 elif str(spe[i][0:2]) == 'Sb': znum_int[i] = 51 elif str(spe[i][0:2]) == 'Te': znum_int[i] = 52 elif str(spe[i][0:2]) == 'I ': znum_int[i] = 53 elif str(spe[i][0:2]) == 'Xe': znum_int[i] = 54 elif str(spe[i][0:2]) == 'Cs': znum_int[i] = 55 elif str(spe[i][0:2]) == 'Ba': znum_int[i] = 56 elif str(spe[i][0:2]) == 'La': znum_int[i] = 57 elif str(spe[i][0:2]) == 'Ce': znum_int[i] = 58 elif str(spe[i][0:2]) == 'Pr': znum_int[i] = 59 elif str(spe[i][0:2]) == 'Nd': znum_int[i] = 60 elif str(spe[i][0:2]) == 'Pm': znum_int[i] = 61 elif str(spe[i][0:2]) == 'Sm': znum_int[i] = 62 elif str(spe[i][0:2]) == 'Eu': znum_int[i] = 63 elif str(spe[i][0:2]) == 'Gd': znum_int[i] = 64 elif str(spe[i][0:2]) == 'Tb': znum_int[i] = 65 elif str(spe[i][0:2]) == 'Dy': znum_int[i] = 66 elif str(spe[i][0:2]) == 'Ho': znum_int[i] = 67 elif str(spe[i][0:2]) == 'Er': znum_int[i] = 68 elif str(spe[i][0:2]) == 'Tm': znum_int[i] = 69 elif str(spe[i][0:2]) == 'Yb': znum_int[i] = 70 elif str(spe[i][0:2]) == 'Lu': znum_int[i] = 71 elif str(spe[i][0:2]) == 'Hf': znum_int[i] = 72 elif str(spe[i][0:2]) == 'Ta': znum_int[i] = 73 elif str(spe[i][0:2]) == 'W ': znum_int[i] = 74 elif str(spe[i][0:2]) == 'Re': znum_int[i] = 75 elif str(spe[i][0:2]) == 'Os': znum_int[i] = 76 elif str(spe[i][0:2]) == 'Ir': znum_int[i] = 77 elif str(spe[i][0:2]) == 'Pt': znum_int[i] = 78 elif str(spe[i][0:2]) == 'Au': znum_int[i] = 79 elif str(spe[i][0:2]) == 'Hg': znum_int[i] = 80 elif str(spe[i][0:2]) == 'Tl': znum_int[i] = 81 elif str(spe[i][0:2]) == 'Pb': znum_int[i] = 82 elif str(spe[i][0:2]) == 'Bi': znum_int[i] = 83 elif str(spe[i][0:2]) == 'Po': znum_int[i] = 84 elif str(spe[i][0:2]) == 'At': znum_int[i] = 85 if spe[0] == 'N 1': znum_int[0] = 0 # here the index to connect name and atomic numbers. global index_atomic_number index_atomic_number = {} for a,b in zip(spe,znum_int): index_atomic_number[a]=b def define_zip_index_for_species(names_ppn_world, number_names_ppn_world): ''' This just give back cl, that is the original index as it is read from files from a data file.''' #connect the specie number in the list, with the specie name global cl cl={} for a,b in zip(names_ppn_world,number_names_ppn_world): cl[a] = b def element_abund_marco(i_decay, stable_isotope_list, stable_isotope_identifier, mass_fractions_array_not_decayed, mass_fractions_array_decayed): ''' Given an array of isotopic abundances not decayed and a similar array of isotopic abundances not decayed, here elements abundances, and production factors for elements are calculated ''' # this way is done in a really simple way. May be done better for sure, in a couple of loops. # I keep this, since I have only to copy over old script. Falk will probably redo it. import numpy as np #import utils as u global elem_abund elem_abund = np.zeros(z_bismuth) global elem_abund_decayed elem_abund_decayed = np.zeros(z_bismuth) global elem_prod_fac elem_prod_fac = np.zeros(z_bismuth) global elem_prod_fac_decayed elem_prod_fac_decayed = np.zeros(z_bismuth) # notice that elem_abund include all contribution, both from stables and unstables in # that moment. for i in range(z_bismuth): dummy = 0. for j in range(len(spe)): if znum_int[j] == i+1 and stable_isotope_identifier[j] > 0.5: dummy = dummy + float(mass_fractions_array_not_decayed[j]) elem_abund[i] = dummy for i in range(z_bismuth): if index_stable[i] == 1: elem_prod_fac[i] = float(elem_abund[i]/solar_elem_abund[i]) elif index_stable[i] == 0: elem_prod_fac[i] = 0. if i_decay == 2: for i in range(z_bismuth): dummy = 0. for j in range(len(mass_fractions_array_decayed)): if znum_int[cl[stable_isotope_list[j].capitalize()]] == i+1: #print znum_int[cl[stable[j].capitalize()]],cl[stable[j].capitalize()],stable[j] dummy = dummy + float(mass_fractions_array_decayed[j]) elem_abund_decayed[i] = dummy for i in range(z_bismuth): if index_stable[i] == 1: elem_prod_fac_decayed[i] = float(elem_abund_decayed[i]/solar_elem_abund[i]) elif index_stable[i] == 0: elem_prod_fac_decayed[i] = 0. def stable_specie(): ''' provide the list of stable species, and decay path feeding stables ''' import numpy as np stable_raw=[] stable_raw = ['H 1', 'H 2',\ 'HE 3', 'HE 4',\ 'LI 6', 'LI 7',\ 'BE 9',\ 'B 10', 'B 11',\ 'C 12', 'C 13',\ 'N 14', 'N 15',\ 'O 16', 'O 17', 'O 18',\ 'F 19',\ 'NE 20', 'NE 21', 'NE 22',\ 'NA 23',\ 'MG 24', 'MG 25', 'MG 26',\ 'AL 27',\ 'SI 28', 'SI 29', 'SI 30',\ 'P 31',\ 'S 32', 'S 33', 'S 34', 'S 36',\ 'CL 35', 'CL 37',\ 'AR 36', 'AR 38', 'AR 40',\ 'K 39', 'K 40', 'K 41',\ 'CA 40', 'CA 42', 'CA 43', 'CA 44', 'CA 46', 'CA 48',\ 'SC 45',\ 'TI 46', 'TI 47', 'TI 48', 'TI 49', 'TI 50',\ 'V 50', 'V 51',\ 'CR 50', 'CR 52', 'CR 53', 'CR 54',\ 'MN 55',\ 'FE 54', 'FE 56', 'FE 57', 'FE 58',\ 'CO 59',\ 'NI 58', 'NI 60', 'NI 61', 'NI 62', 'NI 64',\ 'CU 63', 'CU 65',\ 'ZN 64', 'ZN 66', 'ZN 67', 'ZN 68', 'ZN 70',\ 'GA 69', 'GA 71',\ 'GE 70', 'GE 72', 'GE 73', 'GE 74', 'GE 76',\ 'AS 75',\ 'SE 74', 'SE 76', 'SE 77', 'SE 78', 'SE 80', 'SE 82',\ 'BR 79', 'BR 81',\ 'KR 78', 'KR 80', 'KR 82', 'KR 83', 'KR 84', 'KR 86',\ 'RB 85', 'RB 87',\ 'SR 84', 'SR 86', 'SR 87', 'SR 88',\ 'Y 89',\ 'ZR 90', 'ZR 91', 'ZR 92', 'ZR 94', 'ZR 96',\ 'NB 93',\ 'MO 92', 'MO 94', 'MO 95', 'MO 96', 'MO 97', 'MO 98', 'MO100',\ 'RU 96', 'RU 98', 'RU 99', 'RU100', 'RU101', 'RU102', 'RU104',\ 'RH103',\ 'PD102', 'PD104', 'PD105', 'PD106', 'PD108', 'PD110',\ 'AG107', 'AG109',\ 'CD106', 'CD108', 'CD110', 'CD111', 'CD112', 'CD113', 'CD114', 'CD116',\ 'IN113', 'IN115',\ 'SN112', 'SN114', 'SN115', 'SN116', 'SN117', 'SN118', 'SN119', 'SN120', 'SN122', 'SN124',\ 'SB121', 'SB123',\ 'TE120', 'TE122', 'TE123', 'TE124', 'TE125', 'TE126', 'TE128', 'TE130',\ 'I 127',\ 'XE124', 'XE126', 'XE128', 'XE129', 'XE130', 'XE131', 'XE132', 'XE134', 'XE136',\ 'CS133',\ 'BA130', 'BA132', 'BA134', 'BA135', 'BA136', 'BA137', 'BA138',\ 'LA138', 'LA139',\ 'CE136', 'CE138', 'CE140', 'CE142',\ 'PR141',\ 'ND142', 'ND143', 'ND144', 'ND145', 'ND146', 'ND148', 'ND150',\ 'SM144', 'SM147', 'SM148', 'SM149', 'SM150', 'SM152', 'SM154',\ 'EU151', 'EU153',\ 'GD152', 'GD154', 'GD155', 'GD156', 'GD157', 'GD158', 'GD160',\ 'TB159',\ 'DY156', 'DY158', 'DY160', 'DY161', 'DY162', 'DY163', 'DY164',\ 'HO165',\ 'ER162', 'ER164', 'ER166', 'ER167', 'ER168', 'ER170',\ 'TM169',\ 'YB168', 'YB170', 'YB171', 'YB172', 'YB173', 'YB174', 'YB176',\ 'LU175', 'LU176',\ 'HF174', 'HF176', 'HF177', 'HF178', 'HF179', 'HF180',\ 'TA180', 'TA181',\ 'W 180', 'W 182', 'W 183', 'W 184', 'W 186',\ 'RE185', 'RE187',\ 'OS184', 'OS186', 'OS187', 'OS188', 'OS189', 'OS190', 'OS192',\ 'IR191', 'IR193',\ 'PT190', 'PT192', 'PT194', 'PT195', 'PT196', 'PT198',\ 'AU197',\ 'HG196', 'HG198', 'HG199', 'HG200', 'HG201', 'HG202', 'HG204',\ 'TL203', 'TL205',\ 'PB204', 'PB206', 'PB207', 'PB208',\ 'BI209'] jj=-1 global count_size_stable count_size_stable=[] global stable stable=[] global jdum jdum=np.zeros(len(stable_raw)) global jjdum jjdum=np.zeros(len(spe)) for i in range(len(stable_raw)): dum_str = stable_raw[i] for j in range(len(spe)): if stable_raw[i].capitalize() == spe[j]: stable.append(stable_raw[i]) jdum[i]=1 jjdum[j]=1 jj=jj+1 count_size_stable.append(int(jj)) #print stable # back_ind is an index to go back, to use the order of stable # useful for example for decayed yields. global back_ind back_ind={} for a,b in zip(stable,count_size_stable): back_ind[a]=b #print 'in stable:',back_ind['SE 74'] # definition of decay paths global decay_raw decay_raw=[] decay_raw=[['H 1','N 1'],\ ['H 2'],\ ['HE 3'],\ ['HE 4','B 8'],\ ['LI 6'],\ ['LI 7','BE 7'],\ ['BE 9'],\ ['B 10','BE 10'],\ ['B 11','C 11','BE 11'],\ ['C 12'],\ ['C 13','N 13','O 13'],\ ['N 14','C 14','O 14'],\ ['N 15','C 15','O 15','F 15'],\ ['O 16'],\ ['O 17','F 17'],\ ['O 18','F 18','NE 18'],\ ['F 19','O 19','NE 19'],\ ['NE 20','F 20','NA 20'],\ ['NE 21','F 21','NA 21'],\ ['NE 22','NA 22','MG 22','F 22'],\ ['NA 23','MG 23','NE 23'],\ ['MG 24','AL 24','NA 24','NE 24',],\ ['MG 25','NA 25','AL 25'],\ ['MG 26','SI 26','AL 26','NA 26','AL*26'],\ ['AL 27','SI 27','MG 27'],\ ['SI 28','AL 28','MG 28'],\ ['SI 29','P 29','AL 29','MG 29'],\ ['SI 30','S 30','P 30','AL 30','MG 30'],\ ['P 31','S 31','SI 31'],\ ['S 32','P 32','SI 32'],\ ['S 33','CL 33','P 33','SI 33'],\ ['S 34','CL 34','P 34'],\ ['S 36','P 36'],\ ['CL 35','S 35','AR 35'],\ ['CL 37','S 37','P 37','AR 37','K 37'],\ ['AR 36','CL 36'],\ ['AR 38','CL 38','S 38','P 38','K 38'],\ ['AR 40','CL 40','S 40'],\ ['K 39','AR 39','CL 39','S 39','CA 39'],\ ['K 40'],\ ['K 41','AR 41','CL 41','S 41','CA 41','SC 41'],\ ['CA 40','SC 40'],\ ['CA 42','K 42','AR 42','CL 42','S 42','SC 42','TI 42'],\ ['CA 43','K 43','AR 43','CL 43','SC 43','TI 43','V 43'],\ ['CA 44','K 44','AR 44','CL 44','SC 44','TI 44'],\ ['CA 46','K 46','AR 46'],\ ['CA 48','K 48','AR 48'],\ ['SC 45','CA 45','K 45','AR 45','CL 45','TI 45','V 45'],\ ['TI 46','SC 46','V 46','CR 46'],\ ['TI 47','SC 47','CA 47','K 47','AR 47','V 47','CR 47'],\ ['TI 48','SC 48','V 48','CR 48'],\ ['TI 49','SC 49','CA 49','K 49','V 49','CR 49','MN 49'],\ ['TI 50','SC 50','CA 50','K 50'],\ ['V 50'],\ ['V 51','CR 51','TI 51','SC 51','CA 51','MN 51'],\ ['CR 50','MN 50'],\ ['CR 52','MN 52','FE 52','V 52','TI 52','SC 52','CA 52'],\ ['CR 53','MN 53','FE 53','V 53','TI 53','SC 53'],\ ['CR 54','MN 54','V 54','TI 54','SC 54'],\ ['MN 55','FE 55','CR 55','V 55','TI 55','CO 55'],\ ['FE 54','CO 54'],\ ['FE 56','NI 56','CO 56','MN 56','CR 56'],\ ['FE 57','NI 57','CO 57','MN 57','CR 57'],\ ['FE 58','CO 58','MN 58','CR 58'],\ ['CO 59','FE 59','MN 59','CR 59','NI 59','CU 59'],\ ['NI 58','CU 58'],\ ['NI 60','CO 60','FE 60','MN 60','CR 60','CU 60','ZN 60'],\ ['NI 61','CO 61','FE 61','MN 61','CU 61','ZN 61'],\ ['NI 62','CO 62','FE 62','CU 62','ZN 62'],\ ['NI 64','CO 64','FE 64','CU 64'],\ ['CU 63','NI 63','CO 63','FE 63','MN 63','ZN 63','GA 63'],\ ['CU 65','NI 65','CO 65','FE 65','ZN 65','GA 65','GE 65'],\ ['ZN 64','CU 64','GA 64','GE 64'],\ ['ZN 66','CU 66','NI 66','CO 66','FE 66','GA 66','GE 66'],\ ['ZN 67','CU 67','NI 67','CO 67','FE 67','GA 67','GE 67','AS 77'],\ ['ZN 68','NI 68','CO 68','GA 68','GE 68','CU 68','AS 68','SE 68'],\ ['ZN 70','CU 70','NI 70','CO 70'],\ ['GA 69','ZN 69','CU 69','NI 69','GE 69','AS 69','SE 69'],\ ['GA 71','ZN 71','CU 71','NI 71','GE 71','AS 71','SE 71','BR 71'],\ ['GE 70','GA 70','AS 70','SE 70','BR 70'],\ ['GE 72','GA 72','ZN 72','CU 72','NI 72','AS 72','SE 72','BR 72','KR 72'],\ ['GE 73','GA 73','ZN 73','CU 73','NI 73','AS 73','SE 73','BR 73','KR 73'],\ ['GE 74','GA 74','ZN 74','CU 74','NI 74','AS 74'],\ ['GE 76','GA 76','ZN 76','CU 76'],\ ['AS 75','GE 75','GA 75','ZN 75','CU 75','SE 75','BR 75','KR 75','RB 75'],\ ['SE 74','AS 74','BR 74','KR 74'],\ ['SE 76','AS 76','BR 76','KR 76','RB 76','SR 76'],\ ['SE 77','AS 77','GE 77','BR 77','GA 77','ZN 77','KR 77','RB 77','SR 77'],\ ['SE 78','AS 78','GE 78','GA 78','ZN 78','BR 78'],\ ['SE 80','AS 80','GE 80','GA 80','ZN 80'],\ ['SE 82','AS 82','GE 82','GA 82'],\ ['BR 79','SE 79','AS 79','GE 79','GA 79','ZN 79','KR 79','RB 79','SR 79','Y 79'],\ ['BR 81','SE 81','KR 81','AS 81','GE 81','GA 81','RB 81','SR 81','Y 81','ZR 81'],\ ['KR 78','RB 78','SR 78','Y 78'],\ ['KR 80','BR 80','RB 80','SR 80','ZR 80'],\ ['KR 82','BR 82','RB 82','SR 82','Y 82','ZR 82'],\ ['KR 83','BR 83','SE 83','AS 83','GE 83','RB 83','SR 83','Y 83','ZR 83','NB 83'],\ ['KR 84','BR 84','SE 84','AS 84','GE 84','RB 84'],\ ['KR 86','BR 86','SE 86','AS 86'],\ ['RB 85','KR 85','SR 85','KR*85','BR 85','SE 85','AS 85','Y 85','ZR 85','NB 85','MO 85'],\ ['RB 87','KR 87','BR 87','SE 87','AS 87'],\ ['SR 84','Y 84','ZR 84','NB 84','MO 84'],\ ['SR 86','RB 86','Y 86','ZR 86','NB 86','MO 86'],\ ['SR 87','Y 87','ZR 87','NB 87','MO 87','TC 87'],\ ['SR 88','RB 88','KR 88','BR 88','SE 88','Y 88','ZR 88','NB 88','MO 88','TC 88'],\ ['Y 89','SR 89','RB 89','KR 89','BR 89','ZR 89','NB 89','MO 89','TC 89','RU 89'],\ ['ZR 90','Y 90','SR 90','RB 90','KR 90','BR 90','NB 90','MO 90','TC 90','RU 90'],\ ['ZR 91','Y 91','SR 91','RB 91','KR 91','BR 91','SE 91','NB 91','MO 91','TC 91','RU 91','RH 91'],\ ['ZR 92','Y 92','SR 92','RB 92','KR 92','BR 92','NB 92'],\ ['ZR 94','Y 94','SR 94','RB 94','KR 94'],\ ['ZR 96','Y 96','SR 96'],\ ['NB 93','ZR 93','Y 93','SR 93','RB 93','KR 93','MO 93','TC 93','RU 93','RH 93','PD 93'],\ ['MO 92','TC 92','RU 92','RH 92','PD 92'],\ ['MO 94','NB 94','TC 94','RU 94','RH 94','PD 94'],\ ['MO 95','NB 95','ZR 95','Y 95','SR 95','TC 95','RU 95','RH 95','PD 95','AG 95'],\ ['MO 96','NB 96','TC 96'],\ ['MO 97','NB 97','ZR 97','Y 97','SR 97','TC 97','RU 97','RH 97','PD 97','AG 97','CD 97'],\ ['MO 98','NB 98','ZR 98','Y 98','SR 98'],\ ['MO100','NB100','ZR100','Y 100'],\ ['RU 96','RH 96','PD 96','AG 96'],\ ['RU 98','TC 98','RH 98','PD 98','AG 98','CD 98'],\ ['RU 99','TC 99','MO 99','NB 99','ZR 99','Y 99','RH 99','PD 99','AG 99','CD 99','IN 99'],\ ['RU100','TC100','RH100','PD100','AG100','CD100','IN100','SN100'],\ ['RU101','TC101','MO101','NB101','ZR101','Y 101','RH101','PD101','AG101','CD101','IN101','SN101'],\ ['RU102','MO102','TC102','NB102','ZR102','Y 102','RH102'],\ ['RU104','TC104','MO104','NB104'],\ ['RH103','RU103','TC103','MO103','NB103','ZR103','Y 103','PD103','AG103','CD103','IN103','SN103'],\ ['PD102','AG102','CD102','IN102','SN102'],\ ['PD104','RH104','AG104','CD104','IN104','SN104','SB104'],\ ['PD105','RH105','RU105','TC105','MO105','NB105','ZR105','AG105','CD105','IN105','SN105','SB105'],\ ['PD106','RH106','RU106','TC106','MO106','NB106','AG106'],\ ['PD108','RH108','RU108','TC108','MO108','NB108'],\ ['PD110','RH110','RU110','TC110','MO110','NB110'],\ ['AG107','PD107','RH107','RU107','TC107','MO107','CD107','IN107','SN107','SB107'],\ ['AG109','PD109','RH109','RU109','TC109','MO109','NB109','CD109','IN109','SN109','SB109','TE109'],\ ['CD106','IN106','SN106','SB106'],\ ['CD108','AG108','IN108','SN108','SB108'],\ ['CD110','AG110','IN110','SN110','SB110','TE110'],\ ['CD111','AG111','PD111','RH111','RU111','TC111','IN111','SN111','SB111','TE111','I 111'],\ ['CD112','AG112','PD112','RH112','RU112','TC112'],\ ['CD113','AG113','PD113','RH113','RU113'],\ ['CD114','AG114','PD114','RH114','RU114'],\ ['CD116','AG116','PD116','RH116'],\ ['IN113','SN113','SB113','TE113','I 113','XE113'],\ ['IN115','CD115','AG115','PD115','RH115','RU115'],\ ['SN112','IN112','SB112','TE112','I 112','XE112'],\ ['SN114','IN114','SB114','TE114','I 114','XE114','CS114'],\ ['SN115','SB115','TE115','I 115','XE115','CS115','BA115'],\ ['SN116','IN116','SB116','TE116','I 116','XE116','CS116','BA116'],\ ['SN117','IN117','AG117','PD117','RH117','SB117','CD117','TE117','I 117','XE117','CS117','BA117'],\ ['SN118','IN118','CD118','AG118','PD118','SB118','TE118','XE118','CS118','BA118'],\ ['SN119','IN119','CD119','AG119','PD119','SB119','TE119','XE119','CS119','BA119'],\ ['SN120','SB120','IN120','CD120','AG120','PD120'],\ ['SN122','IN122','CD122','AG122'],\ ['SN124','IN124','CD124'],\ ['SB121','SN121','IN121','CD121','AG121','TE121','I 121','XE121','CS121','XE121','CS121','BA121','LA121','CE121'], ['SB123','SN123','IN123','CD123','AG123'],\ ['TE120','I 120','XE120','CS120','BA120','LA120'],\ ['TE122','SB122','I 122','XE122','CS122','BA122','LA122'],\ ['TE123','I 123','XE123','CS123','BA123','LA123','CE123'],\ ['TE124','SB124','I 124'],\ ['TE125','SB125','SN125','IN125','CD125','I 125','XE125','CS125','BA125','LA125','CE125','PR125'],\ ['TE126','SB126','SN126','IN126','CD126'],\ ['TE128','SB128','SN128','IN128','CD128'],\ ['TE130','SB130','SN130','IN130'],\ ['I 127','TE127','SB127','SN127','IN127','CD127','XE127','CS127','BA127','LA127','CE127','PR127','ND127'],\ ['XE124','CS124','BA124','LA124','CE124','PR124'],\ ['XE126','I 126','CS126','BA126','LA126','CE126','PR126'],\ ['XE128','I 128','CS128','BA128','LA128','CE128','PR128'],\ ['XE129','I 129','TE129','SB129','SN129','IN129','CD129','CS129','BA129','LA129','CE129','PR129','ND129'],\ ['XE130','I 130','CS130'],\ ['XE131','I 131','TE131','SB131','SN131','IN131','CS131','BA131','LA131','CE131','PR131','ND131','PM131','SM131'],\ ['XE132','I 132','TE132','SB132','SN132','IN132','CS132'], ['XE134','I 134','TE134','SB134','SN134'],\ ['XE136','I 136','TE136','SB136'],\ ['CS133','XE133','I 133','TE133','SB133','SN133','BA133','LA133','CE133','PR133','ND133','PM133','SM133'],\ ['BA130','LA130','CE130','PR130','ND130','PM130'],\ ['BA132','LA132','CE132','PR132','ND132','PM132','SM132'],\ ['BA134','CS134','LA134','CE134','PR134','ND134','PM134','SM134','EU134'],\ ['BA135','CS135','XE135','I 135','TE135','SB135','SN135','LA135','CE135','PR135','ND135','PM135','SM135','EU135','GD135'],\ ['BA136','CS136','LA136'],\ ['BA137','CS137','XE137','I 137','TE137','SB137','LA137','CE137','PR137','ND137','PM137','SM137','EU137','GD137'],\ ['BA138','CS138','XE138','I 138','TE138'],\ ['LA138'],\ ['LA139','BA139','CS139','XE139','I 139','CE139','PR139','ND139','PM139','SM139','EU139','GD139','TB139','DY139'],\ ['CE136','PR136','ND136','PM136','SM136','EU136'],\ ['CE138','PR138','ND138','PM138','SM138','EU138','GD138'],\ ['CE140','LA140','BA140','CS140','XE140','I 140','PR140','ND140','PM140','SM140','EU140','GD140','TB140'],\ ['CE142','LA142','BA142','CS142','XE142','I 142'],\ ['PR141','CE141','LA141','BA141','CS141','XE141','I 141','ND141','PM141','SM141','EU141','GD141','TB141','DY141'],\ ['ND142','PR142','PM142','SM142','EU142','GD142','TB142','DY142','HO142','SM146','EU146','GD146','TB146','DY146','HO146','ER146',\ 'GD150','TB150','DY150','HO150','ER150','DY154','HO154','ER154'],\ ['ND143','PR143','CE143','LA143','BA143','CS143','XE143','PM143','SM143','EU143','GD143','TB143','DY143'],\ ['ND144','PR144','CE144','LA144','BA144','CS144','XE144','PM144'],\ ['ND145','PR145','CE145','LA145','BA145','CS145','PM145','SM145','EU145','GD145','TB145','DY145','HO145','ER145'],\ ['ND146','PR146','CE146','LA146','BA146','CS146','PM146'],\ ['ND148','PR148','CE148','LA148','BA148'],\ ['ND150','PR150','CE150','LA150','BA150'],\ ['SM144','EU144','GD144','TB144','DY144','HO144','GD148','TB148','DY148','HO148','ER148','TM148'],\ ['SM147','PM147','ND147','PR147','CE147','LA147','BA147','EU147','GD147','TB147','DY147','HO147','ER147'],\ ['SM148','PM148','EU148'],\ ['SM149','PM149','ND149','PR149','CE149','LA149','EU149','GD149','TB149','DY149','HO149','ER149','TM149','YB149'],\ ['SM150','PM150','EU150'],\ ['SM152','PM152','ND152','PR152','CE152','EU152'],\ ['SM154','PM154','ND154','PR154'],\ ['EU151','SM151','PM151','ND151','PR151','CE151','GD151','TB151'],\ ['EU153','SM153','PM153','GD153','PR153','GD153','TB153','DY153','HO153'],\ ['GD152','TB152','DY152','HO152'],\ ['GD154','EU154','TB154'],\ ['GD155','EU155','SM155','PM155','ND155','TB155','DY155','HO155','ER155','TM155'],\ ['GD156','EU156','SM156','PM156','ND156','TB156'],\ ['GD157','EU157','SM157','PM157','TB157','DY157','HO157','ER157','TM157','YB157'],\ ['GD158','EU158','SM158','PM158','TB158'],\ ['GD160','EU160','SM160'],\ ['TB159','GD159','EU159','SM159','PM159','DY159','HO159','ER159','TM159','YB159','LU159'],\ ['DY156','HO156','ER156','TM156','YB156','HO156','ER156','TM156','YB156'],\ ['DY158','HO158','ER158','TM158','YB158','LU158'],\ ['DY160','TB160','HO160','ER160','TM160','YB160','LU160','HF160'],\ ['DY161','TB161','GD161','EU161','SM161','PM161','ND161','HO161','ER161','TM161','YB161','LU161','HF161','TA161'],\ ['DY162','TB162','GD162','EU162','SM162','PM162','HO162'],\ ['DY163','TB163','HO163','GD163','EU163','SM163','PM163','ND163','ER163','TM163','YB163','LU163','HF163','TA163'],\ ['DY164','TB164','HO164','GD164','EU164','SM164','PM164','ND164','HO164'],\ ['HO165','DY165','ER165','TB165','GD165','EU165','SM165','PM165','HO165','TM165','YB165','LU165','HF165','TA165','W 165'],\ ['ER162','TM162','YB162','LU162','HF162','TA162'],\ ['ER164','TM164','YB164','LU164','HF164','TA164','W 164'],\ ['ER166','HO166','DY166','TB166','GD166','EU166','SM166','PM166','ND166','TM166','YB166','LU166','HF166','TA166','W 166','RE166'],\ ['ER167','HO167','DY167','TB167','GD167','EU167','SM167','PM167','ND167','TM167','YB167','LU167','HF167','TA167','W 167','RE167'],\ ['ER168','HO168','DY168','TB168','GD168','EU168','SM168','PM168','ND168','TM168'],\ ['ER170','HO170','DY170','TB170','GD170','EU170','SM170','PM170','ND170'],\ ['TM169','ER169','HO169','DY169','TB169','GD169','EU169','SM169','PM169','ND169','YB169','LU169','HF169','TA169','W 169','RE169'],\ ['YB168','LU168','HF168','TA168','W 168','RE168'],\ ['YB170','TM170','LU170','HF170','TA170','W 170','RE170','OS170'],\ ['YB171','TM171','ER171','HO171','DY171','TB171','GD171','EU171','SM171','PM171','ND171','LU171','HF171','TA171','W 171','RE171','OS171'],\ ['YB172','TM172','ER172','HO172','DY172','TB172','GD172','EU172','SM172','PM172','ND172','LU172','HF172','TA172','W 172','RE172','OS172','IR172'],\ ['YB173','TM173','ER173','HO173','DY173','TB173','GD173','EU173','SM173','PM173','ND173','LU173','HF173','TA173','W 173','RE173','OS173','IR173'],\ ['YB174','TM174','ER174','HO174','DY174','TB174','GD174','EU174','SM174','PM174','ND174','LU174'],\ ['YB176','TM176','ER176','HO176','DY176','TB176','GD176','EU176','SM176','PM176','ND176'],\ ['LU175','YB175','TM175','ER175','HO175','DY175','TB175','GD175','EU175','SM175','PM175','ND175','HF175','TA175','W 175','RE175','OS175','IR175'],\ ['LU176','HF176','TA176','W 176','RE176','OS176','IR176'],\ ['HF174','TA174','W 174','RE174','OS174','IR174'],\ ['HF176','TA176','W 176','RE176','OS176','IR176'],\ ['HF177','LU177','YB177','LU177','YB177','TM177','ER177','HO177','DY177','TB177','GD177','TA177','W 177','RE177','OS177','IR177'],\ ['HF178','LU178','YB178','LU178','YB178','TM178','ER178','HO178','DY178','TB178','GD178','TA178','W 178','RE178','OS178','IR178'],\ ['HF179','LU179','YB179','LU179','YB179','TM179','ER179','HO179','DY179','TB179','GD179','TA179','W 179','RE179','OS179','IR179','PT179'],\ ['HF180','LU180','YB180','LU180','YB180','TM180','ER180','HO180','DY180','TB180','GD180','TA180','W 180','RE180','OS180','IR180','PT180','AU180'],\ ['TA180'],\ ['TA181','HF181','LU181','YB181','LU181','YB181','TM181','ER181','HO181','DY181','TB181','GD181','W 181','RE181','OS181','IR181','PT181','AU181'],\ ['W 180','RE180','OS180','IR180','PT180','AU180'],\ ['W 182','TA182','HF182','LU182','YB182','LU182','YB182','TM182','ER182','HO182','DY182','TB182','GD182','RE182','OS182','IR182','PT182','AU182'],\ ['W 183','TA183','HF183','LU183','YB183','LU183','YB183','TM183','ER183','HO183','DY183','TB183','GD183','RE183','OS183','IR183','PT183','AU183'],\ ['W 184','TA184','HF184','LU184','YB184','LU184','YB184','TM184','ER184','HO184','DY184','TB184','GD184','RE184'],\ ['W 186','TA186','HF186','LU186','YB186','LU186','YB186','TM186','ER186','HO186','DY186','TB186','GD186'],\ ['RE185','W 185','TA185','HF185','LU185','YB185','LU185','YB185','TM185','ER185','HO185','DY185','TB185','GD185','OS185','IR185','PT185','AU185','HG185','TL185'],\ ['RE187','W 187','TA187','HF187','LU187','YB187','LU187','YB187','TM187','ER187','HO187','DY187','TB187','GD187'],\ ['OS184','IR184','PT184','AU184','HG184','TL184'],\ ['OS186','RE186','IR186','PT186','AU186','HG186','TL186'],\ ['OS187','IR187','PT187','AU187','HG187','TL187','PB187'],\ ['OS188','RE188','W 188','TA188','HF188','LU188','YB188','LU188','YB188','TM188','ER188','HO188','DY188','TB188','GD188','IR188','PT188','AU188','HG188','TL188','PB188'],\ ['OS189','RE189','W 189','TA189','HF189','LU189','YB189','LU189','YB189','TM189','ER189','HO189','DY189','TB189','GD189','IR189','PT189','AU189','HG189','TL189','PB189'],\ ['OS190','RE190','W 190','TA190','HF190','LU190','YB190','LU190','YB190','TM190','ER190','HO190','DY190','TB190','GD190','IR190'],\ ['OS192','RE190','W 192','TA192','HF192','LU192','YB192','LU192','YB192','TM192','ER192','HO192','DY192','TB192','GD192'],\ ['IR191','OS191','RE191','W 191','TA191','HF191','LU191','YB191','LU191','YB191','TM191','ER191','HO191','DY191','TB191','GD191','PT191','AU191','HG191','TL191','PB191'],\ ['IR193','OS193','RE193','W 193','TA193','HF193','LU193','YB193','LU193','YB193','TM193','ER193','HO193','DY193','TB193','GD193','PT193','AU193','HG193','TL193','PB193'],\ ['PT190','AU190','HG190','TL190','PB190'],\ ['PT192','IR192','AU192','HG192','TL192','PB192'],\ ['PT194','IR194','OS194','RE194','W 194','TA194','HF194','AU194','HG194','TL194','PB194','BI194'], ['PT195','IR195','OS195','RE195','W 195','TA195','HF195','AU195','HG195','TL195','PB195','BI195'],\ ['PT196','IR196','OS196','RE196','W 196','TA196','HF196','AU196'],\ ['PT198','IR198','OS198','RE198','W 198','TA198','HF198'],\ ['AU197','PT197','IR197','OS197','RE197','W 197','TA197','HF197','HG197','TL197','PB197','BI197'],\ ['HG196','TL196','PB196','BI196'],\ ['HG198','AU198','TL198','PB198','BI198'],\ ['HG199','AU199','PT199','IR199','OS199','RE199','W 199','TA199','HF199','TL199','PB199','BI199'],\ ['HG200','AU200','PT200','IR200','OS200','RE200','W 200','TA200','HF200','TL200','PB200','BI200'],\ ['HG201','AU201','PT201','IR201','OS201','RE201','W 201','TA201','HF201','TL201','PB201','BI201','PO201'],\ ['HG202','AU202','PT202','IR202','OS202','RE202','W 202','TA202','HF202','TL202','PB202','BI202','PO202'],\ ['HG204','AU204','PT204','IR204','OS204','RE204','W 204','TA204','HF204'],\ ['TL203','HG203','AU203','PT203','IR203','OS203','RE203','W 203','TA203','HF203','PB203','BI203','PO203'],\ ['TL205','HG205','AU205','PT205','IR205','OS205','RE205','W 205','TA205','HF205','PB205','BI205','PO205'],\ ['PB204','TL204','BI204','PO204'],\ ['PB206','TL206','HG206','AU206','PT206','IR206','OS206','RE206','W 206','TA206','HF206','BI206','PO210'],\ ['PB207','TL207','HG207','AU207','PT207','IR207','OS207','RE207','W 207','TA207','HF207','BI207','PO211','BI211'],\ ['PB208','TL208','HG208','AU208','PT208','IR208','OS208','RE208','W 208','TA208','HF208','BI208','PO212','BI212'],\ ['BI209','PB209','TL209','HG209','AU209','PT209','IR209','OS209','RE209','W 209','TA209','HF209']] #print decay_raw def give_zip_element_z_and_names(element_name): ''' create 2 indexes that, given the name of the element/specie, give the atomic number.''' import numpy as np global z_bismuth z_bismuth = 83 global z_for_elem z_for_elem = [] global index_stable index_stable = [] i_for_stable = 1 i_for_unstable = 0 for i in range(z_bismuth): z_for_elem.append(int(i+1)) # the only elements below bismuth with no stable isotopes are Tc and Pm if i+1 == 43 or i+1 == 61: index_stable.append(i_for_unstable) else: index_stable.append(i_for_stable) dummy_index = np.zeros(len(element_name)) for i in range(len(element_name)): if element_name[i] == 'Neutron': dummy_index[i] = 0 elif element_name[i] == 'H': dummy_index[i] = 1 elif element_name[i] == 'He': dummy_index[i] = 2 elif element_name[i] == 'Li': dummy_index[i] = 3 elif element_name[i] == 'Be': dummy_index[i] = 4 elif element_name[i] == 'B': dummy_index[i] = 5 elif element_name[i] == 'C': dummy_index[i] = 6 elif element_name[i] == 'N': dummy_index[i] = 7 elif element_name[i] == 'O': dummy_index[i] = 8 elif element_name[i] == 'F': dummy_index[i] = 9 elif element_name[i] == 'Ne': dummy_index[i] = 10 elif element_name[i] == 'Na': dummy_index[i] = 11 elif element_name[i] == 'Mg': dummy_index[i] = 12 elif element_name[i] == 'Al': dummy_index[i] = 13 elif element_name[i] == 'Si': dummy_index[i] = 14 elif element_name[i] == 'P': dummy_index[i] = 15 elif element_name[i] == 'S': dummy_index[i] = 16 elif element_name[i] == 'Cl': dummy_index[i] = 17 elif element_name[i] == 'Ar': dummy_index[i] = 18 elif element_name[i] == 'K': dummy_index[i] = 19 elif element_name[i] == 'Ca': dummy_index[i] = 20 elif element_name[i] == 'Sc': dummy_index[i] = 21 elif element_name[i] == 'Ti': dummy_index[i] = 22 elif element_name[i] == 'V': dummy_index[i] = 23 elif element_name[i] == 'Cr': dummy_index[i] = 24 elif element_name[i] == 'Mn': dummy_index[i] = 25 elif element_name[i] == 'Fe': dummy_index[i] = 26 elif element_name[i] == 'Co': dummy_index[i] = 27 elif element_name[i] == 'Ni': dummy_index[i] = 28 elif element_name[i] == 'Cu': dummy_index[i] = 29 elif element_name[i] == 'Zn': dummy_index[i] = 30 elif element_name[i] == 'Ga': dummy_index[i] = 31 elif element_name[i] == 'Ge': dummy_index[i] = 32 elif element_name[i] == 'As': dummy_index[i] = 33 elif element_name[i] == 'Se': dummy_index[i] = 34 elif element_name[i] == 'Br': dummy_index[i] = 35 elif element_name[i] == 'Kr': dummy_index[i] = 36 elif element_name[i] == 'Rb': dummy_index[i] = 37 elif element_name[i] == 'Sr': dummy_index[i] = 38 elif element_name[i] == 'Y': dummy_index[i] = 39 elif element_name[i] == 'Zr': dummy_index[i] = 40 elif element_name[i] == 'Nb': dummy_index[i] = 41 elif element_name[i] == 'Mo': dummy_index[i] = 42 elif element_name[i] == 'Tc': dummy_index[i] = 43 elif element_name[i] == 'Ru': dummy_index[i] = 44 elif element_name[i] == 'Rh': dummy_index[i] = 45 elif element_name[i] == 'Pd': dummy_index[i] = 46 elif element_name[i] == 'Ag': dummy_index[i] = 47 elif element_name[i] == 'Cd': dummy_index[i] = 48 elif element_name[i] == 'In': dummy_index[i] = 49 elif element_name[i] == 'Sn': dummy_index[i] = 50 elif element_name[i] == 'Sb': dummy_index[i] = 51 elif element_name[i] == 'Te': dummy_index[i] = 52 elif element_name[i] == 'I': dummy_index[i] = 53 elif element_name[i] == 'Xe': dummy_index[i] = 54 elif element_name[i] == 'Cs': dummy_index[i] = 55 elif element_name[i] == 'Ba': dummy_index[i] = 56 elif element_name[i] == 'La': dummy_index[i] = 57 elif element_name[i] == 'Ce': dummy_index[i] = 58 elif element_name[i] == 'Pr': dummy_index[i] = 59 elif element_name[i] == 'Nd': dummy_index[i] = 60 elif element_name[i] == 'Pm': dummy_index[i] = 61 elif element_name[i] == 'Sm': dummy_index[i] = 62 elif element_name[i] == 'Eu': dummy_index[i] = 63 elif element_name[i] == 'Gd': dummy_index[i] = 64 elif element_name[i] == 'Tb': dummy_index[i] = 65 elif element_name[i] == 'Dy': dummy_index[i] = 66 elif element_name[i] == 'Ho': dummy_index[i] = 67 elif element_name[i] == 'Er': dummy_index[i] = 68 elif element_name[i] == 'Tm': dummy_index[i] = 69 elif element_name[i] == 'Yb': dummy_index[i] = 70 elif element_name[i] == 'Lu': dummy_index[i] = 71 elif element_name[i] == 'Hf': dummy_index[i] = 72 elif element_name[i] == 'Ta': dummy_index[i] = 73 elif element_name[i] == 'W': dummy_index[i] = 74 elif element_name[i] == 'Re': dummy_index[i] = 75 elif element_name[i] == 'Os': dummy_index[i] = 76 elif element_name[i] == 'Ir': dummy_index[i] = 77 elif element_name[i] == 'Pt': dummy_index[i] = 78 elif element_name[i] == 'Au': dummy_index[i] = 79 elif element_name[i] == 'Hg': dummy_index[i] = 80 elif element_name[i] == 'Tl': dummy_index[i] = 81 elif element_name[i] == 'Pb': dummy_index[i] = 82 elif element_name[i] == 'Bi': dummy_index[i] = 83 elif element_name[i] == 'Po': dummy_index[i] = 84 elif element_name[i] == 'At': dummy_index[i] = 85 #if spe[0] == 'N 1': # znum_int[0] = 0 # here the index to connect name and atomic numbers. global index_z_for_elements index_z_for_elements = {} for a,b in zip(element_name,dummy_index): index_z_for_elements[a]=b # subroutine to handle grain data, from presolar grain database # files available are specified in the beginning # routine can convert between ratios and deltas, also take inverse of a fraction # check if wanted values are available # grain data files live in validation/grain_data folder # maintained by RT # subroutine: october 2012, RT def graindata_handler(isosx, isosy=None, graintype_in='all', deltax=True, deltay=True, iniabufile_in='../../frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn'): ''' This file reads in grain data and returns them. Possible values for input are: Parameters ---------- isosx: list Isotopes to handle, format ['Si',29,'Si',28]. isosy: list, optional Same format as isox, but not requred. If given, it's read at the same time and given back. required for plotting (2D plotting). The default is None. graintype_in : string, If graintype_in is 'all', then all grains are returned. Also graintype_in could be one of 'sic', 'oxides', 'silicates', 'graphites' or 'misc' for Silicon carbide, Oxides, Silicates, Graphites, and miscellaneous (at the moment Si3N4 grains), respectively. Or - a list combining several grains ['sic','oxides']. - Specific subgroups, e.g., [['sic','M'],['oxides','1','2']] would return mainstream SiC grains as well as group 1 and group 2 oxides. Subgroups are: - SiC: - M >> Mainstream - X >> X-grains - Y >> Y-grains - Z >> Z-grains - AB >> AB grains - N >> nova grains - U >> unclassified grains - graphites: - LD >> low density graphites - HD >> high density graphites - Oxides - 1, 2, 3, or 4 >> respective group of grains - U >> unknown - Silicates - 1, 2, 3, or 4 >> respective group of grains - U >> unknown - Misc - Si3N4 >> silicon nitride grains - IDP - Interplanetary dust particles -> Pepin data only at the moment (no official file) -> Nova - Private - private.txt -> database file in validation folder is used The default is 'all'. deltax : boolean, optional Do you want delta values or not (on x-axis). X axis is the axis of choice if no second isotope is chosen. The default is True. deltay : boolean, optional Same as deltax, but for second axis if chosen. The default is True. iniabufile_in : string File with initial abundances. The default is '../../frames/mppnp/USEEPP/iniab2.0E-02GN93.ppn'. ''' # compatibility if isosx != None and len(isosx) == 2: dumb = [] dumb = isosx[0].split('-') dumb.append(isosx[1].split('-')[0]) dumb.append(isosx[1].split('-')[1]) isosx = dumb if isosy != None and len(isosy) == 2: dumb = [] dumb = isosy[0].split('-') dumb.append(isosy[1].split('-')[0]) dumb.append(isosy[1].split('-')[1]) isosy = dumb # do path scriptpathtmp = __file__ if len(scriptpathtmp.split('/')) == 1: # in folder where nugridse is scriptpathtmp = os.path.abspath('.') + '/utils.py' # to get the current dir svnpathtmp = '/' for i in range(len(scriptpathtmp.split('/'))-3): # -3 to go to folders up! if scriptpathtmp.split('/')[i] != '': svnpathtmp += scriptpathtmp.split('/')[i] + '/' validationpath = svnpathtmp + 'validation/grain_data/' # filenames file_sic = validationpath + 'SiC-All.txt' file_graphite = validationpath + 'graphite-All.txt' file_oxsi = validationpath + 'oxide-silicate-all.txt' file_misc = validationpath + 'miscellaneous-SiN.txt' file_idp = validationpath + 'idp.txt' file_private = validationpath + 'private.txt' # private database file - or compilation of grains (manual filtering...) # process input allgraintypes = [['SiC','M','X','Y','Z','AB','N','U'], ['Oxides','1','2','3','4','U'], ['Silicates','1','2','3','4','U'], ['Graphites','LD','HD','U'], ['Misc','Si3N4'], ['IDP', 'IDP'], ['private','M','C','U/C','X','Y','Z','AB','N','U','1','2','3','4','LD','HD','Si3N4']] if type(graintype_in) == str: if graintype_in.lower() == 'all': graintype = allgraintypes[0:5] # exclude private else: graintype = list() for i in range(len(allgraintypes)): if allgraintypes[i][0].lower() == graintype_in.lower(): graintype.append(allgraintypes[i]) elif type(graintype_in) == list and type(graintype_in[0]) == str: graintype = list() for i in range(len(graintype_in)): for j in range(len(allgraintypes)): if graintype_in[i].lower() == allgraintypes[j][0].lower(): graintype.append(allgraintypes[j]) else: # okay, it's a list in a list, as we want it graintype = list() for i in range(len(graintype_in)): if graintype_in[i][1].lower() != 'all': graintype.append(graintype_in[i]) else: for j in range(len(allgraintypes)): if graintype_in[i][0].lower() == allgraintypes[j][0].lower(): graintype.append(allgraintypes[j]) break # make everything string, if not already for i in range(len(graintype)): for j in range(len(graintype[i])): graintype[i][j] = str(graintype[i][j]) # read in data for the different grains, each line in list is one grain type graindata = list() graindataerr = [] if isosy != None: graindatay = list() graindatayerr = [] graintype_list = list() # as graindata for grain_i in range(len(graintype)): # file if graintype[grain_i][0].lower() == 'sic': fname = file_sic elif graintype[grain_i][0].lower() == 'oxides' or graintype[grain_i][0] == 'silicates': fname = file_oxsi elif graintype[grain_i][0].lower() == 'graphites': fname = file_graphite elif graintype[grain_i][0].lower() == 'misc': fname = file_misc elif graintype[grain_i][0].lower() == 'idp': fname = file_idp elif graintype[grain_i][0].lower() == 'private': fname = file_private else: print 'Problem w/ database filename' return -1 grainclass = graintype[grain_i][0] graintype_tmp = graintype[grain_i][1:len(graintype[grain_i])] # read data if isosy == None: graindata_tmp,graindata_tmperr = _graindata_reader(isosx,isosy,grainclass,graintype_tmp,deltax,deltay,fname,iniabufile_in) if graindata_tmp != -1: # in case header does not exist, or other problem for i in range(len(graindata_tmp)): graintype_list.append([graintype[grain_i][0], graintype_tmp[i]]) graindata.append(graindata_tmp[i]) graindataerr.append(graindata_tmperr[i]) else: graindata_tmp,graindata_tmperr,graindatay_tmp, graindatay_tmperr = _graindata_reader(isosx,isosy,grainclass,graintype_tmp,deltax,deltay,fname,iniabufile_in) if graindata_tmp != -1 and graindatay_tmp != -1: # in case header does not exist, or other problem for i in range(len(graindata_tmp)): if graindata_tmp[i] != []: # empty list graintype_list.append([graintype[grain_i][0], graintype_tmp[i]]) graindata.append(graindata_tmp[i]) graindataerr.append(graindata_tmperr[i]) graindatay.append(graindatay_tmp[i]) graindatayerr.append(graindatay_tmperr[i]) # check if no data at all # give back the read data and the labels for the data if isosy == None: return graintype_list, graindata, graindataerr else: return graintype_list, graindata, graindataerr, graindatay, graindatayerr ### private routine to read data ### def _graindata_reader(isos, isos2, gclass, gtype, deltax, deltay, fname, iniabufile): ''' Private function to actually read grain files. Only called from graindata_handler! Parameters ---------- isos : list Isotopes, style: ['Fe',58,'Fe',56]. isos2 : list If second isotope is required. gclass : string Grain class, e.g., 'sic' gtype : list Grain type -> list! ''' # isotopes iso1 = str(isos[1]) + isos[0] iso2 = str(isos[3]) + isos[2] if isos2 != None: iso3 = str(isos2[1]) + isos2[0] iso4 = str(isos2[3]) + isos2[2] # read file f = open(fname,'r') f_in = f.read().split('\n') f.close() # data header = f_in[0].split('\t') data = list() for i in range(1,len(f_in)): datatmp = list() for j in range(len(f_in[i].split('\t'))): datatmp.append(f_in[i].split('\t')[j]) data.append(datatmp) # remove " from header header2 = header header = [] for i in header2: header.append(i.replace('"','')) # check if value is availabe as ratio or as delta, normal or 1 / (ratio) for i in range(len(header)): if header[i].lower() == iso1.lower() + '/' + iso2.lower(): # ratio index = i one_over = False deltadb = False break elif header[i].lower() == iso2.lower() + '/' + iso1.lower(): # 1 / ratio index = i one_over = True deltadb = False break elif header[i].lower() == 'd(' + iso1.lower() + '/' + iso2.lower() + ')': # delta ratio index = i one_over = False deltadb = True break elif header[i].lower() == 'd(' + iso2.lower() + '/' + iso1.lower() + ')': # delta 1 / ratio index = i one_over = True deltadb = True break else: index = -1 if index == -1: print 'Data entry does not exist for ' + iso1 + '/' + iso2 if isos2 == None: return -1,-1 # can be read by other file w/o error! else: return -1,-1,-1,-1 if isos2 != None: for i in range(len(header)): if header[i].lower() == iso3.lower() + '/' + iso4.lower(): # ratio index2 = i one_over2 = False deltadb2 = False break elif header[i].lower() == iso4.lower() + '/' + iso3.lower(): # 1 / ratio index2 = i one_over2 = True deltadb2 = False break elif header[i].lower() == 'd(' + iso3.lower() + '/' + iso4.lower() + ')': # 1 / ratio index2 = i one_over2 = False deltadb2 = True break elif header[i].lower() == 'd(' + iso4.lower() + '/' + iso3.lower() + ')': # 1 / ratio index2 = i one_over2 = True deltadb2 = True break else: index2 = -1 if index2 == -1: print 'Data entry does not exist for ' + iso3 + '/' + iso4 return -1,-1,-1,-1 # can be read by other file w/o error! # get solar system ratio of wanted isotopes inut = iniabu(iniabufile) ss_ratio = inut.isoratio_init(isos) if isos2 != None: ss_ratio2 = inut.isoratio_init(isos2) # now make return list dataret = list() datareterr = [] if isos2 != None: dataret2 = list() dataret2err = [] for gtype_i in gtype: datatmp = list() # one list for each graintype datatmperr = [] if isos2 != None: datatmp2 = list() datatmp2err = [] for i in range(len(data)): if data[i][0].replace(' ','') != '': # oxide and silicate checker if gclass == 'oxides': if data[i][1].lower() == 'oxide': ifchecksilox = True else: ifchecksilox = False elif gclass == 'silicates': if data[i][1].lower() == 'silicate': ifchecksilox = True else: ifchecksilox = False else: ifchecksilox = True if data[i][0].lower() == gtype_i.lower() and ifchecksilox: # 0 is the column where the grain type is saved, also silicates and oxides separation if isos2 == None: ifcheck = data[i][index].replace(' ','') != '' # replace to remove spaces from empty string w/ just spaces else: ifcheck = data[i][index].replace(' ','') != '' and data[i][index2].replace(' ','') != '' if ifcheck: try: cellvalue = float(data[i][index]) # in case of database error, give it back cellvalueerr = data[i][index+1].replace(' ','') if cellvalueerr == '': cellvalueerr = 0. else: cellvalueerr = float(data[i][index+1])/cellvalue # relative uncertainty # run through if loop if deltax and deltadb and one_over != True: # delta requested, delta available, this ratio tmpvalue = cellvalue datatmp.append(cellvalue) elif deltax and deltadb and one_over: # delta requested, delta available, one over ratio tmpvalue2 = (cellvalue/1000. + 1) * (1. / ss_ratio) # now ratio tmpvalue = ((1/tmpvalue2) / ss_ratio - 1.) * 1000. datatmp.append(tmpvalue) elif deltax != True and deltadb != True and one_over != True: # ratio requested, ratio available, this ratio tmpvalue = cellvalue datatmp.append(cellvalue) elif deltax != True and deltadb != True and one_over: # ratio requested, ratio available, one over ratio tmpvalue = 1. / cellvalue datatmp.append(tmpvalue) elif deltax and deltadb != True and one_over != True: # delta requested, ratio available, this ratio tmpvalue = (cellvalue / ss_ratio - 1.) * 1000. datatmp.append(tmpvalue) elif deltax and deltadb != True and one_over: # delta requested, ratio available, one over ratio tmpvalue = ((1/cellvalue) / ss_ratio - 1.) * 1000. datatmp.append(tmpvalue) elif deltax != True and deltadb and one_over != True: # ratio requested, delta available, this ratio tmpvalue = (cellvalue / 1000. + 1.) * ss_ratio datatmp.append(tmpvalue) elif deltax != True and deltadb and one_over: # ratio requested, delta available, one over ratio tmpvalue = ((1 / cellvalue) / 1000. + 1.) * ss_ratio datatmp.append(tmpvalue) tmpvalueerr = cellvalueerr * tmpvalue datatmperr.append(tmpvalueerr) except ValueError: print 'error in database, file: ' + fname + ' line ' + str(i+2) + ' column ' + str(index+1) print data[i][index] print data[i][index+1] if isos2 != None: if ifcheck: try: cellvalue = float(data[i][index2]) # in case of database error, give it back cellvalueerr = data[i][index2+1].replace(' ','') if cellvalueerr == '': cellvalueerr = 0. else: cellvalueerr = float(data[i][index2+1])/cellvalue # relative uncertainty # run trough if loop if deltay and deltadb2 and one_over2 != True: # delta requested, delta available, this ratio tmpvalue = cellvalue datatmp2.append(cellvalue) elif deltay and deltadb2 and one_over2: # delta requested, delta available, one over ratio tmpvalue2 = (cellvalue/1000. + 1) * (1. / ss_ratio) # now ratio tmpvalue = ((1/tmpvalue2) / ss_ratio - 1.) * 1000. datatmp2.append(tmpvalue) elif deltay != True and deltadb2 != True and one_over2 != True: # ratio requested, ratio available, this ratio tmpvalue = cellvalue datatmp2.append(cellvalue) elif deltay != True and deltadb2 != True and one_over2: # ratio requested, ratio available, one over ratio tmpvalue = 1. / cellvalue datatmp2.append(tmpvalue) elif deltay and deltadb2 != True and one_over2 != True: # delta requested, ratio available, this ratio tmpvalue = (cellvalue / ss_ratio2 - 1.) * 1000. datatmp2.append(tmpvalue) elif deltay and deltadb2 != True and one_over2: # delta requested, ratio available, one over ratio tmpvalue = ((1/cellvalue) / ss_ratio2 - 1.) * 1000. datatmp2.append(tmpvalue) elif deltay != True and deltadb2 and one_over2 != True: # ratio requested, delta available, this ratio tmpvalue = (cellvalue / 1000. + 1.) * ss_ratio2 datatmp2.append(tmpvalue) elif deltay != True and deltadb2 and one_over2: # ratio requested, delta available, one over ratio tmpvalue = ((1 / cellvalue) / 1000. + 1.) * ss_ratio2 datatmp2.append(tmpvalue) tmpvalueerr = cellvalueerr * tmpvalue datatmp2err.append(tmpvalueerr) except ValueError: print 'error in database, file: ' + fname + ' line ' + str(i+2) + ' column ' + str(index2+1) dataret.append(datatmp) datareterr.append(datatmperr) if isos2 != None: dataret2.append(datatmp2) dataret2err.append(datatmp2err) # return the data if isos2 == None: return dataret,datareterr else: return dataret,datareterr,dataret2,dataret2err def get_z_from_el(element): ''' Very simple function that gives the atomic number AS A STRING when given the element symbol. Uses predefined a dictionnary. Parameter : element : string For the other way, see get_el_from_z ''' dict_name={'Ru': '44', 'Re': '75', 'Ra': '88', 'Rb': '37', 'Rn': '86', 'Rh': '45', 'Be': '4', 'Ba': '56', 'Bi': '83', 'Br': '35', 'H': '1', 'P': '15', 'Os': '76', 'Hg': '80', 'Ge': '32', 'Gd': '64', 'Ga': '31', 'Pr': '59', 'Pt': '78', 'C': '6', 'Pb': '82', 'Pa': '91', 'Pd': '46', 'Cd': '48', 'Po': '84', 'Pm': '61', 'Ho': '67', 'Hf': '72', 'K': '19', 'He': '2', 'Mg': '12', 'Mo': '42', 'Mn': '25', 'O': '8', 'S': '16', 'W': '74', 'Zn': '30', 'Eu': '63', 'Zr': '40', 'Er': '68', 'Ni': '28', 'Na': '11', 'Nb': '41', 'Nd': '60', 'Ne': '10', 'Fr': '87', 'Fe': '26', 'B': '5', 'F': '9', 'Sr': '38', 'N': '7', 'Kr': '36', 'Si': '14', 'Sn': '50', 'Sm': '62', 'V': '23', 'Sc': '21', 'Sb': '51', 'Se': '34', 'Co': '27', 'Cl': '17', 'Ca': '20', 'Ce': '58', 'Xe': '54', 'Lu': '71', 'Cs': '55', 'Cr': '24', 'Cu': '29', 'La': '57', 'Li': '3', 'Tl': '81', 'Tm': '69', 'Th': '90', 'Ti': '22', 'Te': '52', 'Tb': '65', 'Tc': '43', 'Ta': '73', 'Yb': '70', 'Dy': '66', 'I': '53', 'U': '92', 'Y': '39', 'Ac': '89', 'Ag': '47', 'Ir': '77', 'Al': '13', 'As': '33', 'Ar': '18', 'Au': '79', 'At': '85', 'In': '49'} return int(dict_name[element]) def get_el_from_z(z): ''' Very simple Vfunction that gives the atomic number AS A STRING when given the element symbol. Uses predefined a dictionnary. Parameter : z : string or number For the other way, see get_z_from_el ''' if(type(z)==float): z=int(z) if(type(z)==int): z=str(z) dict_z={'24': 'Cr', '25': 'Mn', '26': 'Fe', '27': 'Co', '20': 'Ca', '21': 'Sc', '22': 'Ti', '23': 'V', '28': 'Ni', '29': 'Cu', '4': 'Be', '8': 'O', '59': 'Pr', '58': 'Ce', '55': 'Cs', '54': 'Xe', '57': 'La', '56': 'Ba', '51': 'Sb', '50': 'Sn', '53': 'I', '52': 'Te', '88': 'Ra', '89': 'Ac', '82': 'Pb', '83': 'Bi', '80': 'Hg', '81': 'Tl', '86': 'Rn', '87': 'Fr', '84': 'Po', '85': 'At', '3': 'Li', '7': 'N', '39': 'Y', '38': 'Sr', '33': 'As', '32': 'Ge', '31': 'Ga', '30': 'Zn', '37': 'Rb', '36': 'Kr', '35': 'Br', '34': 'Se', '60': 'Nd', '61': 'Pm', '62': 'Sm', '63': 'Eu', '64': 'Gd', '65': 'Tb', '66': 'Dy', '67': 'Ho', '68': 'Er', '69': 'Tm', '2': 'He', '6': 'C', '91': 'Pa', '90': 'Th', '92': 'U', '11': 'Na', '10': 'Ne', '13': 'Al', '12': 'Mg', '15': 'P', '14': 'Si', '17': 'Cl', '16': 'S', '19': 'K', '18': 'Ar', '48': 'Cd', '49': 'In', '46': 'Pd', '47': 'Ag', '44': 'Ru', '45': 'Rh', '42': 'Mo', '43': 'Tc', '40': 'Zr', '41': 'Nb', '1': 'H', '5': 'B', '9': 'F', '77': 'Ir', '76': 'Os', '75': 'Re', '74': 'W', '73': 'Ta', '72': 'Hf', '71': 'Lu', '70': 'Yb', '79': 'Au', '78': 'Pt'} return dict_z[z]