Python 3 Tkinter Radiation & Thermal Electricity Shielding Management System Using Matplotlib & Numpy Library GUI Desktop App
"""
A simple example of an animated plot... In 3D!
see sources in Monte function
TODO
-why is there a tendency to drift right?
-is there a way to make a sort of heatmap of flux on the surface plot?
-Make parameters enterable - Tkinter?
-writeup - include flux vs. distance, flux vs. number of particles, flux vs. isotropy factor
-write about the assumptions - uses standard monte carlo but also two-term Henyey-Greenstein phase functions
-how it differs from the 741
-real attenuation absorption scattering coeff.
"""
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation
from matplotlib import cm
from random import randint,random
from tkinter import *
fields = ('mu_s (cm-1)', 'mu_a (cm-1)', 'Isotropy', 'Number Particles', 'World Size (cm)','Flux Width x (cm)','Flux Width y (cm)','Flux Distance z (cm)',
'Russian Roulette Factor','Initial E-Gamma (MeV)','Cutoff E-Gamma (MeV)','Flux v Distance Range (cm)','ALPHA')
def get_mu_a(e_gamma_0,mu_s):
MATCHER=[
( 1.00E-003 , 3.71E+003 ),
( 1.04E-003 , 3.39E+003 ),
( 1.07E-003 , 3.09E+003 ),
( 1.50E-003 , 1.25E+003 ),
( 2.00E-003 , 5.60E+002 ),
( 2.15E-003 , 4.58E+002 ),
( 2.30E-003 , 3.80E+002 ),
( 2.47E-003 , 3.10E+002 ),
( 2.64E-003 , 2.61E+002 ),
( 2.82E-003 , 2.16E+002 ),
( 3.00E-003 , 1.84E+002 ),
( 3.61E-003 , 1.07E+002 ),
( 4.00E-003 , 8.16E+001 ),
( 5.00E-003 , 4.22E+001 ),
( 6.00E-003 , 2.46E+001 ),
( 8.00E-003 , 1.04E+001 ),
( 1.00E-002 , 5.38E+000 ),
( 1.50E-002 , 1.70E+000 ),
( 2.00E-002 , 8.23E-001 ),
( 3.00E-002 , 3.79E-001 ),
( 4.00E-002 , 2.69E-001 ),
( 5.00E-002 , 2.26E-001 ),
( 6.00E-002 , 2.05E-001 ),
( 8.00E-002 , 1.82E-001 ),
( 1.00E-001 , 1.69E-001 ),
( 1.50E-001 , 1.49E-001 ),
( 2.00E-001 , 1.36E-001 ),
( 3.00E-001 , 1.18E-001 ),
( 4.00E-001 , 1.05E-001 ),
( 5.00E-001 , 9.60E-002 ),
( 6.00E-001 , 8.87E-002 ),
( 8.00E-001 , 7.79E-002 ),
( 1.00E+000 , 7.01E-002 ),
( 1.25E+000 , 6.27E-002 ),
( 1.50E+000 , 5.70E-002 ),
( 2.00E+000 , 4.90E-002 ),
( 3.00E+000 , 3.93E-002 ),
( 4.00E+000 , 3.37E-002 ),
( 5.00E+000 , 3.00E-002 ),
( 6.00E+000 , 2.74E-002 ),
( 8.00E+000 , 2.40E-002 ),
( 1.00E+001 , 2.19E-002 ),
( 1.50E+001 , 1.92E-002 ),
( 2.00E+001 , 1.79E-002 ),
]
for e,mu_s in MATCHER:
if e==e_gamma_0:
mu=mu_s
return mu
def runit(entries):
mu_s = float(entries['mu_s (cm-1)'].get())
mu_a = float(entries['mu_a (cm-1)'].get())
e_gamma_0 = float(entries['Initial E-Gamma (MeV)'].get())
mu_s=get_mu_a(e_gamma_0,mu_s)
num_particles = int(entries['Number Particles'].get())
fluxsurface_len_x = int(entries['Flux Width x (cm)'].get())
fluxsurface_len_y = int(entries['Flux Width y (cm)'].get())
flux_dist_from_src_z = float(entries['Flux Distance z (cm)'].get())
g = float(entries['Isotropy'].get())
lim = float(entries['World Size (cm)'].get())
m = float(entries['Russian Roulette Factor'].get())
lines_or_points=int(v0.get())
rand_init_dir=int(v1.get())
plot_flux_surface=int(v2.get())
sampling=v5.get()
cutoff_e = float(entries['Cutoff E-Gamma (MeV)'].get())
fig = plt.figure()
ax = p3.Axes3D(fig)
ax.set_xlim3d([-lim, lim])
ax.set_xlabel('X')
ax.set_ylim3d([-lim, lim])
ax.set_ylabel('Y')
ax.set_zlim3d([-lim, lim])
ax.set_zlabel('Z')
ax.set_title('3D Monte Carlo Walking')
flux_surface_area = float(fluxsurface_len_x)*float(fluxsurface_len_y)
def Monte(length, dims=6,g=0.0,mu_a=0.03,mu_s=0.3,e_gamma_0=5.0,rand_init_dir=False, lines_or_points=True,sampling='H'):
'''
from http://en.wikipedia.org/wiki/Monte_Carlo_method_for_photon_transport
AND
'Successive order, multiple scattering of two-term Henyey-Greenstein phase functions'
AND
Wang, Jaques MCML Monte Carlo method for light transport 1995
'''
e_gamma=e_gamma_0
mu_t = mu_a+mu_s
lineData = np.empty((dims, length))
lineData[:, 0] = 0.0,0.0,0.0, 0,0,0
if rand_init_dir:
x_dir,y_dir,z_dir = randint(0,10),randint(0,10),randint(0,10)
mag = np.sqrt(x_dir**2+y_dir**2+z_dir**2)
mu_x,mu_y,mu_z = x_dir/mag,y_dir/mag,z_dir/mag
else:
mu_x,mu_y,mu_z = 0.0,0.0,1.0
single_p_crossing_factor=0
x_cross,y_cross=0,0
for index in range(1, length):
alpha = e_gamma/0.511
zeta1,zeta2,zeta3,zeta4=random(),random(),random(),random()
step_scalar = -np.log(zeta1)/mu_t
if (sampling == 'H') or (sampling == 'B'):
ct = (1/(2*g))*(1+g**2-((1-g**2)/(1-g+2*g*zeta2))**2) if g != 0 else 2*zeta2-1
elif sampling =='K':
while 1:
ct = -1 + 2*random()
trm1 = (1 + ct**2)/2.0
trm2 = 1.0/(1 + alpha*(1 - ct))
trm3 = 2*alpha**2*(1 - ct)**2*trm2/trm1
ang_diff_scat = trm1*trm2**2*(1 + trm3)
if ang_diff_scat<random():
break
phi = 2*np.pi*zeta3
theta = np.arccos(ct)
c = np.sqrt(1-mu_z**2)
st = np.sin(theta)
sp = np.sin(phi)
cp = np.cos(phi)
if mu_z >= 0.999:
mu_x,mu_y,mu_z = st*cp,st*sp,ct
elif mu_z <= -0.999:
mu_x,mu_y,mu_z = st*cp,-st*sp,-ct
else:
new_mu_x = (st/c)*(mu_x*mu_z*cp-mu_y*sp)+mu_x*ct
new_mu_y = (st/c)*(mu_y*mu_z*cp-mu_x*sp)+mu_y*ct
new_mu_z = -c*st*cp+mu_z*ct
mu_x,mu_y,mu_z = new_mu_x,new_mu_y,new_mu_z
x,y,z = mu_x*step_scalar,mu_y*step_scalar,mu_z*step_scalar
old_z = lineData[2,index-1]
avg_x,avg_y = lineData[0,index-1]+x/2.0,lineData[1,index-1]+y/2.0
new_z = old_z+z
if ((-fluxsurface_len_x/2<avg_x<fluxsurface_len_x/2) and (-fluxsurface_len_y/2<avg_y<fluxsurface_len_y/2) and ((old_z < flux_dist_from_src_z and new_z > flux_dist_from_src_z) or (old_z > flux_dist_from_src_z and new_z < flux_dist_from_src_z))):
single_p_crossing_factor+=1*abs(1/np.dot([mu_x,mu_y,mu_z],[0,0,1]))
x_cross,y_cross=lineData[0, index-1] + x,lineData[1, index-1] + y
step = x,y,z
de = (mu_a/mu_t)*e_gamma
e_gamma=e_gamma-de
if e_gamma<cutoff_e:
if zeta4 < (1.0/m):
e_gamma=m*e_gamma
lineData[:3, index] = lineData[:3, index-1] + step
lineData[3, index]=single_p_crossing_factor
lineData[4, index]=x_cross
lineData[5, index]=y_cross
else:
e_gamma=0
lineData[:3, index] = lineData[:3, index-1] + step
lineData[3, index]=single_p_crossing_factor
lineData[4, index]=x_cross
lineData[5, index]=y_cross
final = lineData[:, :index]
break
else:
lineData[:3, index] = lineData[:3, index-1] + step
lineData[3, index]=single_p_crossing_factor
lineData[4, index]=x_cross
lineData[5, index]=y_cross
final = lineData[:, :index]
continue
return final if lines_or_points else final[:, -1]
if not lines_or_points:
total_crossing_factor=0
gx=[]
gy=[]
for j in range(num_particles):
ptdata = [Monte(25, 6,g,mu_a,mu_s,e_gamma_0,rand_init_dir,lines_or_points,sampling) for index in range(1)]
for i,d in enumerate(ptdata):
total_crossing_factor = total_crossing_factor+d[3]
ax.scatter(d[0],d[1],d[2])
flux = total_crossing_factor/(num_particles*flux_surface_area)
if plot_flux_surface:
x = np.arange(-fluxsurface_len_x/2, fluxsurface_len_x/2, 0.1)
y = np.arange(-fluxsurface_len_y/2, fluxsurface_len_y/2, 0.1)
X, Y = np.meshgrid(x, y)
zs = np.array([flux_dist_from_src_z for x,y in zip(np.ravel(X), np.ravel(Y))])
Z = zs.reshape(X.shape)
ax.plot_surface(X, Y, Z)
def update_lines(num, dataLines, lines) :
for line, data in zip(lines, dataLines):
line.set_data(data[0:2, :num])
line.set_3d_properties(data[2,:num])
return 1
if lines_or_points:
data = [Monte(25, 6,g,mu_a,mu_s,e_gamma_0,rand_init_dir,True,sampling) for index in range(num_particles)]
lines = [ax.plot(dat[0, 0:1], dat[1, 0:1], dat[2, 0:1])[0] for i,dat in enumerate(data)]
line_ani = animation.FuncAnimation(fig, update_lines, 25, fargs=(data, lines),interval=50, blit=False)
plt.legend()
plt.show()
def fluxit(entries):
mu_s = float(entries['mu_s (cm-1)'].get())
mu_a = float(entries['mu_a (cm-1)'].get())
e_gamma_0 = float(entries['Initial E-Gamma (MeV)'].get())
mu_s=get_mu_a(e_gamma_0,mu_s)
num_particles = int(entries['Number Particles'].get())
fluxsurface_len_x = int(entries['Flux Width x (cm)'].get())
fluxsurface_len_y = int(entries['Flux Width y (cm)'].get())
flux_dist_from_src_z = float(entries['Flux Distance z (cm)'].get())
g = float(entries['Isotropy'].get())
lim = float(entries['World Size (cm)'].get())
m = float(entries['Russian Roulette Factor'].get())
lines_or_points=int(v0.get())
rand_init_dir=int(v1.get())
plot_flux_surface=int(v2.get())
sampling=v5.get()
cutoff_e = float(entries['Cutoff E-Gamma (MeV)'].get())
fluxrange = int(entries['Flux v Distance Range (cm)'].get())
fig, ax = plt.subplots()
'''
ax = p3.Axes3D(fig)
ax.set_xlim3d([-lim, lim])
ax.set_xlabel('X')
ax.set_ylim3d([-lim, lim])
ax.set_ylabel('Y')
ax.set_zlim3d([-lim, lim])
ax.set_zlabel('Z')
ax.set_title('Flux v D')
'''
flux_surface_area = float(fluxsurface_len_x)*float(fluxsurface_len_y)
def FluxMonte(length, dims=6,g=0.0,mu_a=0.03,mu_s=0.3,e_gamma_0=5.0,rand_init_dir=False, dis=1.0, sampling='H'):
'''
from http://en.wikipedia.org/wiki/Monte_Carlo_method_for_photon_transport
AND
'Successive order, multiple scattering of two-term Henyey-Greenstein phase functions'
AND
Wang, Jaques MCML Monte Carlo method for light transport 1995
'''
e_gamma=e_gamma_0
mu_t = mu_a+mu_s
lineData = np.empty((dims, length))
lineData[:, 0] = 0.0,0.0,0.0, 0,0,0
if rand_init_dir:
x_dir,y_dir,z_dir = randint(0,10),randint(0,10),randint(0,10)
mag = np.sqrt(x_dir**2+y_dir**2+z_dir**2)
mu_x,mu_y,mu_z = x_dir/mag,y_dir/mag,z_dir/mag
else:
mu_x,mu_y,mu_z = 0.0,0.0,1.0
single_p_crossing_factor=0
x_cross,y_cross=0,0
for index in range(1, length):
alpha = e_gamma/0.511
zeta1,zeta2,zeta3,zeta4=random(),random(),random(),random()
step_scalar = -np.log(zeta1)/mu_t
if (sampling == 'H') or (sampling == 'B'):
ct = (1/(2*g))*(1+g**2-((1-g**2)/(1-g+2*g*zeta2))**2) if g != 0 else 2*zeta2-1
elif sampling =='K':
while 1:
ct = -1 + 2*random()
trm1 = (1 + ct**2)/2.0
trm2 = 1.0/(1 + alpha*(1 - ct))
trm3 = 2*alpha**2*(1 - ct)**2*trm2/trm1
ang_diff_scat = trm1*trm2**2*(1 + trm3)
if ang_diff_scat<random():
break
phi = 2*np.pi*zeta3
theta = np.arccos(ct)
c = np.sqrt(1-mu_z**2)
st = np.sin(theta)
sp = np.sin(phi)
cp = np.cos(phi)
if mu_z == 1:
mu_x,mu_y,mu_z = st*cp,st*sp,ct
elif mu_z == -1:
mu_x,mu_y,mu_z = st*cp,-st*sp,-ct
else:
new_mu_x = (st/c)*(mu_x*mu_z*cp-mu_y*sp)+mu_x*ct
new_mu_y = (st/c)*(mu_y*mu_z*cp-mu_x*sp)+mu_y*ct
new_mu_z = -c*st*cp+mu_z*ct
mu_x,mu_y,mu_z = new_mu_x,new_mu_y,new_mu_z
x,y,z = mu_x*step_scalar,mu_y*step_scalar,mu_z*step_scalar
old_z = lineData[2,index-1]
avg_x,avg_y = lineData[0,index-1]+x/2.0,lineData[1,index-1]+y/2.0
new_z = old_z+z
if ((-fluxsurface_len_x/2<avg_x<fluxsurface_len_x/2) and (-fluxsurface_len_y/2<avg_y<fluxsurface_len_y/2) and ((old_z < dis and new_z > dis) or (old_z > dis and new_z < dis))):
single_p_crossing_factor+=1*abs(1/np.dot([mu_x,mu_y,mu_z],[0,0,1]))
x_cross,y_cross=lineData[0, index-1] + x,lineData[1, index-1] + y
step = x,y,z
de = (mu_a/mu_t)*e_gamma
e_gamma=e_gamma-de
if e_gamma<cutoff_e:
if zeta4 < (1.0/m):
e_gamma=m*e_gamma
lineData[:3, index] = lineData[:3, index-1] + step
lineData[3, index]=single_p_crossing_factor
lineData[4, index]=x_cross
lineData[5, index]=y_cross
else:
e_gamma=0
lineData[:3, index] = lineData[:3, index-1] + step
lineData[3, index]=single_p_crossing_factor
lineData[4, index]=x_cross
lineData[5, index]=y_cross
final = lineData[:, :index]
break
else:
lineData[:3, index] = lineData[:3, index-1] + step
lineData[3, index]=single_p_crossing_factor
lineData[4, index]=x_cross
lineData[5, index]=y_cross
final = lineData[:, :index]
continue
return final[:, -1]
to_plot=[]
if sampling == 'B':
to_plot = ['H','K']
if sampling == 'H':
to_plot = ['H']
if sampling == 'K':
to_plot = ['K']
for p in to_plot:
dis_a=[]
flux_a=[]
for dis in range(fluxrange*4):
dis = dis/4.0
total_crossing_factor=0
gx=[]
gy=[]
for j in range(num_particles):
ptdata = [FluxMonte(25, 6,g,mu_a,mu_s,e_gamma_0,rand_init_dir,dis,p) for index in range(1)]
for i,d in enumerate(ptdata):
total_crossing_factor = total_crossing_factor+d[3]
flux = total_crossing_factor/(num_particles*flux_surface_area)
dis_a.append(dis)
flux_a.append(flux)
plt.plot(dis_a,flux_a,label=p)
ax.set_xlabel('Distance (cm)')
ax.set_ylabel('Scalar Flux')
ax.set_title('Scalar Flux v Distance')
plt.legend()
plt.show()
def makeform(root, fields):
entries = {}
for field in fields:
row = Frame(root)
lab = Label(row, width=22, text=field+": ", anchor='w')
ent = Entry(row)
if field == 'mu_s (cm-1)':
ent.insert(0,"0.16")
elif field =='mu_a (cm-1)':
ent.insert(0,"0.38")
elif field == 'Number Particles':
ent.insert(0,"500")
elif field == 'Isotropy':
ent.insert(0,"0.8")
elif field == 'World Size (cm)':
ent.insert(0,"40")
elif field == 'Flux Width x (cm)':
ent.insert(0,"1")
elif field == 'Flux Width y (cm)':
ent.insert(0,"1")
elif field == 'Flux Distance z (cm)':
ent.insert(0,"4")
elif field == 'Russian Roulette Factor':
ent.insert(0,"1.4")
elif field =='Cutoff E-Gamma (MeV)':
ent.insert(0,"0.03")
elif field == 'Initial E-Gamma (MeV)':
ent.insert(0,"1.25")
elif field == 'Flux v Distance Range (cm)':
ent.insert(0,"5")
elif field == 'ALPHA':
ent.insert(0,"2.0")
row.pack(side=TOP, fill=X, padx=5, pady=5)
lab.pack(side=LEFT)
ent.pack(side=RIGHT, expand=YES, fill=X)
entries[field] = ent
return entries
if __name__ == '__main__':
root = Tk()
ents = makeform(root, fields)
MODES0 = [
("Walking Animation", "1"),
("Flux Scatter Plot", "0"),
]
v0 = StringVar()
v0.set("0")
for text, mode in MODES0:
b0 = Radiobutton(root, text=text,
variable=v0, value=mode)
b0.pack(anchor=W)
MODES1 = [
("Random Init Dir", "1"),
("Z-Directed", "0"),
]
v1 = StringVar()
v1.set("0")
for text, mode in MODES1:
b1 = Radiobutton(root, text=text,
variable=v1, value=mode)
b1.pack(anchor=E)
MODES2 = [
("Plot Flux Surface", "1"),
("Don't Plot Flux Surface", "0"),
]
v2 = StringVar()
v2.set("0")
for text, mode in MODES2:
b2 = Radiobutton(root, text=text,
variable=v2, value=mode)
b2.pack(anchor=W)
MODES5 = [
("Klein Nishina", "K"),
("Henyey Greenstein", "H"),
("Both (Compare - only for flux v distance)", "B"),
]
v5 = StringVar()
v5.set("H")
for text, mode in MODES5:
b5 = Radiobutton(root, text=text,
variable=v5, value=mode)
b5.pack(anchor=E)
root.bind('<Return>', (lambda event, e=ents: fetch(e)))
b2 = Button(root, text='3d Visual',
command=(lambda e=ents: runit(e)))
b2.pack(side=LEFT, padx=5, pady=5)
b4 = Button(root, text='Flux v Distance',
command=(lambda e=ents: fluxit(e)))
b4.pack(side=LEFT, padx=5, pady=5)
b3 = Button(root, text='Quit', command=root.quit)
b3.pack(side=LEFT, padx=5, pady=5)
root.mainloop()