diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 6f6ae99b..bcb4302e 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -23,8 +23,11 @@ jobs: python -m pip install --upgrade pip pip install -r requirements.txt - - name: Run tests + - name: Run unit tests run: python -m pytest -vv + - name: Run integration tests + run: ./tests/integration_tests.sh + - name: Run linter run: python -m black neucbot/ tests/ --check diff --git a/neucbot.py b/neucbot.py index 26ae01a6..360b5c39 100755 --- a/neucbot.py +++ b/neucbot.py @@ -7,13 +7,13 @@ import subprocess import shutil import math -from neucbot import ensdf -from neucbot import elements from neucbot import alpha +from neucbot import ensdf +from neucbot import elements +from neucbot import material class constants: - N_A = 6.0221409e+23 MeV_to_keV= 1.e3 mb_to_cm2 = 1.e-27 year_to_s = 31536000 @@ -28,12 +28,6 @@ class constants: force_recalculation = False ofile = sys.stdout -class material: - def __init__(self,e,a,f): - self.ele = str(e) - self.A = float(a) - self.frac = float(f) - def isoDir(ele,A): return './Data/Isotopes/'+ele.capitalize()+'/'+ele.capitalize()+str(int(A))+'/' @@ -68,100 +62,6 @@ def loadChainAlphaList(fname): alpha_list.append([ene, intensity*br/100]) return alpha_list -def readTargetMaterial(fname): - f = open(fname) - mat_comp = [] - tokens = [line.split() for line in f.readlines()] - for line in tokens: - if len(list(line)) < 3: - continue - if line[0][0] == '#': - continue - element = elements.Element(line[0]) - A = int(line[1]) - frac = float(line[2]) - - if A == 0: - natIso_list = element.isotopes() - print(natIso_list) - for A_i in natIso_list: - abund = element.abundance(A_i) - print(A_i) - print(abund) - mater = material(element.symbol,A_i,frac*abund/100.) - mat_comp.append(mater) - else: - mater = material(element.symbol,A,frac) - mat_comp.append(mater) - - # Normalize - norm = 0 - for mat in mat_comp: - norm += mat.frac - for mat in mat_comp: - mat.frac /= norm - - return mat_comp - -def calcStoppingPower(e_alpha_MeV,mat_comp): - # Stopping power as units of keV/(mg/cm^2) or MeV/(g/cm^2) - e_alpha = e_alpha_MeV - sp_total = 0 - # First, reduce the material to combine all isotopes with the same Z - mat_comp_reduced = {} - for mat in mat_comp: - if mat.ele in mat_comp_reduced: - mat_comp_reduced[mat.ele] += mat.frac - else: - mat_comp_reduced[mat.ele] = mat.frac - - # Then, for each element, get the stopping power at this alpha energy - for mat in mat_comp_reduced: - spDir = './Data/StoppingPowers/' - spFile = spDir + mat.lower() + '.dat' - spf = open(spFile) - - tokens = [line.split() for line in spf.readlines()] - first = True - sp_found = False - e_curr = 0 - e_last = 0 - sp_curr = 0 - sp_last = 0 - sp_alpha = 0 - for line in tokens: - if line[0][0] == '#': - continue - e_curr = float(line[0]) - if str(line[1]) == 'keV': - e_curr /= 1000 - elif str(line[1]) == 'MeV': - e_curr *= 1 - sp_curr = float(line[3])+float(line[2]) - - # Alpha energy is below the list. Use the lowest energy in the list - if e_curr > e_alpha and first: - first = False - sp_found = True - sp_alpha = sp_curr - break - # If this entry is above the alpha energy, the alpha is between this - # entry and the previous one - if e_curr > e_alpha: - first = False - sp_alpha = (sp_curr-sp_last)*(e_alpha-e_last)/(e_curr-e_last) + sp_last - sp_found = True - break - # Otherwise, keep looking for the entry - first = False - sp_last = sp_curr - e_last = e_curr - # if the alpha energy is too high for the list, use the highest energy on the list - if not sp_found: - sp_alpha = sp_last - sp_total += sp_alpha * mat_comp_reduced[mat]/100 - return sp_total - def runTALYS(e_a, ele, A): iso = str(ele)+str(int(A)) inpdir = isoDir(ele,A) + 'TalysInputs/' @@ -208,14 +108,6 @@ def runTALYS(e_a, ele, A): blank_f.write("EMPTY") blank_f.close() - -def getMatTerm(mat,mat_comp): - # mat_comp structure: [ele,A,frac] - A = mat.A - conc = mat.frac/100. - mat_term = (constants.N_A * conc)/A - return mat_term - def getIsotopeDifferentialNSpec(e_a, ele, A): target = ele+str(int(A)) path = isoDir(ele,A) + 'NSpectra/' @@ -382,20 +274,20 @@ def run_alpha(alpha_list, mat_comp, e_alpha_step): stopping_power = 0 if stopping_power == 0: - stopping_power = calcStoppingPower(e_a, mat_comp) - for mat in mat_comp: - mat_term = getMatTerm(mat,mat_comp) + stopping_power = mat_comp.stopping_power(e_a) + for mat in mat_comp.materials: + mat_term = mat.material_term() # Get alpha n spectrum for this alpha and this target - spec_raw = getIsotopeDifferentialNSpec(e_a, mat.ele, mat.A) + spec_raw = getIsotopeDifferentialNSpec(e_a, mat.element.symbol, mat.mass_number) spec = rebin(spec_raw,constants.delta_bin,constants.min_bin,constants.max_bin) # Add this spectrum to the total spectrum delta_ea = e_alpha_step if e_a - e_alpha_step < 0: delta_ea = e_a prefactors = (intensity/100.)*mat_term*delta_ea/stopping_power - xsect = prefactors * readTotalNXsect(e_a,mat.ele,mat.A) + xsect = prefactors * readTotalNXsect(e_a,mat.element.symbol, mat.mass_number) total_xsect += xsect - matname = str(mat.ele)+str(mat.A) + matname = str(mat.element.symbol)+str(mat.mass_number) if matname in xsects: xsects[matname] += xsect else: @@ -441,7 +333,7 @@ def main(): if arg == '-m': mat_file = sys.argv[sys.argv.index(arg)+1] print('load target material', mat_file, file = sys.stdout) - mat_comp = readTargetMaterial(mat_file) + mat_comp = material.Composition.from_file(mat_file) if arg == '-s': alpha_step_size = float(sys.argv[sys.argv.index(arg)+1]) print('step size', alpha_step_size, file = sys.stdout) @@ -473,9 +365,9 @@ def main(): constants.ofile = open(ofile,'w') #sys.stdout = open(ofile,'w') - if len(alpha_list) == 0 or len(mat_comp) == 0: + if len(alpha_list) == 0 or mat_comp.empty(): if len(alpha_list)==0: print('No alpha list or chain specified', file = sys.stdout) - if len(mat_comp)==0: print('No target material specified', file = sys.stdout) + if mat_comp.empty(): print('No target material specified', file = sys.stdout) print('', file = sys.stdout) help_message() return 0 @@ -488,8 +380,8 @@ def main(): print(alph[0],'&', alph[1],'\\\\', file = sys.stdout) if constants.download_data: - for mat in mat_comp: - ele = mat.ele + for mat in mat_comp.materials: + ele = mat.element.symbol if not os.path.exists('./Data/Isotopes/'+ele.capitalize()): if constants.download_version == 2: print('\tDownloading (datset V2) data for',ele, file = sys.stdout) diff --git a/neucbot/elements.py b/neucbot/elements.py index 6aed9303..9421998a 100644 --- a/neucbot/elements.py +++ b/neucbot/elements.py @@ -39,5 +39,6 @@ def __init__(self, element): def isotopes(self): return [iso for iso in self.isos.keys()] + # This should return a value between 0 and 1 def abundance(self, isotope): - return float(self.isos.get(isotope).get("abundance")) + return float(self.isos.get(isotope).get("abundance")) / 100.0 diff --git a/neucbot/material.py b/neucbot/material.py new file mode 100644 index 00000000..8e2f0fb7 --- /dev/null +++ b/neucbot/material.py @@ -0,0 +1,152 @@ +from bisect import bisect + +from neucbot import elements + +N_A = 6.0221409e23 + + +class Isotope: + # self.fraction should be a number between 0 and 1 + def __init__(self, element, mass_number, fraction): + self.element = element + self.mass_number = int(mass_number) + self.fraction = float(fraction) + + def material_term(self): + return (N_A * self.fraction) / self.mass_number + + +class StoppingPowerList: + def __init__(self, element_symbol): + self.element_symbol = element_symbol + self.stopping_powers = {} + + def load_file(self): + file_path = f"./Data/StoppingPowers/{self.element_symbol.lower()}.dat" + file = open(file_path) + + for data in [ + line.split() for line in file.readlines() if not line.startswith("#") + ]: + energy = float(data[0]) + units = str(data[1]) + stopping_power = float(data[2]) + float(data[3]) + + if units == "keV": + energy /= 1000 + + self.stopping_powers[energy] = stopping_power + + # Use binary search to find energy in log(N) time + def for_alpha(self, alpha_energy): + energy_intervals = list(self.stopping_powers.keys()) + min_energy = energy_intervals[0] + max_energy = energy_intervals[-1] + + if alpha_energy < min_energy: + return self.stopping_powers[min_energy] + elif alpha_energy > max_energy: + return self.stopping_powers[max_energy] + + range_end = bisect(energy_intervals, alpha_energy) + range_start = range_end - 1 + + energy_start = energy_intervals[range_start] + energy_end = energy_intervals[range_end] + energy_diff = (alpha_energy - energy_start) / (energy_end - energy_start) + + stop_power_start = self.stopping_powers[energy_start] + stop_power_end = self.stopping_powers[energy_end] + + return (stop_power_end - stop_power_start) * energy_diff + stop_power_start + + +class Composition: + @classmethod + def from_file(cls, file_path): + file = open(file_path) + + composition = cls() + + for material in [ + line.split() for line in file.readlines() if not line.startswith("#") + ]: + if len(material) < 3: + continue + + element = elements.Element(material[0]) + mass_number = int(material[1]) + fraction = float(material[2]) + + # If a single mass number isn't specified, use all isotopes + # along with their natural abundances + if mass_number == 0: + for isotope in element.isotopes(): + composition.add( + Isotope( + element, + isotope, + fraction * element.abundance(isotope) / 100.0, + ) + ) + + # Otherwise, if a single mass number is provided, + # use the fraction provided + else: + composition.add( + Isotope( + element, + mass_number, + fraction / 100.0, + ) + ) + + composition.normalize() + composition.populate_stopping_powers() + + return composition + + def __init__(self): + self.materials = [] + self.fractions = {} + self.stopping_powers = {} + + def normalize(self): + norm = 0 + + for material in self.materials: + norm += material.fraction + + for material in self.materials: + material.fraction /= norm + + # Computes the fraction of this element in the overall material, + # grouping isotopes with the same Z + symbol = material.element.symbol + self.fractions[symbol] = self.fractions.get(symbol, 0) + material.fraction + + def populate_stopping_powers(self): + # Populates stopping powers from ./Data/StoppingPowers once so that + # they don't need to be populated for every energy step + for element in self.fractions: + stop_power_list = StoppingPowerList(element) + stop_power_list.load_file() + + self.stopping_powers[element] = stop_power_list + + def empty(self): + return len(self.materials) == 0 + + def add(self, material): + self.materials.append(material) + + # Expects an alpha energy in units of MeV + def stopping_power(self, e_alpha): + total_stopping_power = 0 + + for element, fraction in self.fractions.items(): + element_stop_power = self.stopping_powers[element].for_alpha(e_alpha) + + total_stopping_power += element_stop_power * fraction + + return total_stopping_power diff --git a/tests/integration_tests.sh b/tests/integration_tests.sh new file mode 100755 index 00000000..71d6f0d9 --- /dev/null +++ b/tests/integration_tests.sh @@ -0,0 +1,14 @@ +#!/bin/bash +set -e + +python3 ./neucbot.py -m Materials/Acrylic.dat -c Chains/Th232Chain.dat -d v2 -o tmp-acrylic-th232-chain.txt +diff tmp-acrylic-th232-chain.txt tests/integration_tests/acrylic-th232-chain.txt +rm tmp-acrylic-th232-chain.txt + +python3 ./neucbot.py -m Materials/Acrylic.dat -l AlphaLists/Rn220Alphas.dat -d v2 -o tmp-acrylic-rn220-alphalist.txt +diff tmp-acrylic-rn220-alphalist.txt tests/integration_tests/acrylic-rn220-alphalist.txt +rm tmp-acrylic-rn220-alphalist.txt + +python3 ./neucbot.py -m Materials/Acrylic.dat -l AlphaLists/Bi212Alphas.dat -d v2 -o tmp-acrylic-bi212-alphalist.txt +diff tmp-acrylic-bi212-alphalist.txt tests/integration_tests/acrylic-bi212-alphalist.txt +rm tmp-acrylic-bi212-alphalist.txt diff --git a/tests/integration_tests/acrylic-bi212-alphalist.txt b/tests/integration_tests/acrylic-bi212-alphalist.txt new file mode 100644 index 00000000..83cd1551 --- /dev/null +++ b/tests/integration_tests/acrylic-bi212-alphalist.txt @@ -0,0 +1,160 @@ + +# Total neutron yield = 1.865130270096359e-07 n/decay + C12 0.0 + C13 1.6378752779824458e-07 + H1 0.0 + H2 0.0 + O16 0.0 + O17 2.3211544641210924e-09 + O18 2.0404344747270445e-08 +# Integral of spectrum = 1.9061304301735257e-07 n/decay +100 5.330112308624654e-11 +200 3.9376208188415496e-11 +300 4.001384730671926e-11 +400 4.781442246128734e-11 +500 5.4826660748663663e-11 +600 6.20450363652541e-11 +700 5.975786053045866e-11 +800 8.511980524683539e-12 +900 4.452296612431278e-12 +1000 4.98009158798087e-12 +1100 5.500154431958868e-12 +1200 6.026645257835461e-12 +1300 5.346850388458856e-12 +1400 3.440267261296956e-12 +1500 3.260934424284713e-12 +1600 3.796896868019909e-12 +1700 4.374509239229058e-12 +1800 4.987107784149935e-12 +1900 5.629368473586987e-12 +2000 6.281517981547197e-12 +2100 6.914639101858719e-12 +2200 7.529871436310418e-12 +2300 8.066211860147639e-12 +2400 7.806738498734402e-12 +2500 6.870577790003208e-12 +2600 6.909735908227458e-12 +2700 7.413018084995017e-12 +2800 7.955659030975638e-12 +2900 8.518056188995167e-12 +3000 9.118796030490186e-12 +3100 9.779083826652588e-12 +3200 1.0514187890563138e-11 +3300 1.1328626721938092e-11 +3400 1.2219237040734537e-11 +3500 1.317096673318369e-11 +3600 1.4097596413606278e-11 +3700 1.4780840208264885e-11 +3800 1.5077794163312583e-11 +3900 1.527433917119779e-11 +4000 1.58266200925234e-11 +4100 1.6805129023908367e-11 +4200 1.807594153495321e-11 +4300 1.970133694487426e-11 +4400 2.179906625729911e-11 +4500 2.4308943955656753e-11 +4600 2.706539426520013e-11 +4700 2.993863566478047e-11 +4800 3.285735060947402e-11 +4900 3.577702848373924e-11 +5000 3.866067359601289e-11 +5100 4.1473943368956715e-11 +5200 4.418478773381399e-11 +5300 4.67637611254364e-11 +5400 4.918423995649324e-11 +5500 5.142185594627517e-11 +5600 5.345168292594853e-11 +5700 5.524357516838596e-11 +5800 5.675633593376265e-11 +5900 5.793153668537925e-11 +6000 5.868702385216847e-11 +6100 5.891019885085282e-11 +6200 5.845450472005078e-11 +6300 5.714653421176892e-11 +6400 5.481177558284542e-11 +6500 5.132184497391116e-11 +6600 4.665407021853484e-11 +6700 4.0942822011253944e-11 +6800 3.4496253537047574e-11 +6900 2.776167769142191e-11 +7000 2.124231167988514e-11 +7100 1.539212634547179e-11 +7200 1.0525605224514981e-11 +7300 6.773140612309679e-12 +7400 4.091467805052524e-12 +7500 2.315467892060277e-12 +7600 1.2255803499407363e-12 +7700 6.058758118350932e-13 +7800 2.7942109380778477e-13 +7900 1.201030537879883e-13 +8000 4.807532501571943e-14 +8100 1.7909837360464826e-14 +8200 6.2071257913824475e-15 +8300 2.00129572348099e-15 +8400 6.009226550880638e-16 +8500 1.6875068446204716e-16 +8600 4.4977214486755855e-17 +8700 1.1942163399994318e-17 +8800 3.607736573157667e-18 +8900 1.5221698212304653e-18 +9000 9.27396022247372e-19 +9100 6.856386168727167e-19 +9200 5.404398007419637e-19 +9300 4.340625453813333e-19 +9400 3.510777624213833e-19 +9500 2.851676853956289e-19 +9600 2.324561123162766e-19 +9700 1.9011250841552484e-19 +9800 1.5596651760234166e-19 +9900 1.2833313598653127e-19 +10000 1.0589381623224982e-19 +10100 8.761334660964267e-20 +10200 7.267498438648091e-20 +10300 6.043126126044168e-20 +10400 5.0367935852633504e-20 +10500 4.2074087358795934e-20 +10600 3.522086761182016e-20 +10700 2.954399136100848e-20 +10800 2.483021916129612e-20 +10900 7.513478780018326e-21 +11000 5.920862684572933e-21 +11100 1.066830831516095e-22 +11200 7.144994383952812e-23 +11300 9.552234608023408e-26 +11400 5.121074335434301e-26 +11500 1.0109459281577453e-29 +11600 4.16074090138763e-30 +11700 2.2060096194836707e-34 +11800 6.762212596291996e-35 +11900 1.091986797032694e-39 +12000 2.413324683226788e-40 +12100 3.609447760960841e-45 +12200 5.7321587757368785e-46 +12300 8.907884187023513e-47 +12400 1.3546178392057362e-47 +12500 2.0156472855371275e-48 +12600 0.0 +12700 0.0 +12800 0.0 +12900 0.0 +13000 0.0 +13100 0.0 +13200 0.0 +13300 0.0 +13400 0.0 +13500 0.0 +13600 0.0 +13700 0.0 +13800 0.0 +13900 0.0 +14000 0.0 +14100 0.0 +14200 0.0 +14300 0.0 +14400 0.0 +14500 0.0 +14600 0.0 +14700 0.0 +14800 0.0 +14900 0.0 +15000 0.0 diff --git a/tests/integration_tests/acrylic-rn220-alphalist.txt b/tests/integration_tests/acrylic-rn220-alphalist.txt new file mode 100644 index 00000000..dafdfe7f --- /dev/null +++ b/tests/integration_tests/acrylic-rn220-alphalist.txt @@ -0,0 +1,165 @@ + +# Total neutron yield = 2.175050799970592e-07 n/decay + C12 0.0 + C13 1.9090874958585992e-07 + H1 0.0 + H2 0.0 + O16 0.0 + O17 2.717884256663278e-09 + O18 2.3878446154536136e-08 +# Integral of spectrum = 2.234829523074957e-07 n/decay +100 8.399058799577613e-11 +200 3.93741047927989e-11 +300 4.0395498090183434e-11 +400 4.9071723453597974e-11 +500 5.699599925241164e-11 +600 6.560463298572045e-11 +700 7.31936685837022e-11 +800 7.951891038182988e-11 +900 6.166312735573753e-11 +1000 7.45012343798205e-12 +1100 5.584268309102103e-12 +1200 6.162889659774142e-12 +1300 6.781056908689634e-12 +1400 7.419698984109826e-12 +1500 6.607347306530609e-12 +1600 4.45358457644802e-12 +1700 4.373822275867504e-12 +1800 4.985630933689642e-12 +1900 5.6279035271055364e-12 +2000 6.285560756322328e-12 +2100 6.94321787409264e-12 +2200 7.592443683666304e-12 +2300 8.245116512407125e-12 +2400 8.897822401124406e-12 +2500 9.431116732064359e-12 +2600 9.061928365471975e-12 +2700 7.975750554343153e-12 +2800 7.985738803242626e-12 +2900 8.515791638669295e-12 +3000 9.116024964440568e-12 +3100 9.77630006431952e-12 +3200 1.05125994509859e-11 +3300 1.1333327687724134e-11 +3400 1.2244739316642701e-11 +3500 1.3259541764784031e-11 +3600 1.4393933073646168e-11 +3700 1.564145105682864e-11 +3800 1.6902617343609524e-11 +3900 1.793661874660029e-11 +4000 1.8590941914546e-11 +4100 1.915767480846843e-11 +4200 2.0091484012344917e-11 +4300 2.1441253304247503e-11 +4400 2.3049653538008326e-11 +4500 2.496752226617796e-11 +4600 2.730533697682226e-11 +4700 2.9993606469046606e-11 +4800 3.2859061160756566e-11 +4900 3.576888231772832e-11 +5000 3.865288989798765e-11 +5100 4.147052277607339e-11 +5200 4.419017104831832e-11 +5300 4.67828826305539e-11 +5400 4.922057351461002e-11 +5500 5.1475810563170764e-11 +5600 5.3522469881803106e-11 +5700 5.5337018080809e-11 +5800 5.689858432084474e-11 +5900 5.81867358772636e-11 +6000 5.917702234354171e-11 +6100 5.9833846425872e-11 +6200 6.009988175980547e-11 +6300 5.9882958303298e-11 +6400 5.90459076183967e-11 +6500 5.7410391995560456e-11 +6600 5.478519061114634e-11 +6700 5.1022338236564706e-11 +6800 4.608856379241022e-11 +6900 4.01254606124186e-11 +7000 3.346644382725309e-11 +7100 2.6592280701837744e-11 +7200 2.0032566873200198e-11 +7300 1.4247227026106653e-11 +7400 9.532410745648528e-12 +7500 5.982374102975212e-12 +7600 3.513043062047288e-12 +7700 1.926449999977861e-12 +7800 9.848577615913349e-13 +7900 4.687443687300258e-13 +8000 2.0746758779389327e-13 +8100 8.531100002742838e-14 +8200 3.256652793266868e-14 +8300 1.1534012123034703e-14 +8400 3.78861721986466e-15 +8500 1.1542706660677161e-15 +8600 3.2665222166086473e-16 +8700 8.636319554945377e-17 +8800 2.178208055494739e-17 +8900 5.6242607094318516e-18 +9000 1.7829328601150386e-18 +9100 8.504138471321703e-19 +9200 5.696646056517925e-19 +9300 4.387660784538719e-19 +9400 3.517045431359482e-19 +9500 2.8518242618837565e-19 +9600 2.323972469575384e-19 +9700 1.9005499374463495e-19 +9800 1.5591819002928417e-19 +9900 1.282932422405795e-19 +10000 1.0586088463990778e-19 +10100 8.758609872535646e-20 +10200 7.265238224255218e-20 +10300 6.041246694288412e-20 +10400 5.035227126110674e-20 +10500 4.206100218102141e-20 +10600 3.5209913807634883e-20 +10700 2.953480308376916e-20 +10800 2.4822496882515984e-20 +10900 7.511142063722308e-21 +11000 5.9190212770536474e-21 +11100 1.0664990436636127e-22 +11200 7.142772267500401e-23 +11300 9.549263832045272e-26 +11400 5.119481664688423e-26 +11500 1.0106315206916594e-29 +11600 4.159446897457836e-30 +11700 2.205323543329344e-34 +11800 6.760109526218397e-35 +11900 1.0916471855932736e-39 +12000 2.4125741314145106e-40 +12100 3.608325210991356e-45 +12200 5.73037605574594e-46 +12300 8.905113806118437e-47 +12400 1.3541965486594494e-47 +12500 2.0150204126867424e-48 +12600 0.0 +12700 0.0 +12800 0.0 +12900 0.0 +13000 0.0 +13100 0.0 +13200 0.0 +13300 0.0 +13400 0.0 +13500 0.0 +13600 0.0 +13700 0.0 +13800 0.0 +13900 0.0 +14000 0.0 +14100 0.0 +14200 0.0 +14300 0.0 +14400 0.0 +14500 0.0 +14600 0.0 +14700 0.0 +14800 0.0 +14900 0.0 +15000 0.0 +15100 0.0 +15200 0.0 +15300 0.0 +15400 0.0 +15500 0.0 diff --git a/tests/integration_tests/acrylic-th232-chain.txt b/tests/integration_tests/acrylic-th232-chain.txt new file mode 100644 index 00000000..31cc20ce --- /dev/null +++ b/tests/integration_tests/acrylic-th232-chain.txt @@ -0,0 +1,181 @@ + +# Total neutron yield = 1.318008946796577e-06 n/decay + C12 0.0 + C13 1.1595483411751048e-06 + H1 0.0 + H2 0.0 + O16 0.0 + O17 1.681566318550402e-08 + O18 1.4164494243596873e-07 +# Integral of spectrum = 1.3469312787505442e-06 n/decay +0 3.849284879079153e-13 +100 3.99263896445411e-10 +200 2.626860489058724e-10 +300 2.122811436291354e-10 +400 2.4330157085886814e-10 +500 2.136018515274356e-10 +600 2.36454856040748e-10 +700 2.610336277179039e-10 +800 2.636230655267057e-10 +900 2.2963767238716684e-10 +1000 2.1383501714644674e-10 +1100 2.2008150122772638e-10 +1200 2.3125371337781055e-10 +1300 1.6016454932718687e-10 +1400 1.4753299515158148e-10 +1500 1.518650494631864e-10 +1600 1.563403734213266e-10 +1700 1.6252037027259294e-10 +1800 1.6566788375651406e-10 +1900 1.5323344558408191e-10 +2000 1.4401047170011939e-10 +2100 1.212825250135785e-10 +2200 1.2304483900713374e-10 +2300 1.2589379694554784e-10 +2400 1.2869605174722963e-10 +2500 1.3040046341793176e-10 +2600 1.3103764090659264e-10 +2700 1.3081041497769275e-10 +2800 1.012015957924326e-10 +2900 6.073299557657584e-11 +3000 5.648336848216792e-11 +3100 5.7709985851616314e-11 +3200 6.050919545922343e-11 +3300 6.36504231487108e-11 +3400 6.683343138071284e-11 +3500 6.942580949804323e-11 +3600 7.182106938726258e-11 +3700 7.6046280272992e-11 +3800 8.216136282416267e-11 +3900 8.940998938029766e-11 +4000 9.780836489304162e-11 +4100 1.0770898676266788e-10 +4200 1.1907453522129532e-10 +4300 1.3123577094033347e-10 +4400 1.43679722714824e-10 +4500 1.566199894577459e-10 +4600 1.7009308492159586e-10 +4700 1.8359572551573189e-10 +4800 1.9707733664831161e-10 +4900 2.1057042509066386e-10 +5000 2.2372854276830387e-10 +5100 2.3614710685974033e-10 +5200 2.475410682807095e-10 +5300 2.5778860232915414e-10 +5400 2.6686340045923035e-10 +5500 2.7473724504710184e-10 +5600 2.813065969883059e-10 +5700 2.8634283390600047e-10 +5800 2.89474072848681e-10 +5900 2.9021433988454443e-10 +6000 2.880852874733361e-10 +6100 2.8282326258814873e-10 +6200 2.74450733404262e-10 +6300 2.631361264698348e-10 +6400 2.4914450903915477e-10 +6500 2.32931852968441e-10 +6600 2.151054763693959e-10 +6700 1.962576921498402e-10 +6800 1.7688330464323465e-10 +6900 1.5740742475969613e-10 +7000 1.3822634473985118e-10 +7100 1.197282320008027e-10 +7200 1.0230146903941748e-10 +7300 8.63278604544093e-11 +7400 7.215434832434901e-11 +7500 6.004540164216144e-11 +7600 5.0132128860529355e-11 +7700 4.237955403496196e-11 +7800 3.658882715717464e-11 +7900 3.2436858120392024e-11 +8000 2.953886183950548e-11 +8100 2.751112568340469e-11 +8200 2.6012890795601455e-11 +8300 2.4759219987871305e-11 +8400 2.3512765007397322e-11 +8500 2.207280261092996e-11 +8600 2.028054200247841e-11 +8700 1.8046774938256726e-11 +8800 1.5388408546552442e-11 +8900 1.2446662283465413e-11 +9000 9.462236561053242e-12 +9100 6.707818298197342e-12 +9200 4.4050394069503076e-12 +9300 2.66536257066155e-12 +9400 1.4795188515502056e-12 +9500 7.508109573489653e-13 +9600 3.4736357568458737e-13 +9700 1.461917313569647e-13 +9800 5.5868361894366303e-14 +9900 1.935970498861533e-14 +10000 6.07595252916376e-15 +10100 1.7255963797025647e-15 +10200 4.432852019127042e-16 +10300 1.0307348990841779e-16 +10400 2.18131497166582e-17 +10500 4.318425398855969e-18 +10600 9.046740631850162e-19 +10700 2.8386954191429216e-19 +10800 1.6373418306861454e-19 +10900 4.69185403657041e-20 +11000 3.5724744937513105e-20 +11100 6.610797422691547e-22 +11200 4.305277152477696e-22 +11300 7.321893826745056e-25 +11400 3.1899841699314516e-25 +11500 8.494706281411055e-28 +11600 7.240419901481241e-29 +11700 2.5707293745456596e-30 +11800 1.2566957458750554e-31 +11900 5.4975869906767066e-33 +12000 2.1717755334422506e-34 +12100 7.722098087328418e-36 +12200 2.471337009242062e-37 +12300 7.117864066642182e-39 +12400 1.8450179388981175e-40 +12500 4.3036084859358695e-42 +12600 9.033282434738447e-44 +12700 1.7062462358090595e-45 +12800 2.899847120351541e-47 +12900 4.434674250399212e-49 +13000 6.101824653235464e-51 +13100 7.553888710868867e-53 +13200 8.407838639582765e-55 +13300 8.233513944741679e-57 +13400 0.0 +13500 0.0 +13600 0.0 +13700 0.0 +13800 0.0 +13900 0.0 +14000 0.0 +14100 0.0 +14200 0.0 +14300 0.0 +14400 0.0 +14500 0.0 +14600 0.0 +14700 0.0 +14800 0.0 +14900 0.0 +15000 0.0 +15100 0.0 +15200 0.0 +15300 0.0 +15400 0.0 +15500 0.0 +15600 0.0 +15700 0.0 +15800 0.0 +15900 0.0 +16000 0.0 +16100 0.0 +16200 0.0 +16300 0.0 +16400 0.0 +16500 0.0 +16600 0.0 +16700 0.0 +16800 0.0 +16900 0.0 +17000 0.0 diff --git a/tests/test_elements.py b/tests/test_elements.py index 7ee10942..680edae6 100644 --- a/tests/test_elements.py +++ b/tests/test_elements.py @@ -16,9 +16,9 @@ def test_isotopes(self): def test_abundance(self): hydrogen = elements.Element("H") - assert hydrogen.abundance("1") == 99.985 - assert hydrogen.abundance("2") == 0.015 + assert hydrogen.abundance("1") == pytest.approx(0.99985) + assert hydrogen.abundance("2") == pytest.approx(0.00015) carbon = elements.Element("C") - assert carbon.abundance("12") == 98.90 - assert carbon.abundance("13") == 1.10 + assert carbon.abundance("12") == pytest.approx(0.9890) + assert carbon.abundance("13") == pytest.approx(0.0110) diff --git a/tests/test_material.py b/tests/test_material.py new file mode 100644 index 00000000..7bdcc4c3 --- /dev/null +++ b/tests/test_material.py @@ -0,0 +1,79 @@ +import pytest + +from neucbot import elements, material + + +class TestComposition: + def test_from_file_isotopes_specified(self): + comp = material.Composition.from_file("./tests/test_material/WithIsotopes.dat") + + assert len(comp.materials) == 3 + assert comp.fractions.get("C") == pytest.approx(0.5) + assert comp.fractions.get("H") == pytest.approx(0.25) + assert comp.fractions.get("O") == pytest.approx(0.25) + + def test_from_file_no_isotopes_specified(self): + comp = material.Composition.from_file("./tests/test_material/NoIsotopes.dat") + + assert len(comp.materials) == 7 + assert comp.fractions.get("C") == pytest.approx(0.5) + assert comp.fractions.get("H") == pytest.approx(0.25) + assert comp.fractions.get("O") == pytest.approx(0.25) + + def test_normalize(self): + comp = material.Composition() + + comp.add(material.Isotope(elements.Element("C"), 12, 0.4)) + comp.add(material.Isotope(elements.Element("H"), 12, 0.2)) + comp.add(material.Isotope(elements.Element("O"), 12, 0.2)) + + comp.normalize() + + assert comp.fractions.get("C") == pytest.approx(0.5) + assert comp.fractions.get("H") == pytest.approx(0.25) + assert comp.fractions.get("O") == pytest.approx(0.25) + + def test_stopping_power_single_element_material(self): + comp = material.Composition.from_file("./tests/test_material/CarbonOnly.dat") + assert len(comp.materials) == 1 + assert comp.fractions == {"C": 1.0} + + # Since this material is 100% carbon 12, stopping power is just computed + # as the stopping power of carbon at this energy level + assert comp.stopping_power(11.0) == 486.0835 + + def test_stopping_power_multi_element_material(self): + comp = material.Composition.from_file("./tests/test_material/WithIsotopes.dat") + + assert len(comp.materials) == 3 + assert comp.fractions.get("C") == pytest.approx(0.5) + assert comp.fractions.get("H") == pytest.approx(0.25) + assert comp.fractions.get("O") == pytest.approx(0.25) + + # 50% stopping power of C = 0.5 * 486.0835 + # 25% stopping power of H = 0.25 * 1371.204 + # 25% stopping power of O = 0.25 * 462.8758 + assert comp.stopping_power(11.0) == 701.5617 + + +class TestStoppingPowerList: + def test_load_file(self): + stop_list = material.StoppingPowerList("C") + stop_list.load_file() + + assert len(stop_list.stopping_powers.items()) == 79 + assert stop_list.stopping_powers[0.01] == 511.72 + + def test_for_alpha(self): + stop_list = material.StoppingPowerList("C") + stop_list.load_file() + + # Energy less than lowest energy in the list + assert stop_list.for_alpha(0.001) == 511.72 + + # Energies in between two entries in the list + assert stop_list.for_alpha(0.525) == 1927.124 + assert stop_list.for_alpha(4.25) == 906.4551 + + # Energy higher than the highest in the list + assert stop_list.for_alpha(11.0) == 486.0835 diff --git a/tests/test_material/CarbonOnly.dat b/tests/test_material/CarbonOnly.dat new file mode 100644 index 00000000..491a6606 --- /dev/null +++ b/tests/test_material/CarbonOnly.dat @@ -0,0 +1 @@ +c 12 100.0 diff --git a/tests/test_material/NoIsotopes.dat b/tests/test_material/NoIsotopes.dat new file mode 100644 index 00000000..e336886d --- /dev/null +++ b/tests/test_material/NoIsotopes.dat @@ -0,0 +1,3 @@ +c 0 50.0 +o 0 25.0 +h 0 25.0 diff --git a/tests/test_material/WithIsotopes.dat b/tests/test_material/WithIsotopes.dat new file mode 100644 index 00000000..19ab6c06 --- /dev/null +++ b/tests/test_material/WithIsotopes.dat @@ -0,0 +1,3 @@ +c 12 50.0 +o 16 25.0 +h 1 25.0