# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # [phitar](https://jangcom.github.io/phitar/): Testing of the e-type = 2 source subsection
#
# by <NAME>
#
# * Last: 2020-06-06
# * First: 2020-04-27

# ## Introduction
#
# Options for setting source energy distributions will be added to phitar v1.04. One of the options is a Gaussian energy distribution setter that uses the e-type = 2 subsection of PHITS. Before beginning the update, here we examine if the e-type = 2 subsection leads to equivalent results of the e0 command used until phitar v1.03. In addition, whether the e-type = 2 can be used as the e0 command, in which case the former can completely replace the latter, is also tested.

import re
import subprocess
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
from IPython.display import set_matplotlib_formats
from IPython.display import SVG

# +
# I/O
out_path = './out'  # Outbound
phitar_fig_path = '../../../figs'  # Inbound
phits_path = '../phits'

# pandas
flue_header = ['Energy (MeV)', 'Monte Carlo fluence (cm$^{-2}$ source$^{-1}$)']

# MPL
mpl.rcParams.update({'font.sans-serif': 'Arial',
                     'pdf.fonttype': 42})
plt.style.use('bmh')
xlabpad = 15
ylabpad = 15

# IPython.display
set_matplotlib_formats('svg')


# -

def disp_rows_of_int(lines, rows, rows_of_int):
    """Display the rows of interest of a PHITS input file."""
    for i in range(len(lines)):
        for row in rows_of_int:
            if i in rows[row]['range']:
                if i == rows[row]['range'][0]:
                    print(rows[row]['header'])
                print(lines[i])


def save_curr_fig(out_bname, out_path=out_path):
    """Save a figure in multiple formats."""
    out_bname_full = '{}/{}'.format(out_path, out_bname)
    for fmt in ['pdf', 'svg']:
        plt.savefig('{}.{}'.format(out_bname_full, fmt), bbox_inches='tight')
    rtn = subprocess.run('inkscape {}.svg -M {}.emf'.format(out_bname_full, out_bname_full))


def plot_flue(ax, xdata, ydata,
              xmin=30, xmax=40,
              clr='r', mrk='o', ls='None', lab='e0', loc=0):
    """Plot particle fluence against energy."""
    ax.set_xlim(xmin, xmax)
    ax.xaxis.set_minor_locator(AutoMinorLocator(2))
    ax.yaxis.set_minor_locator(AutoMinorLocator(2))
    ax.tick_params(axis='both', which='major', direction='in', length=5)
    ax.tick_params(axis='both', which='minor', direction='in', length=3)
    ax.set_xlabel(flue_header[0], labelpad=xlabpad)
    ax.set_ylabel(flue_header[1], labelpad=ylabpad)
    ax.plot(xdata, ydata, color=clr, marker=mrk, linestyle=ls, label=lab)
    ax.legend(loc=loc)


# ## Comparison using T-Cross tallies

# ### Geometry and tally settings

# Part of a PHITS input file 'geom.inp'

# +
# I/O
f = '{}/geom.inp'.format(phits_path)
fh = open(f)
lines = [v.rstrip() for v in fh]
fh.close()

rows = {
    'params': {
        'header': '',
        'range': [22, 23],
    },
    't_cross': {
        'header': '',
        'range': range(175, 200),
    },
    't_gshow': {
        'header': '',
        'range': [203, 218, 223],
    },
}
rows_of_int = [
    'params',
    't_cross',
    't_gshow',
]

disp_rows_of_int(lines, rows, rows_of_int)
# -

# Run the PHITS input file with icntl = 7 and obtain a G-Show file.

# !phits ../phits/geom.inp

# Obtain a descriptive figure by decorating the G-Show file.

SVG(filename='{}/etype2.svg'.format(phitar_fig_path))

# * Electron beams begin propagation at z = -50 mm
#   toward the +z direction with an incidence angle of 0 deg.
# * The energy spectra of "intact" electron beams are measured
#   by using a T-Cross tally placed at an "ambient" virtual plane at z = -40 mm.
# * The T-Cross tally output is set to be "flux",
#   but is equivalent to "current" because the incidence angle is 0 deg
#   (p. 210, PHITS v3.02 JP man).

# ### Monoenergetic

# Part of a PHITS input file 'e0.inp'

# +
# I/O
f = '{}/e0.inp'.format(phits_path)
fh = open(f)
lines = [v.rstrip() for v in fh]
fh.close()

rows = {
    'params': {
        'header': '',
        'range': [22, 23],
    },
    'source': {
        'header': '',
        'range': range(47, 58),
    },
}
rows_of_int = [
    'params',
    'source',
]

disp_rows_of_int(lines, rows, rows_of_int)
# -

# Run the PHITS input with icntl = 0.

# !phits ../phits/e0.inp

# Check the T-Cross result.

t_cross_bname_full = '{}/e0-cross-eng-intact'.format(phits_path)
t_cross = '{}.ang'.format(t_cross_bname_full)
t_cross_trimmed = '{}_trimmed.ang'.format(t_cross_bname_full)

# Read in the T-Cross result file.
t_cross_fh = open(t_cross)
t_cross_lines = [v for v in t_cross_fh]
t_cross_fh.close()

# Retrieve only the data part and write them to a file.
data_rows = range(251, 1252)
t_cross_data = [t_cross_lines[v] for v in data_rows]
t_cross_data[0] = re.sub('#(.*)', ' \\1', t_cross_data[0])
t_cross_trimmed_fh = open(t_cross_trimmed, 'w')
for line in t_cross_data:
    t_cross_trimmed_fh.write(line)
t_cross_trimmed_fh.close()

df = pd.read_csv(t_cross_trimmed, sep=r'\s+')
df

plot_flue(plt.subplot(), df.loc[:, 'e-upper'], df.loc[:, 'electron'])
save_curr_fig('e0')

# The total electron fluence is

tot_elec_flue = df.electron.sum()
print('{:.3f} cm^-2 source-^1'.format(tot_elec_flue))

# The T-Cross area is

t_cross_area = (t_cross_lines[1263].split())[2]
t_cross_area = float(t_cross_area)
print('{} cm^2'.format(t_cross_area))

# Multiplying the two fluence by the t-cross area gives

tot_elec_flue * t_cross_area

# which is almost equal to 1. This means that the fluence had been calculated as the number of particles passed through the T-Cross area divided by the T-Cross area:

the_tot_elec_flue = 1 / t_cross_area
print('{:.3f} cm^-2 source^-1'.format(the_tot_elec_flue))

# We have just confirmed how the T-Cross tally is calculated (p. 210, PHITS v3.02 Jp man). Now, change the electron energy distribution from monoenergetic to Gaussian-energetic and, see the effect.

# ### Gaussian energy distribution

# Part of a PHITS input file 'etype2.inp'

# +
# I/O
f = '{}/etype2.inp'.format(phits_path)
fh = open(f)
lines = [v.rstrip() for v in fh]
fh.close()

rows = {
    'params': {
        'header': '',
        'range': [22, 23],
    },
    'source': {
        'header': '',
        'range': range(47, 56),
    },
}
rows_of_int = [
    'params',
    'source',
]

disp_rows_of_int(lines, rows, rows_of_int)
# -

# Run the PHITS input with icntl = 0.

# !phits ../phits/etype2.inp

# Check the T-Cross result.

t_cross_bname_full = '{}/etype2-cross-eng-intact'.format(phits_path)
t_cross = '{}.ang'.format(t_cross_bname_full)
t_cross_trimmed = '{}_trimmed.ang'.format(t_cross_bname_full)

# Read in the T-Cross result file.
t_cross_fh = open(t_cross)
t_cross_lines = [v for v in t_cross_fh]
t_cross_fh.close()

# Retrieve only the data part and write them to a file.
data_rows = range(251, 1252)
t_cross_data = [t_cross_lines[v] for v in data_rows]
t_cross_data[0] = re.sub('#(.*)', ' \\1', t_cross_data[0])
t_cross_trimmed_fh = open(t_cross_trimmed, 'w')
for line in t_cross_data:
    t_cross_trimmed_fh.write(line)
t_cross_trimmed_fh.close()

df2 = pd.read_csv(t_cross_trimmed, sep=r'\s+')
df2

plot_flue(plt.subplot(), df2.loc[:, 'e-upper'], df2.loc[:, 'electron'],
          clr='b', mrk='None', ls='-', lab='e-type = 2')
save_curr_fig('etype2')

# The total electron fluence is

tot_elec_flue = df2.electron.sum()
print('{:.3f} cm^-2 source-^1'.format(tot_elec_flue))

# which is almost the same as the monoenergetic one (see [Monoenergetic](#Monoenergetic)).

# ### Additional test: Use of e-type = 2 as e0

# The main commands of the e-type = 2 subsection include (p. 100, PHITS v3.02 JP man)
#
# * eg0: Center of Gaussian distribution (MeV)
# * eg1: FWHM of Gaussian distribution (MeV)
# * eg2: Minimum cutoff for Gaussian distribution (MeV)
# * eg3: Maximum cutoff for Gaussian distribution (MeV)
#
# One can guess that by setting eg1 = 0, the energy distribution would become monoenergetic, in which case the e-type = 2 can completely replace the e0 command.

# Prepare a PHITS input file having e-type = 2 and eg1 = 0.

# Part of a PHITS input file 'etype2_fwhm0.inp'

# +
# I/O
f = '{}/etype2_fwhm0.inp'.format(phits_path)
fh = open(f)
lines = [v.rstrip() for v in fh]
fh.close()

rows = {
    'params': {
        'header': '',
        'range': [22, 23],
    },
    'source': {
        'header': '',
        'range': range(47, 63),
    },
}
rows_of_int = [
    'params',
    'source',
]

disp_rows_of_int(lines, rows, rows_of_int)
# -

# Run the PHITS input with icntl = 0.

# !phits ../phits/etype2_fwhm0.inp

# Check the T-Cross result.

t_cross_bname_full = '{}/etype2_fwhm0-cross-eng-intact'.format(phits_path)
t_cross = '{}.ang'.format(t_cross_bname_full)
t_cross_trimmed = '{}_trimmed.ang'.format(t_cross_bname_full)

# Read in the T-Cross result file.
t_cross_fh = open(t_cross)
t_cross_lines = [v for v in t_cross_fh]
t_cross_fh.close()

# Retrieve only the data part and write them to a file.
data_rows = range(251, 1252)
t_cross_data = [t_cross_lines[v] for v in data_rows]
t_cross_data[0] = re.sub('#(.*)', ' \\1', t_cross_data[0])
t_cross_trimmed_fh = open(t_cross_trimmed, 'w')
for line in t_cross_data:
    t_cross_trimmed_fh.write(line)
t_cross_trimmed_fh.close()

df3 = pd.read_csv(t_cross_trimmed, sep=r'\s+')
df3

fig, axs = plt.subplots(1, 2, figsize=(8, 3.5))
plot_flue(axs[0], df.loc[:, 'e-upper'], df.loc[:, 'electron'], loc=5)
plot_flue(axs[1], df3.loc[:, 'e-upper'], df3.loc[:, 'electron'],
          clr='b', mrk='x', ls='None', lab='e-type = 2 with FWHM = 0', loc=5)
fig.tight_layout(pad=1.5)
save_curr_fig('etype2_as_e0')

# The total electron fluence is

tot_elec_flue = df3.electron.sum()
print('{:.3f} cm^-2 source-^1'.format(tot_elec_flue))

# which is exactly the same as the monoenergetic one (see [Monoenergetic](#Monoenergetic)).

# ## Summary

# * In this test, the e-type = 2 subsection was compared with the e0 command.
# * The obtained electron fluences were found to be equivalent.
# * By setting the FWHM of Gaussian energy distribution (eg1), the e-type = 2 subsection can be used as the e0 command.

StarcoderdataJupyterScriptsDedupFiltered