Add new python script and update plot_tally.py
This commit is contained in:
58
analyze_materials.py
Normal file
58
analyze_materials.py
Normal file
@@ -0,0 +1,58 @@
|
||||
from dataclasses import dataclass
|
||||
from typing import List, Union, Optional
|
||||
import periodictable
|
||||
|
||||
@dataclass
|
||||
class Composition:
|
||||
element: int
|
||||
isotope: int
|
||||
abundance: float
|
||||
library: Optional[str] = None
|
||||
|
||||
def __repr__(self):
|
||||
ele = periodictable.elements[self.element]
|
||||
return f'{self.isotope}{ele} ({self.abundance})'
|
||||
|
||||
@dataclass
|
||||
class Material:
|
||||
ID: int
|
||||
elements: List[Composition]
|
||||
|
||||
def get_elements(istring):
|
||||
ilines = istring.splitlines()
|
||||
alines = []
|
||||
for li in ilines:
|
||||
alines.append(li.split('$', 1)[0])
|
||||
tstring = ' '.join(alines)
|
||||
items = tstring.split()
|
||||
elements = []
|
||||
for i in range(len(items)//2):
|
||||
eil = items[i*2]
|
||||
if '.' in eil:
|
||||
ei, l = eil.split('.', 1)
|
||||
else:
|
||||
ei = eil
|
||||
l = None
|
||||
ei = int(ei)
|
||||
a = float(items[i*2+1])
|
||||
comp = Composition(ei//1000, ei%1000, a, l)
|
||||
elements.append(comp)
|
||||
return elements
|
||||
|
||||
|
||||
txt=open('geometry/materials.i', 'r').read()
|
||||
lns=[]
|
||||
|
||||
while txt.find('\n')>0:
|
||||
nn = txt.find('\n')
|
||||
while len(txt)>(nn+1) and txt[nn+1] in [' ', '\t'] and txt[nn+1:].find('\n')>0:
|
||||
nn = nn+1+txt[nn+1:].find('\n')
|
||||
lns.append(txt[:nn+1])
|
||||
txt = txt[nn+1:]
|
||||
|
||||
flns = [li for li in lns if not li.startswith('c')]
|
||||
mat={}
|
||||
for fi in flns:
|
||||
if fi[0]=='m' and fi[1].isdigit():
|
||||
mat[int(fi.split()[0][1:])] = fi.split(None, 1)[1]
|
||||
|
||||
@@ -61,11 +61,12 @@ if BUNKER:
|
||||
("px_11.15_6x6.png", 6.0)]
|
||||
else:
|
||||
GEO_PATH="plot_geometry"
|
||||
XY_PLANES=[0., -30., 161., 240.]
|
||||
XY_PLANES=[0., -30., 161., 240., -170.]
|
||||
XY_GEO=[("pz_0.1_16x16.png", 32.0), ("pz_-0.3_16x16.png", 32.0),
|
||||
("pz_1.6_16x16.png", 32.0), ("pz_2.4_16x16.png", 32.0)]
|
||||
XZ_PLANES=[0.]
|
||||
XZ_GEO=[("py_0.0_16x16.png", 32.0)]
|
||||
("pz_1.6_16x16.png", 32.0), ("pz_2.4_16x16.png", 32.0),
|
||||
("pz_-17.0_16x16.png", 32.0)]
|
||||
XZ_PLANES=[0., 330.]
|
||||
XZ_GEO=[("py_0.0_16x16.png", 32.0), ("py_3.3_16x16.png", 32.0)]
|
||||
YZ_PLANES=[1500., 1700., 2000., 2180., 2300., 2900., 3500.]
|
||||
YZ_GEO=[("px_15.0_6x6.png", 12.0), ("px_17.0_6x6.png", 12.0),
|
||||
("px_20.0_6x6.png", 12.0), ("px_21.8_6x6.png", 12.0),
|
||||
@@ -176,6 +177,10 @@ def plot_xy(outfile, res, zidx=5, geo=None):
|
||||
ax.pcolormesh(imgx, imgy,
|
||||
where(img[:,:,1], log10(2e-3), nan), alpha=0.2,
|
||||
norm=Rnorm, cmap=Rcmap)
|
||||
# red lines for instrument boundary
|
||||
ax.pcolormesh(imgx, imgy,
|
||||
where(img[:,:,2]-img[:,:,0], log10(1.0), nan),
|
||||
norm=Rnorm, cmap='autumn')
|
||||
elif geometry is not None:
|
||||
ax.pcolormesh(geometry['x'], geometry['y'],
|
||||
log10(geometry['lines']*0+2e-3),
|
||||
@@ -226,6 +231,10 @@ def plot_xz(outfile, res, yidx=5, geo=None):
|
||||
ax.pcolormesh(imgx, imgy,
|
||||
where(img[:,:,1], log10(2e-3), nan), alpha=0.2,
|
||||
norm=Rnorm, cmap=Rcmap)
|
||||
# red lines for instrument boundary
|
||||
ax.pcolormesh(imgx, imgy,
|
||||
where(img[:,:,2]-img[:,:,0], log10(1.0), nan),
|
||||
norm=Rnorm, cmap='autumn')
|
||||
|
||||
xlim(x[0]/100, x[-1]/100)
|
||||
ylim(z[0]/100, z[-1]/100)
|
||||
@@ -271,6 +280,10 @@ def plot_yz(outfile, res, xidx=5, geo=None):
|
||||
ax.pcolormesh(imgx, imgy,
|
||||
where(img[:,:,1], log10(2e-3), nan), alpha=0.2,
|
||||
norm=Rnorm, cmap=Rcmap)
|
||||
# red lines for instrument boundary
|
||||
ax.pcolormesh(imgx, imgy,
|
||||
where(img[:,:,2]-img[:,:,0], log10(1.0), nan),
|
||||
norm=Rnorm, cmap='autumn')
|
||||
|
||||
xlim(y[0]/100, y[-1]/100)
|
||||
ylim(z[0]/100, z[-1]/100)
|
||||
@@ -341,14 +354,22 @@ if __name__=='__main__':
|
||||
output_path=os.path.join('.', 'images', sys.argv.pop(idx))
|
||||
else:
|
||||
output_path=None
|
||||
if '-s' in sys.argv:
|
||||
idx=sys.argv.index('-s')
|
||||
sys.argv.pop(idx)
|
||||
SCALING=float(sys.argv.pop(idx))
|
||||
maxadd=[]
|
||||
scale=SCALING
|
||||
while '-m' in sys.argv:
|
||||
idx=sys.argv.index('-m')
|
||||
sys.argv.pop(idx)
|
||||
maxadd.append(sys.argv.pop(idx))
|
||||
print(sys.argv[1])
|
||||
if ':' in sys.argv[1]:
|
||||
fname, tally=sys.argv[1].split(':',2)
|
||||
fname, tally=sys.argv[1].split(':',1)
|
||||
if ':' in tally:
|
||||
tally, scale=tally.split(':', 1)
|
||||
scale=float(scale)*SCALING
|
||||
tally=int(tally)
|
||||
else:
|
||||
fname=sys.argv[1]
|
||||
@@ -362,23 +383,27 @@ if __name__=='__main__':
|
||||
meshgrid((y[:-1]+y[1:])/2.,
|
||||
(x[:-1]+x[1:])/2.,
|
||||
(z[:-1]+z[1:])/2.)])[[1,0,2]].T
|
||||
I*=SCALING
|
||||
dI*=SCALING
|
||||
I*=scale
|
||||
dI*=scale
|
||||
fpath=os.path.join('.', 'images', get_name(fname)[0]+'-'+get_name(fname)[1])
|
||||
if tally is not None:
|
||||
fpath+='[%s]'%tally
|
||||
if len(sys.argv)>2:
|
||||
for fi in sys.argv[2:]:
|
||||
scale=SCALING
|
||||
print(fi)
|
||||
if ':' in fi:
|
||||
fi, tally=fi.split(':',2)
|
||||
fi, tally=fi.split(':',1)
|
||||
if ':' in tally:
|
||||
tally, scale=tally.split(':', 1)
|
||||
scale=float(scale)*SCALING
|
||||
tally=int(tally)
|
||||
else:
|
||||
tally=None
|
||||
resi=read_tally(fi, use_tally=tally)
|
||||
xi,yi,zi,Ii,dIi=resi
|
||||
Ii*=SCALING
|
||||
dI*=SCALING
|
||||
Ii*=scale
|
||||
dI*=scale
|
||||
if Ii.shape!=I.shape:
|
||||
print("different shape, interpolate to existing grid...")
|
||||
from scipy.interpolate import interpn
|
||||
@@ -401,16 +426,20 @@ if __name__=='__main__':
|
||||
if tally is not None:
|
||||
fpath+='[%s]'%tally
|
||||
for fi in maxadd:
|
||||
scale=SCALING
|
||||
print("create maximum of other data and", fi)
|
||||
if ':' in fi:
|
||||
fi, tally=fi.split(':',2)
|
||||
fi, tally=fi.split(':',1)
|
||||
if ':' in tally:
|
||||
tally, scale=tally.split(':', 1)
|
||||
scale=float(scale)*SCALING
|
||||
tally=int(tally)
|
||||
else:
|
||||
tally=None
|
||||
resi=read_tally(fi, use_tally=tally)
|
||||
xi,yi,zi,Ii,dIi=resi
|
||||
Ii*=SCALING
|
||||
dI*=SCALING
|
||||
Ii*=scale
|
||||
dI*=scale
|
||||
if Ii.shape!=I.shape:
|
||||
print("different shape, interpolate to existing grid...")
|
||||
from scipy.interpolate import interpn
|
||||
|
||||
Reference in New Issue
Block a user