In [2]:
#Imports
from math import *
import cmath
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import sys
import os
#Import custom modules
from physics import *
In [110]:
#Defined functions
B0 = 1710.0 #Magnetic field strength in Gauss
B0 = B0/10**4
yM = 0.5*.0254
def Radius(KE):
'''Radius of electron orbit in m given KE in keV'''
return me*c/(q*B0)*sqrt(1-(1/(KE*1000*q/(me*c**2)+1))**2)*(KE*1000*q/(me*c**2)+1)
def AngleIncidence(KE):
R = Radius(KE)
return asin((R-yM)/R)
In [116]:
#Create arrays over 2D parameter space (KE and AOI)
KEstepsize = 0.5
KE = np.arange(0.1,3.0+KEstepsize,KEstepsize)
ThetaStepsize = 30.0
ThetaStepsize *= 2*pi/360
Theta = np.arange(0,pi/2+ThetaStepsize,ThetaStepsize)
CosTheta = np.round(np.cos(Theta),6)
Nsegments = 100
SegmentPosition = np.linspace(0.0255,0.0365,Nsegments+2)
SegmentingSurfaces = ''
SurfaceSegmentNumber = range(100,100+Nsegments)
SurfaceSegmentNumberString = ''
line = ''
for x in range(Nsegments):
if len(line)<=45:
line += '-' + str(SurfaceSegmentNumber[x]) + ' '
SurfaceSegmentNumberString += '-' + str(SurfaceSegmentNumber[x]) + ' '
else:
line = ''
SurfaceSegmentNumberString += '\n\t-' + str(SurfaceSegmentNumber[x]) + ' '
SegmentingSurfaces += ' {Segment} pz {Position}\n'.format(Segment=100+x,Position=SegmentPosition[x+1])
for i in KE:
for j in CosTheta:
directory = os.curdir + '/MCNP_Decks/MCNP_{KE}MeV_{Theta}Degrees'\
.format(KE=str(round(i,3)),Theta=str(round(np.arccos(j)*360/(2*pi),2)))
file = open(directory, "w")
file.write('''c ************************
c Cell Cards
c ************************
c number material density then logical definition with surfaces
1 1 -2.70 1 -2 6 -7 8 -9 $Al layer
2 2 -1.44 2 -3 6 -7 8 -9 $Cellophane layer
3 3 -1.38 3 -4 6 -7 8 -9 $PET layer
4 4 -4.48 4 -5 6 -7 8 -9 $Phosphor layer
5 5 -0.000020615 -999 #1 #2 #3 #4
6 0 999
c ************************
c Surface Cards
c ************************
c number type position
1 pz 0 $Lanex surface
2 pz 0.0025 $End of Al layer
3 pz 0.008 $End of cellophane layer
4 pz 0.0255 $End of PET layer
5 pz 0.0365 $End of phosphor layer
6 px -0.1 $Lanex left edge
7 px 0.1 $Lanex right edge
8 py -1.5 $Lanex bottom edge
9 py 1.5 $Lanex top edge\n''')
file.write(SegmentingSurfaces)
file.write(''' 999 so 100
c ************************
c Problem Type
c ************************
mode p e
c **********************************
c Material Specification
c **********************************
m1 13000.02p -1 $aluminum -2.70
m2 1001.03e -0.062162 $cellophane -1.44
6012.03e -0.444462
8016.03e -0.493376
m3 1001.03e -0.143711 $polyethylene -.9
6000.03e -0.856289
m4 64000.03e 2 $Phosphor -4.48
8016.03e 2
16000.03e 1
m5 1001.03e 2 $water vapor -0.000020615
8016.03e 1
c
imp:p 1 4r 0 $ 1, 3
c
c *********************************
c Physics Simplification
c *********************************
phys:p 100 0 0 .9 0 $ pair production on change -1 to # for bias
phys:e 100 0 0 1 0 20 0 1 1 0 $ bremstrahlung biased
c *************************
c Tally
c *************************
F11:e 5 $particle counting
F16:e 4 $Energy deposition averaged over a cell\n''')
file.write('FS16 {Segments}T\n'.format(Segments=SurfaceSegmentNumberString))
file.write('''c **************************
c Energy Tally Card
c **************************
c e0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
c 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4
c 2.5 2.6 2.7 2.8 2.9 3.0
c ************************
c Source Definition
c ************************\n''')
file.write('sdef DIR={mu} VEC=0 0 1 POS=0 0 -0.1 CEL=5 ERG={KE} PAR=3\n'.format(KE=str(i),mu=str(j)))
file.write('''c sp1 -5 1
c
c Number of histories to transport
nps 20000\n''')
file.close()
In [139]:
#Create arrays over 1D parameter space (KE and AOI associated with particle tracking simulations)
KEstepsize = 0.01
KEsim = np.arange(1.00,3.00+KEstepsize,KEstepsize)
CosThetaSim = []
for x in KEsim:
CosThetaSim.append(np.cos(AngleIncidence(x*1000)))
Nsegments = 100
SegmentPosition = np.linspace(0.0255,0.0365,Nsegments+2)
SegmentingSurfaces = ''
SurfaceSegmentNumber = range(100,100+Nsegments)
SurfaceSegmentNumberString = ''
line = ''
for x in range(Nsegments):
if len(line)<=45:
line += '-' + str(SurfaceSegmentNumber[x]) + ' '
SurfaceSegmentNumberString += '-' + str(SurfaceSegmentNumber[x]) + ' '
else:
line = ''
SurfaceSegmentNumberString += '\n\t-' + str(SurfaceSegmentNumber[x]) + ' '
SegmentingSurfaces += ' {Segment} pz {Position}\n'.format(Segment=100+x,Position=SegmentPosition[x+1])
for i in range(len(KEsim)):
directory = os.curdir + '/MCNP_Decks/MCNP_{KE}MeV_{Theta}Degrees'\
.format(KE=str(round(KEsim[i],3)),Theta=str(round(np.arccos(CosThetaSim[i])*360/(2*pi),1)))
file = open(directory, "w")
file.write('''c ************************
c Cell Cards
c ************************
c number material density then logical definition with surfaces
1 1 -2.70 1 -2 6 -7 8 -9 $Al layer
2 2 -1.44 2 -3 6 -7 8 -9 $Cellophane layer
3 3 -1.38 3 -4 6 -7 8 -9 $PET layer
4 4 -4.48 4 -5 6 -7 8 -9 $Phosphor layer
5 5 -0.000020615 -999 #1 #2 #3 #4
6 0 999
c ************************
c Surface Cards
c ************************
c number type position
1 pz 0 $Lanex surface
2 pz 0.0025 $End of Al layer
3 pz 0.008 $End of cellophane layer
4 pz 0.0255 $End of PET layer
5 pz 0.0365 $End of phosphor layer
6 px -0.1 $Lanex left edge
7 px 0.1 $Lanex right edge
8 py -1.5 $Lanex bottom edge
9 py 1.5 $Lanex top edge\n''')
file.write(SegmentingSurfaces)
file.write(''' 999 so 100
c ************************
c Problem Type
c ************************
mode p e
c **********************************
c Material Specification
c **********************************
m1 13000.02p -1 $aluminum -2.70
m2 1001.03e -0.062162 $cellophane -1.44
6012.03e -0.444462
8016.03e -0.493376
m3 1001.03e -0.143711 $polyethylene -.9
6000.03e -0.856289
m4 64000.03e 2 $Phosphor -4.48
8016.03e 2
16000.03e 1
m5 1001.03e 2 $water vapor -0.000020615
8016.03e 1
c
imp:p 1 4r 0 $ 1, 3
c
c *********************************
c Physics Simplification
c *********************************
phys:p 100 0 0 .9 0 $ pair production on change -1 to # for bias
phys:e 100 0 0 1 0 20 0 1 1 0 $ bremstrahlung biased
c *************************
c Tally
c *************************
F11:e 5 $particle counting
F16:e 4 $Energy deposition averaged over a cell\n''')
file.write('FS16 {Segments}T\n'.format(Segments=SurfaceSegmentNumberString))
file.write('''c **************************
c Energy Tally Card
c **************************
c e0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
c 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4
c 2.5 2.6 2.7 2.8 2.9 3.0
c ************************
c Source Definition
c ************************\n''')
file.write('sdef DIR={mu} VEC=0 0 1 POS=0 0 -0.1 CEL=5 ERG={KE} PAR=3\n'\
.format(KE=str(KEsim[i]),mu=str(CosThetaSim[i])))
file.write('''c sp1 -5 1
c
c Number of histories to transport
nps 20000\n''')
file.close()