220 lines
7.1 KiB
Python
220 lines
7.1 KiB
Python
#!/usr/bin/env python
|
|
#-*- coding: utf-8 -*-
|
|
"""
|
|
Generate nice graphs from mesh-tally 1 mctal files.
|
|
|
|
At the current stage only simple text processing is done that doesn't work with
|
|
multiple tallies or energy dependent meshes, yet.
|
|
"""
|
|
|
|
import gzip
|
|
from numpy import *
|
|
from pylab import *
|
|
from matplotlib.colors import LogNorm, BoundaryNorm, ListedColormap
|
|
from matplotlib.ticker import LogFormatterMathtext, LogLocator
|
|
|
|
COLORMAP='tab20c'
|
|
SCALING=1.0
|
|
|
|
# generate colormap to emphasize radiation levels properly
|
|
Rcref=['#3182BD', '#6BAED6', '#9ECAE1', '#C6DBEF',
|
|
'#31A354', '#74C476', '#A1D99B', '#C7E9C0',
|
|
'#756BB1', '#9E9AC8', '#BCBDDC', '#DADAEB',
|
|
'#E6550D', '#FD8D3C', '#FDAE6B', '#FDD0A2',
|
|
'#000000', '#111111', '#222222', '#333333',
|
|
'#444444', '#555555', '#666666', '#777777',
|
|
'#888888', '#999999', '#aaaaaa', '#bbbbbb',
|
|
'#cccccc', '#dddddd', '#eeeeee', '#ffffff',
|
|
]
|
|
Rbounds=[
|
|
1e-7, 2.5e-7, 5e-7, 7.5e-7,
|
|
1e-6, 2.5e-6, 5e-6, 7.5e-6,
|
|
1e-5, 2.5e-5, 5e-5, 7.5e-5,
|
|
1e-4, 2.5e-4, 5e-4, 7.5e-4,
|
|
1e-3, 2.5e-3, 5e-3, 7.5e-3,
|
|
1e-2, 2.5e-2, 5e-2, 7.5e-2,
|
|
1e-1, 2.5e-1, 5e-1, 7.5e-1,
|
|
1, 2.5, 5, 7.5, 10
|
|
]
|
|
Rbounds=[log10(bi) for bi in Rbounds]
|
|
Rcmap=ListedColormap(Rcref)
|
|
Rnorm=BoundaryNorm(Rbounds, Rcmap.N)
|
|
|
|
def read_tally(fname, use_tally=None):
|
|
if fname.endswith('.gz'):
|
|
txt=gzip.open(fname, 'r').read()
|
|
else:
|
|
txt=open(fname, 'r').read()
|
|
|
|
tally_ranges=[]
|
|
tidx=0
|
|
while 'tally' in txt[tidx+1:]:
|
|
tidx=txt.index('tally ', tidx+1)
|
|
tnr=int(txt[tidx:].split(None,3)[1])
|
|
if len(tally_ranges)>0:
|
|
tally_ranges[-1][2]=tidx-1
|
|
tally_ranges.append([tnr, tidx, -1])
|
|
if use_tally is None:
|
|
txt=txt[tally_ranges[0][1]:tally_ranges[0][2]]
|
|
else:
|
|
tidx=[ti[0] for ti in tally_ranges].index(use_tally)
|
|
txt=txt[tally_ranges[tidx][1]:tally_ranges[tidx][2]]
|
|
header, data=txt.split('vals')
|
|
|
|
|
|
# get dimensions of data
|
|
fstart=header.index('\nf ')
|
|
hi=header[fstart+3:].splitlines()[0]
|
|
N,nE,nX,nY,nZ=map(int, hi.strip().split())
|
|
|
|
grid=header[fstart+3+len(hi):header.index('\nd ')]
|
|
grid=array(grid.replace('\n', '').strip().split(), dtype=float)
|
|
|
|
x=grid[:nX+1]
|
|
y=grid[nX+1:nX+nY+2]
|
|
z=grid[nX+nY+2:nX+nY+nZ+3]
|
|
|
|
data=array(data.replace('\n','').strip().split(), dtype=float)
|
|
I,dI=data.reshape(nZ, nY, nX, 2).transpose(3,2,1,0)
|
|
|
|
return x,y,z,I,dI
|
|
|
|
def plot_xy(outfile, res, zidx=5):
|
|
x,y,z,I,dI=res
|
|
fig=gcf()
|
|
ax=gca()
|
|
p=ax.pcolormesh(x/100., y/100.,
|
|
log10(I[:,:,zidx].T+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
ax.set_aspect('equal')
|
|
fig.colorbar(p, orientation='horizontal', shrink=0.9, spacing='proportional',
|
|
label='Dose rate [Sv/h]',
|
|
format="10$^{%i}$")
|
|
xlabel('x [m]')
|
|
ylabel('y [m]')
|
|
title('XY plane at z=%.1fm'%((z[zidx]+z[zidx+1])/200.))
|
|
fig.subplots_adjust()
|
|
#fig.savefig(outfile)
|
|
#pc=p.get_facecolors()
|
|
#pc[:,3]=where(I[:,:,zidx].T.flatten()>0, 0.9-0.25*dI[:,:,zidx].T.flatten(), 0.)
|
|
ax.add_patch(Circle((-0.151,0.021), radius=5.50, color='black', fill=False, lw=2))
|
|
ax.add_patch(Circle((-0.151,0.021), radius=11.50, color='black', fill=False, lw=2))
|
|
ax.add_patch(Circle((-0.151,0.021), radius=15.00, color='black', fill=False, lw=2))
|
|
# shielding walls outside bunker
|
|
ax.add_line(Line2D([21.12-0.151, 21.12-0.151], [-1.5,1.5], color='black', lw=2))
|
|
ax.add_line(Line2D([14.845,19.426], [-0.541,-0.713], color='black', lw=2))
|
|
ax.add_line(Line2D([19.426,19.405], [-0.713,-1.248], color='black', lw=2))
|
|
ax.add_line(Line2D([19.405,21.205], [-1.248,-1.460], color='black', lw=2))
|
|
ax.add_line(Line2D([14.723, 22.740], [2.232, 3.084], color='black', lw=2))
|
|
ax.add_line(Line2D([14.801, 20.894], [-1.210, -1.724], color='black', lw=2))
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300)
|
|
|
|
def plot_xz(outfile, res, yidx=5):
|
|
x,y,z,I,dI=res
|
|
fig=gcf()
|
|
ax=gca()
|
|
p=ax.pcolormesh(x/100., z/100.,
|
|
log10(I[:,yidx].T+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
ax.set_aspect('equal')
|
|
fig.colorbar(p, orientation='horizontal', shrink=0.9, spacing='proportional',
|
|
label='Dose rate [Sv/h]',
|
|
format="10$^{%i}$")
|
|
xlabel('x [m]')
|
|
ylabel('z [m]')
|
|
ypos=(y[yidx]+y[yidx+1])/200.
|
|
title('XZ plane at y=%.1fm'%ypos)
|
|
fig.subplots_adjust()
|
|
#fig.savefig(outfile)
|
|
#pc=p.get_facecolors()
|
|
#pc[:,3]=where(I[:,yidx].T.flatten()>0, 0.9-0.25*dI[:,yidx].T.flatten(), 0.)
|
|
ax.plot([sqrt(5.5**2-ypos**2)-0.151, sqrt(5.5**2-ypos**2)-0.151],
|
|
[z.min()/100., z.max()/100.], color='black', lw=2)
|
|
ax.plot([sqrt(11.5**2-ypos**2)-0.151, sqrt(11.5**2-ypos**2)-0.151],
|
|
[z.min()/100., z.max()/100.], color='black', lw=2)
|
|
ax.plot([sqrt(15.0**2-ypos**2)-0.151, sqrt(15.0**2-ypos**2)-0.151],
|
|
[z.min()/100., z.max()/100.], color='black', lw=2)
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300)
|
|
|
|
def plot_yz(outfile, res, xidx=5):
|
|
x,y,z,I,dI=res
|
|
fig=gcf()
|
|
ax=gca()
|
|
p=ax.pcolormesh(y/100., z/100.,
|
|
log10(I[xidx].T+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
ax.set_aspect('equal')
|
|
fig.colorbar(p, orientation='horizontal', shrink=0.9, spacing='proportional',
|
|
label='Dose rate [Sv/h]',
|
|
format="10$^{%i}$")
|
|
xlabel('y [m]')
|
|
ylabel('z [m]')
|
|
title('YZ plane at x=%.1fm'%((x[xidx]+x[xidx+1])/200.))
|
|
fig.subplots_adjust()
|
|
#fig.savefig(outfile)
|
|
#pc=p.get_facecolors()
|
|
#pc[:,3]=where(I[xidx].T.flatten()>0, 0.9-0.25*dI[xidx].T.flatten(), 0.)
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300)
|
|
|
|
def get_name(fname):
|
|
# extract tag-name and version info from filename
|
|
stripped=os.path.basename(fname)[9:].split('.')[0]
|
|
tag, version=stripped.split('-', 1)
|
|
return tag, version
|
|
|
|
|
|
if __name__=='__main__':
|
|
import sys, os
|
|
if '-o' in sys.argv:
|
|
idx=sys.argv.index('-o')
|
|
sys.argv.pop(idx)
|
|
prefix=sys.argv.pop(idx)
|
|
else:
|
|
prefix='estia_dosemap'
|
|
print sys.argv[1]
|
|
if ':' in sys.argv[1]:
|
|
fname, tally=sys.argv[1].split(':',2)
|
|
tally=int(tally)
|
|
else:
|
|
fname=sys.argv[1]
|
|
tally=None
|
|
res=read_tally(fname, use_tally=tally)
|
|
x,y,z,I,dI=res
|
|
I*=SCALING
|
|
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:]:
|
|
print fi
|
|
if ':' in fi:
|
|
fi, tally=fi.split(':',2)
|
|
tally=int(tally)
|
|
else:
|
|
tally=None
|
|
resi=read_tally(fi, use_tally=tally)
|
|
xi,yi,zi,Ii,dIi=resi
|
|
Ii*=SCALING
|
|
I+=Ii
|
|
dI=sqrt(dI**2+dIi**2)
|
|
res=(x,y,z,I,dI)
|
|
fpath+='+'+get_name(fi)[1]
|
|
if tally is not None:
|
|
fpath+='[%s]'%tally
|
|
os.mkdir(fpath)
|
|
fpre=os.path.join(fpath, prefix)
|
|
xi, yi, zi=res[3].shape
|
|
for zidx in range(zi):
|
|
print fpre+'_xy_%02i.png'%zidx
|
|
figure(figsize=(12,8))
|
|
plot_xy(fpre+'_xy_%02i.png'%zidx, res, zidx=zidx)
|
|
for yidx in range(yi):
|
|
print fpre+'_xz_%02i.png'%yidx
|
|
figure(figsize=(12,8))
|
|
plot_xz(fpre+'_xz_%02i.png'%yidx, res, yidx=yidx)
|
|
for xidx in range(xi):
|
|
print fpre+'_yz_%02i.png'%xidx
|
|
figure(figsize=(8,10))
|
|
plot_yz(fpre+'_yz_%02i.png'%xidx, res, xidx=xidx)
|
|
|