Upload New File

import os

import numpy as np

import matplotlib.pyplot as plt

import pandas as pd

from tabulate import tabulate

def interpolate(energy, dataframe):

    length = len(dataframe)

    i = 1

    while i < length: 

        if energy == dataframe.loc[i,'Energy Level (MeV)']:

           x = dataframe.loc[i, 'Energy Level (MeV)']

           return x

        elif (energy > dataframe.loc[i,'Energy Level (MeV)']) and (energy < dataframe.loc[i+1,'Energy Level (MeV)']):

         top_bound = dataframe.loc[i+1,'Energy Level (MeV)']

         bottom_bound = dataframe.loc[i,'Energy Level (MeV)']

         decimal = (top_bound-energy)/(top_bound-bottom_bound)

         top_attenuation = dataframe.loc[i+1, 'Mass Attenuation']

         bottom_attenuation = dataframe.loc[i, 'Mass Attenuation']

         x = (-1*decimal*(top_attenuation-bottom_attenuation))+top_attenuation

         return x

        else:

           i += 1 

def element_attenuation(energies, element, empty, index):

    attenuation_list = []

    end = len(energies)

    for i in energies:

        dataframe = pd.read_excel('tabulated_attenuation.xlsx', element)

        mass_attenuation = interpolate(i,dataframe)

        attenuation_list.append(mass_attenuation)

    empty.insert(index, element, attenuation_list, True)

def element_iterate(element_list, energies, contribution):

   list = []

   data = []

   empty = pd.DataFrame(data=data)

   for i in element_list:

        j = 0

        name = str(element_list[j])

        list.append(element_attenuation(energies, i, empty,j ))

        j+=1

   k = 0

   while k < len(element_list):

    name = str(element_list[k])

    empty.loc[:, name] *= contribution[k]

    k+=1 

   empty = empty.cumsum(axis='columns')

   mixture_list = empty[str(element_list[0])].tolist()

   return mixture_list

def percent(element_list, energies, contribution, distance, flux_weight, density):

   percent_attenuation = []

   attenuations = element_iterate(element_list,energies, contribution)

   length = len(attenuations)

   i = 0

   while i < len(attenuations):

      mu = attenuations[i] * density

      percent = energy[i] * np.exp(distance * mu * -1)

      weighted_percent = percent * flux_weight[i]

      raw_count = percent*fluence[i]*54000

      percent_attenuation.append(raw_count)

      i+=1

   return percent_attenuation 

#NaCl-UCl3 (Mol%: 63, 37)

salt_list = ['Na','Cl','U']

salt_contribution = [0.08819, 0.37557, 0.53632]

salt_density = 3.375837

salt_density2 = 3.0712227

#SS 316L

stainless_list = ['C','Mn','Si','Cr','Ni','Mo','S','P','Fe']

stainless_cont = [0.00030,0.02000,0.01000,0.17000,0.12000,0.02500,0.00045,0.00030,0.65395]

stainless_density = 8.027

#Tungsten carbide

tungsten_list = ['W', 'C']

tungsten_cont = [0.938674, 0.061326]

#Energy Spectrum of XRS4 (Use this list)

energy = [0.0101,0.0123,0.0154,0.0185,0.022,0.026,0.031,0.0366,0.0438,0.0521,0.0615,0.0733,0.087,0.1034,0.123,0.14607,0.1736,0.2061,0.2451,0.2913,0.3459]

#These are inputs to the percent() function, fluence as an argument will give you the raw count of photons

fluence_percent   = [0.00121317452223333,0.00729311050848389,0.013374026839684,0.0239236822325581,0.0376913199223123,0.0564644355336162,0.077919938317699,0.0966948512280771,0.106891582439115,0.100644661029792,0.213267342396772,0.0727832575637353,0.0584984879130963,0.0452893200629327,0.0331560807948942,0.0231717576562574,0.014622496132929,0.00911736906876152,0.00487067813815931,0.00206215324853904,0.00105027386978412]

fluence = [1050.40277526783,6314.59314224956,11579.6323212009,20713.8393950476,32634.2717592037,48888.5965716274,67465.4124016195,83721.2933693209,92549.9281341009,87141.1567949589,184653.241740663,63017.9206180832,50649.7399399474,39212.847461157,28707.5261190352,20062.8006143416,12660.5943645243,7894.09074906888,4217.17876529466,1785.47394096723,909.358519647794]

kev_energy = [10,12,15,19,22,26,31,37,44,52,62,73,87,103,123,146,174,206,245,291,346]

salt = percent(salt_list,energy,salt_contribution, 0.652, fluence_percent, salt_density)

steel = percent(stainless_list,energy,stainless_cont, 0.652, fluence_percent, stainless_density)

#energy_bin is for graphing purposes, it indicates the unattenuated energy bins, x is the ticks for the x axis of a bar graph

energy_bin = []

for i in range(0, len(energy)):

   photons = energy[i] * fluence[i] * 54000

   energy_bin.append(photons)

x = np.arange(len(energy))

#Graph stuff below here

myList = list(np.around(np.array(energy),3))

#MAKE FONT BIGGER

plt.rcParams.update({'font.size': 16})

plt.xlabel('Average Bin Energy (MeV)')

plt.ylabel('Photon Fluence (cm$^-$$^2$)')

plt.title('Photons Produced from Scintillator Reaction of Detector Plate')

plt.bar(x - 0.35, salt, width=0.25, align='edge', label= 'Through NaCl-UCl3', color='tab:orange')

plt.bar(x, steel, width=0.25, align='center', label= 'Through SS 316L', color='tab:blue')

plt.bar(x + 0.2, energy_bin, width=0.25, align='center', label='Unattenuated', color='tab:green')

plt.xticks(x, myList)

plt.legend()

plt.show()

total_salt = sum(salt)

total_steel = sum(steel)

total = sum(energy_bin)

position = [1,2,3]

totals = [total_salt, total_steel, total]

labels = ['NaCl-UCl$_3$','SS 316L', 'Unattenuated']

bars = plt.bar(position, totals, align='center')

plt.xticks(position, labels)

bars[0].set_color('tab:orange')

bars[1].set_color('tab:blue')

bars[2].set_color('tab:green')

plt.title('Photons per cm$^2$ for All Energy Bins')

plt.yscale('log')

plt.show()

#Print statements for ease of use

print('The total amount of photons per cm$^2$ produced after passing through the salt is ' +str(round(total_salt,0)))

print('The total amount of photons per cm$^2$ produced after passing through the steel is ' +str(round(total_steel,0)))

print('The total amount of photons per cm$^2$ that would be produced after passing through unattenuated is ' +str(round(total,0)))

 No newline at end of file