diff --git a/analyze_materials.py b/analyze_materials.py new file mode 100644 index 0000000..eb87a0d --- /dev/null +++ b/analyze_materials.py @@ -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] + diff --git a/plot_tally.py b/plot_tally.py index dc8f593..dd69c70 100644 --- a/plot_tally.py +++ b/plot_tally.py @@ -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