#!/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])