454 lines
15 KiB
Python
454 lines
15 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 os
|
|
import gzip
|
|
from numpy import *
|
|
from pylab import *
|
|
import matplotlib
|
|
from matplotlib.colors import LogNorm, BoundaryNorm, ListedColormap
|
|
from matplotlib.ticker import LogFormatterMathtext, LogLocator
|
|
from matplotlib.image import imread
|
|
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
|
|
|
#matplotlib.rc('font', size=5)
|
|
|
|
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)
|
|
|
|
BUNKER=False
|
|
FIG_HEIGHT=10.0
|
|
|
|
if BUNKER:
|
|
GEO_PATH="plot_geometry"
|
|
XY_PLANES=[0., -30.]
|
|
XY_GEO=[("pz_0.1_10x10.png", 10.0), ("pz_-0.3_10x10.png", 10.0)]
|
|
XZ_PLANES=[0.]
|
|
XZ_GEO=[("py_0.1_10x10.png", 10.0)]
|
|
YZ_PLANES=[550., 800., 989.9, 1020., 1100.]
|
|
YZ_GEO=[("px_5.6_6x6.png", 6.0), ("px_8.2_6x6.png", 6.0),
|
|
("px_9.9_6x6.png", 6.0), ("px_10.33_6x6.png", 6.0),
|
|
("px_11.15_6x6.png", 6.0)]
|
|
else:
|
|
GEO_PATH="plot_geometry"
|
|
XY_PLANES=[0., -30., 161., 240.]
|
|
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)]
|
|
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),
|
|
("px_23.0_6x6.png", 12.0), ("px_29.0_6x6.png", 12.0),
|
|
("px_35.0_8x8.png", 16.0)]
|
|
|
|
try:
|
|
# image data generated from MCNP plot output for geometry lines
|
|
geometry=load('/home/glavic_a/estia/media/data and simulation/shielding/e2-estia_baseline_v07-0-gb8bc945_lines.npz')
|
|
#geometry=load(r'D:/Documents/Arbeit/OwnCloud/estia/media/data and simulation/shielding\e2-estia_baseline_v07-0-gb8bc945_lines.npz')
|
|
except:
|
|
geometry=None
|
|
|
|
def read_tally(fname, use_tally=None):
|
|
if ';' in fname:
|
|
fname, scale=fname.rsplit(';', 1)
|
|
scale=float(scale)
|
|
else:
|
|
scale=1.0
|
|
if fname.endswith('.gz'):
|
|
txt=gzip.open(fname, 'rt').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*scale,dI*scale
|
|
|
|
def read_mesh_tally(fname, use_tally=None):
|
|
if fname.endswith('.gz'):
|
|
txt=gzip.open(fname, 'rt').readlines()
|
|
else:
|
|
txt=open(fname, 'r').readlines()
|
|
|
|
data=[list(map(float, li.strip().split())) for li in txt[16:]]
|
|
data=array(data).T
|
|
|
|
x=array(txt[10].split(':',1)[1].strip().split(), dtype=float)
|
|
y=array(txt[11].split(':',1)[1].strip().split(), dtype=float)
|
|
z=array(txt[12].split(':',1)[1].strip().split(), dtype=float)
|
|
|
|
nX=len(x)-1
|
|
nY=len(y)-1
|
|
nZ=len(z)-1
|
|
|
|
|
|
_E,_x,_y,_z,I,dI=data.reshape(6, nX, nY, nZ)
|
|
I*=1e-3
|
|
dI*=1e-3
|
|
|
|
return x,y,z,I,dI
|
|
|
|
def plot_xy(outfile, res, zidx=5, geo=None):
|
|
x,y,z,I,dI=res
|
|
xs=x[-1]-x[0]
|
|
ys=y[-1]-y[0]
|
|
I=I[:,:,zidx].T
|
|
fig=figure(figsize=(FIG_HEIGHT*0.9/0.95*xs/ys, FIG_HEIGHT))
|
|
position=[0.05, 0.05, 0.95, 0.9]
|
|
ax=axes(position, aspect='equal')
|
|
p=ax.pcolormesh(x/100., y/100.,
|
|
log10(I+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
X,Y=meshgrid((x[1:]+x[:-1])/200., (y[1:]+y[:-1])/200.)
|
|
ax.contour(X, Y, log10(I+2e-8), [log10(1.5e-6)], colors=['black'])
|
|
xlabel('x [m]')
|
|
ylabel('y [m]')
|
|
title('XY plane at z=%.1fm'%((z[zidx]+z[zidx+1])/200.))
|
|
if geo is not None:
|
|
gf=os.path.join(GEO_PATH, geo[0])
|
|
img=imread(gf)
|
|
img=1.0-img
|
|
x0=(x[0]+x[-1])/200.
|
|
y0=0.
|
|
imgx=linspace(x0-geo[1]/2., x0+geo[1]/2., img.shape[1]+1)
|
|
imgy=linspace(y0+geo[1]/2., y0-geo[1]/2., img.shape[0]+1)
|
|
ax.pcolormesh(imgx, imgy,
|
|
where(img[:,:,0], log10(2e-3), nan),
|
|
norm=Rnorm, cmap=Rcmap)
|
|
ax.pcolormesh(imgx, imgy,
|
|
where(img[:,:,1], log10(2e-3), nan), alpha=0.2,
|
|
norm=Rnorm, cmap=Rcmap)
|
|
elif geometry is not None:
|
|
ax.pcolormesh(geometry['x'], geometry['y'],
|
|
log10(geometry['lines']*0+2e-3),
|
|
norm=Rnorm, cmap=Rcmap
|
|
)
|
|
xlim(x[0]/100, x[-1]/100)
|
|
ylim(y[0]/100, y[-1]/100)
|
|
|
|
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("bottom", size="5%", pad=0.6)
|
|
cbar=fig.colorbar(p, orientation='horizontal', spacing='proportional',
|
|
label='Dose rate [Sv/h]', cax=cax,
|
|
#fraction=0.1, pad=0.1,
|
|
#aspect=50.,
|
|
format="10$^{%i}$")
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300,
|
|
bbox_inches = 'tight', pad_inches = 0)
|
|
|
|
def plot_xz(outfile, res, yidx=5, geo=None):
|
|
x,y,z,I,dI=res
|
|
xs=x[-1]-x[0]
|
|
ys=z[-1]-z[0]
|
|
I=I[:,yidx,:].T
|
|
fig=figure(figsize=(FIG_HEIGHT*0.9/0.95*xs/ys, FIG_HEIGHT))
|
|
position=[0.05, 0.05, 0.95, 0.9]
|
|
ax=axes(position, aspect='equal')
|
|
p=ax.pcolormesh(x/100., z/100.,
|
|
log10(I+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
X,Y=meshgrid((x[1:]+x[:-1])/200., (z[1:]+z[:-1])/200.)
|
|
ax.contour(X, Y, log10(I+2e-8), [log10(1.5e-6)], colors=['black'])
|
|
xlabel('x [m]')
|
|
ylabel('z [m]')
|
|
ypos=(y[yidx]+y[yidx+1])/200.
|
|
title('XZ plane at y=%.1fm'%ypos)
|
|
if geo is not None:
|
|
gf=os.path.join(GEO_PATH, geo[0])
|
|
img=imread(gf)
|
|
img=1.0-img
|
|
x0=(x[0]+x[-1])/200.
|
|
y0=0.
|
|
imgx=linspace(x0-geo[1]/2., x0+geo[1]/2., img.shape[1]+1)
|
|
imgy=linspace(y0+geo[1]/2., y0-geo[1]/2., img.shape[0]+1)
|
|
ax.pcolormesh(imgx, imgy,
|
|
where(img[:,:,0], log10(2e-3), nan),
|
|
norm=Rnorm, cmap=Rcmap)
|
|
ax.pcolormesh(imgx, imgy,
|
|
where(img[:,:,1], log10(2e-3), nan), alpha=0.2,
|
|
norm=Rnorm, cmap=Rcmap)
|
|
|
|
xlim(x[0]/100, x[-1]/100)
|
|
ylim(z[0]/100, z[-1]/100)
|
|
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("bottom", size="5%", pad=0.6)
|
|
cbar=fig.colorbar(p, orientation='horizontal', spacing='proportional',
|
|
label='Dose rate [Sv/h]', cax=cax,
|
|
#fraction=0.1, pad=0.1,
|
|
#aspect=50.,
|
|
format="10$^{%i}$")
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300,
|
|
bbox_inches = 'tight', pad_inches = 0)
|
|
|
|
def plot_yz(outfile, res, xidx=5, geo=None):
|
|
x,y,z,I,dI=res
|
|
xs=y[-1]-y[0]
|
|
ys=z[-1]-z[0]
|
|
I=I[xidx,:,:].T
|
|
fig=figure(figsize=(FIG_HEIGHT*0.9/0.95*xs/ys, FIG_HEIGHT))
|
|
position=[0.05, 0.05, 0.95, 0.9]
|
|
ax=axes(position, aspect='equal')
|
|
p=ax.pcolormesh(y/100., z/100.,
|
|
log10(I+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
X,Y=meshgrid((y[1:]+y[:-1])/200., (z[1:]+z[:-1])/200.)
|
|
ax.contour(X, Y, log10(I+2e-8), [log10(1.5e-6)], colors=['black'])
|
|
xlim(y.max()/100., y.min()/100.)
|
|
xlabel('y [m]')
|
|
ylabel('z [m]')
|
|
title('YZ plane at x=%.1fm'%((x[xidx]+x[xidx+1])/200.))
|
|
if geo is not None:
|
|
gf=os.path.join(GEO_PATH, geo[0])
|
|
img=imread(gf)
|
|
img=1.0-img
|
|
x0=0.
|
|
y0=0.
|
|
imgx=linspace(x0-geo[1]/2., x0+geo[1]/2., img.shape[1]+1)
|
|
imgy=linspace(y0+geo[1]/2., y0-geo[1]/2., img.shape[0]+1)
|
|
ax.pcolormesh(imgx, imgy,
|
|
where(img[:,:,0], log10(2e-3), nan),
|
|
norm=Rnorm, cmap=Rcmap)
|
|
ax.pcolormesh(imgx, imgy,
|
|
where(img[:,:,1], log10(2e-3), nan), alpha=0.2,
|
|
norm=Rnorm, cmap=Rcmap)
|
|
|
|
xlim(y[0]/100, y[-1]/100)
|
|
ylim(z[0]/100, z[-1]/100)
|
|
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("bottom", size="5%", pad=0.6)
|
|
cbar=fig.colorbar(p, orientation='horizontal', spacing='proportional',
|
|
label='Dose rate [Sv/h]', cax=cax,
|
|
#fraction=0.1, pad=0.1,
|
|
#aspect=50.,
|
|
format="10$^{%i}$")
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300,
|
|
bbox_inches = 'tight', pad_inches = 0)
|
|
|
|
def plot_averages(outfile, res):
|
|
# create a graph with averages over some area vs. height
|
|
fig=figure(figsize=(FIG_HEIGHT*1.5, FIG_HEIGHT))
|
|
x,y,z,I,dI=res
|
|
# area above Selene 1 guide, 1m towards Skadi and 2m towards E01
|
|
selene1=1e6*I[x[:-1]<=2300,:,:][:, (y[:-1]>=-100)&(y[:-1]<=200), :].mean(axis=0).mean(axis=0)
|
|
# area between bunker wall and cave, 2m towards Skadi or E01
|
|
guides=1e6*I[(x[:-1]<=3000)&(x[:-1]>=1500),:,:][:, (y[:-1]>=-200)&(y[:-1]<=200), :].mean(axis=0).mean(axis=0)
|
|
# area above experimental cave, to simplify it is a rectangle with xy in [30m:45m] and [-2.5m:6.5m]
|
|
cave=1e6*I[(x[:-1]<=4500)&(x[:-1]>3000),:,:][:, (y[:-1]>=-250)&(y[:-1]<=650), :].mean(axis=0).mean(axis=0)
|
|
|
|
zplot=(z[1:]+z[:-1])/200.
|
|
plot(zplot, selene1, label='Selene 1')
|
|
plot(zplot, guides, label='full guide')
|
|
plot(zplot, cave, label='cave')
|
|
plot(zplot, zplot*0+0.5, label='0.5 µSv/h')
|
|
plot(zplot, zplot*0+1.5, label='1.5 µSv/h')
|
|
plot([2.0, 2.0], [0., 5.], label='reference height')
|
|
legend()
|
|
title('Average dose rate vs. height')
|
|
xlabel('z [m]')
|
|
ylabel('dose rate [µSv/h]')
|
|
ylim(0, 5.0)
|
|
xlim(zplot[0], zplot[-1])
|
|
|
|
if outfile is not None:
|
|
fig.savefig(outfile, transparent=False, dpi=300,
|
|
bbox_inches = 'tight', pad_inches = 0)
|
|
|
|
|
|
def get_name(fname):
|
|
# extract tag-name and version info from filename
|
|
stripped=os.path.basename(fname)[9:].split('.')[0]
|
|
try:
|
|
tag, version=stripped.split('-', 1)
|
|
except ValueError:
|
|
tag=os.path.basename(fname)
|
|
version=''
|
|
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'
|
|
if '-x' in sys.argv:
|
|
idx=sys.argv.index('-x')
|
|
sys.argv.pop(idx)
|
|
output_path=os.path.join('.', 'images', sys.argv.pop(idx))
|
|
else:
|
|
output_path=None
|
|
maxadd=[]
|
|
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)
|
|
tally=int(tally)
|
|
else:
|
|
fname=sys.argv[1]
|
|
tally=None
|
|
res=read_tally(fname, use_tally=tally)
|
|
x,y,z,I,dI=res
|
|
if x.shape==I.shape:
|
|
base_points=None
|
|
else:
|
|
base_points=array([pi.flatten() for pi in
|
|
meshgrid((y[:-1]+y[1:])/2.,
|
|
(x[:-1]+x[1:])/2.,
|
|
(z[:-1]+z[1:])/2.)])[[1,0,2]].T
|
|
I*=SCALING
|
|
dI*=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
|
|
dI*=SCALING
|
|
if Ii.shape!=I.shape:
|
|
print("different shape, interpolate to existing grid...")
|
|
from scipy.interpolate import interpn
|
|
interp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), Ii, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
dinterp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), dIi, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
#print(xi.mean()-x.mean(), yi.mean()-y.mean(),
|
|
#zi.mean()-z.mean(), Ii.max(), interp.max())
|
|
Ii=interp.reshape(*I.shape)
|
|
dIi=dinterp.reshape(*I.shape)
|
|
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
|
|
for fi in maxadd:
|
|
print("create maximum of other data and", 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
|
|
dI*=SCALING
|
|
if Ii.shape!=I.shape:
|
|
print("different shape, interpolate to existing grid...")
|
|
from scipy.interpolate import interpn
|
|
interp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), Ii, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
dinterp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), dIi, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
Ii=interp.reshape(*I.shape)
|
|
dIi=dinterp.reshape(*I.shape)
|
|
dI=where(I>Ii, dI, dIi)
|
|
I=maximum(I, Ii)
|
|
res=(x,y,z,I,dI)
|
|
fpath+='+'+get_name(fi)[1]
|
|
if tally is not None:
|
|
fpath+='[%s]'%tally
|
|
if output_path is not None:
|
|
fpath=output_path
|
|
os.mkdir(fpath)
|
|
fpre=os.path.join(fpath, prefix)
|
|
if x.shape==I.shape:
|
|
x=x[:,0,0]
|
|
y=y[0,:,0]
|
|
z=z[0,0,:]
|
|
plot_averages(fpre+'_average.png', res)
|
|
for i, zi in enumerate(XY_PLANES):
|
|
zidx=where(z>=zi)[0][0]-1
|
|
print(fpre+'_xy_%+05i.png'%zi)
|
|
plot_xy(fpre+'_xy_%+05i.png'%zi, res, zidx=zidx, geo=XY_GEO[i])
|
|
for i, yi in enumerate(XZ_PLANES):
|
|
yidx=where(y>=yi)[0][0]-1
|
|
print(fpre+'_xz_%+05i.png'%yi)
|
|
plot_xz(fpre+'_xz_%+05i.png'%yi, res, yidx=yidx, geo=XZ_GEO[i])
|
|
for i,xi in enumerate(YZ_PLANES):
|
|
xidx=where(x>=xi)[0][0]-1
|
|
print(fpre+'_yz_%+05i.png'%xi)
|
|
plot_yz(fpre+'_yz_%+05i.png'%xi, res, xidx=xidx, geo=YZ_GEO[i])
|