The correct solution? [coding]

Generally when programming, you can ask yourself; is this the most convenient way of doing this? Sometimes you don’t know the answer and often you won’t know it until you find the more suitable coding trick or language. I will here present my insights from learning how to code in R and in which way it simplified my coding.

The purpose of the program

I will here continue on coding for nuclear cross section analysis. Generally when you work in the area of experimental and theoretical data you want to know, in quantifiable terms, what kind of errors there are between the experimental EXFOR data and your own theoretical calculations. My goal with the program in question was to quantify the correctness of my theoretical calculations.

The statistical method I chose for this purpose was Chi-squared per degrees of freedom. Why? In part due to the suitable statistics in the method and in part due to the wide application of the method in the area of nuclear physics. The experimental data is discrete, and this also applies to the data from the nuclear reaction codes EMPIRE and TALYS. The only key requirement in comparing two data sets is that the theoretically calculated cross sections should have an energy set with the relevant experimental energy as a subset.

The codes chosen

So how do you program this? There are several solutions for this problem, but the solutions I will present here is a program in C, and a script in Bash and R. The main difference in the two methods are that C lacks all inbuilt statistical functions, while R is designed for the purpose of statistical computing.

The code: C

/* Chi-squared calculation program
Implementation of the Chi-Square Distribution & Incomplete Gamma Function in C
Written and modified by Alex Björkholm, 2015*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include <stdbool.h>
#include <unistd.h>
#include <dirent.h>
#include <sys/stat.h>

static long double log_igf (long double K, long double X);
static long double KM (long double K, long double X);
long double log_gamma (long double K);
static double chisq_prob (int DOF, double chisq);
static char *findpath(char *cwddir,char *find);
bool match(char *first, char * second);


int main()
{
printf("\n\t################### CHI-SQUARED, DATA #####################");
printf("\n\tProgram for calculation of chi-squared of theoretical data");
printf("\n\tfiles in comparison with (multiple) experimental data files\n");


/* Defining working directory */
char cwddir[1024];
if (getcwd(cwddir, sizeof(cwddir)) != NULL) {
printf("\tCurrent working dir: %s\n", cwddir);
}

/* Parameter definition */
int i,j,f;
int n,m;
char end[1024];
char prminp[20][20];
double prmbase[2][800];      //infoheader, energyrows
double prmref[20][10][800];  //file, infoheader, energyrows
double prmcomp[20][800];     //infoheader, energyrows

char find[20];
char *file_path;
char filebase_path[20][1024];
char fileref_path[20][1024];
int fileref_no;




/* Reading file names from terminal */
file_path = NULL;
find[0]='\0';
do {
printf("\tEnter the theoretical data search string with or without search keys '*'/'?': ");
scanf("%s",find);
file_path = findpath(cwddir,find);

}while ( file_path == NULL );

strcpy(filebase_path[0], file_path);
printf("\tFile path set: %s\n", filebase_path[0]);


printf("\n\tEnter the number of experimental files to reference your theoretical data to: ");
scanf("%i",&fileref_no);

i = 0;
while ( i < fileref_no ) {
file_path = NULL;
find[0]='\0';
do {
printf("\n\tEnter the experimental data search string with or without search keys '*'/'?': ");
scanf("%s",find);
file_path = findpath(cwddir,find);

}while ( file_path == NULL );

strcpy(fileref_path[i],file_path);
i++;
}


printf("\n\tFile base path:\t\t %s", filebase_path[0]);
for ( j = 0; j < fileref_no ; j++ ) {
printf("\n\tFile reference path 0%i:\t %s",j,fileref_path[j]);
}



/* Getting theoretical data from file into array */
/* Clearing arrays */
for ( i = 0; i < 20; i++) {
prminp[i][0]='\0';
}

for ( i = 0; i < 2; i++) {
for ( j = 0; j < 800; j++) {
prmbase[i][j]=0;
}
}

/* Opening filebase_name */
FILE *fpbase;
fpbase = fopen(filebase_path[0],"r"); // read mode
if( fpbase == NULL ) {
perror("\n\tError while opening the file.\n\n");
exit(1);
}

j=0;
/* Extracting the relevant data (first and second column) from filebase_name and storing it into arrays */
while( fgets(end, 1024, fpbase) != NULL ) {
fscanf( fpbase, "%s %s", prminp[0], prminp[1]);
if (( strlen(prminp[0]) != 0 ) && ( prminp[0][0] != '#' )) {
prmbase[0][j] = atof(prminp[0]);
prmbase[1][j] = atof(prminp[1]);
j++;
}
for ( i = 0; i < 20; i++) {
prminp[i][0]='\0';
}
}
fclose(fpbase);


/* Printing base array */
printf("\n\tPrinting relevant info from file: %s\n",filebase_path[0]);
for ( j = 0; j < 800; j++) {
if ( prmbase[0][j] != 0 ) {
for ( i = 0; i < 2; i++) {
printf("\t%0.6E, %i %i\t", prmbase[i][j], i ,j);
}
printf("\n");
}
}



/* Getting experimental data from files into array */
/* Clearing arrays */
for ( f = 0; f < 20; f++) {
for ( i = 0; i < 10; i++) {
for ( j = 0; j < 800; j++) {
prmref[f][i][j]='\0';
}
}
}

/* Opening fileref_name */
for ( n = 0; n < fileref_no; n++) {

FILE *fpref;
fpref = fopen(fileref_path[n],"r"); // read mode
if( fpref == NULL ) {
perror("\tError while opening the file.\n\n");
exit(1);
}

/* Extracting the relevant data (first and second column) from fileref_name and storing it into arrays */
for ( i = 0; i < 20; i++) {
prminp[i][0]='\0';
}

j = 0;
while( fgets(end, 1024, fpref) != NULL ) {
fscanf( fpref, "%s %s %s %s", prminp[0], prminp[1], prminp[2], prminp[3]);
if (( strlen(prminp[0]) != 0 ) && ( prminp[0][0] != '#' )) {
prmref[n][0][j] = atof(prminp[0]); // Energy data
prmref[n][1][j] = atof(prminp[1]); // Cross section data
j += 1;
}
for ( i = 0; i < 20; i++) {
prminp[i][0]='\0';
}

}

fclose(fpref);
}


/* Printing reference array */
for ( f = 0; f < fileref_no; f++) {
printf("\n\tPrinting relevant info from file: %s\n",fileref_path[f]);
for ( j = 0; j < 80; j++) {
if ( prmref[f][0][j] != 0 ) {
for ( i = 0; i < 2; i++) {
printf("\t%0.6E\t", prmref[f][i][j]);
}
printf("\n");
}
}
}



/* Setting data from parameter arrays with corresponding energy into same array */
for ( i = 0; i < 20; i++) {
for ( j = 0; j < 800; j++) {
prmcomp[i][j] = '\0';
}
}

n = 0;
for ( m = 0 ; m < 800; m++) {
if (( prmbase[0][m] != 0 ) && ( prmbase[0][m] != '\0')) { // ignoring null energies in theoretical file
for ( f = 0 ; f < fileref_no; f++) {
for ( j = 0 ; j < 80; j++) {
if ( prmbase[0][m] == prmref[f][0][j] ) {
prmcomp[0][n] = prmbase[0][m];
prmcomp[2][n] = prmbase[1][m];
prmcomp[f+4][n] = prmref[f][1][j];
}
}
}
}
n++;
}



/* Calculating expected value of all experimental values for each element */
for ( j = 0; j < 800; j++) {
double N = 0;
for ( i = 4; i < 20; i++) {
if (( prmcomp[i][j] != '\0' ) && ( prmcomp[i][j] != 0 )){
prmcomp[3][j] = prmcomp[3][j] + prmcomp[i][j];
N = N + 1;
}
}
prmcomp[3][j] = prmcomp[3][j]/N;
}



/* Calculating chi-squared for each element */
double N = 0;
for ( i = 0; i < 800; i++) {
if (( prmcomp[0][i] != '\0' ) && ( prmcomp[3][i] != 0 )) {
prmcomp[1][i] = pow((prmcomp[2][i] - prmcomp[3][i]),2) / prmcomp[3][i];
N++;
}
}

/* Calculating sum of chi-squared per DOF for elements */
double chisq = 0;
for ( i = 0; i < 800; i++) {
chisq = chisq + prmcomp[1][i];
}

double DOF = N - 1.0;
double chisqdof = chisq / DOF;

printf("\n\tThe correlating, relevant, parameters:");
printf("\n\tchisq %0.4E\t chisqdof %0.4E\t DOF %0.2f\n\n", chisq, chisqdof, DOF);



/* Opening filecomp_name */
FILE *fpcomp;
fpcomp = fopen("chisq-calculate-output.tot","wb+"); // write mode, creating new file if needed
if( fpcomp == NULL ) {
perror("\tError while opening the file.\n\n");
exit(1);
}

/* Printing data from comparison parameter array */
fprintf(fpcomp, "# Theoretical input file:\n# %s\n#\n", filebase_path[0]);
fprintf(fpcomp, "# Experimental input file(s):\n");
printf("\tTheoretical input file:\n\t%s\n\n", filebase_path[0]);
printf("\tExperimental input file(s):\n");
for ( i = 0 ; i < fileref_no ; i++) {
fprintf(fpcomp, "# %s\n", fileref_path[i]);
printf("\t%s\n", fileref_path[i]);
}

fprintf(fpcomp, "#\n# Energy [MeV]\tchi-squared\txs(theor.) [mb]\txs(exp. aver.) [mb]\txs(exp. individual) [mb]\n");
printf("\n\tEnergy [MeV]\t\tchi-squared\t\txs(theor.) [mb]\t\txs(exp. aver.) [mb]\txs(exp. ind.) [mb]\n");
for ( j = 0 ; j < 800 ; j++) {
if (( prmcomp[0][j] != '\0' ) && ( prmcomp[0][j] != 0 )){
for ( i = 0 ; i < fileref_no + 4 ; i++) {
printf("\t%0.6E\t", prmcomp[i][j]);
fprintf(fpcomp, "%0.6E\t", prmcomp[i][j]);
}

printf("\n");
fprintf(fpcomp,"\n");
}
j++;
}


printf("\n\tDegrees of freedome:  %0.2f\n", DOF);
printf("\tChi-squared value:  %0.2f\n", chisq);
printf("\tProbability value:  %0.15f\n\n", chisq_prob(DOF, chisq));
fprintf(fpcomp,"# Degrees of freedome:  %0.2f\n", DOF);
fprintf(fpcomp,"# Chi-squared value:  %0.2f\n", chisq);
fprintf(fpcomp,"# Probability value:  %0.15f\n\n", chisq_prob(DOF, chisq));

fclose(fpcomp);

return 0;

}


static char *findpath(char *cwddir,char *find)
{
char *pathout;
char *finddir;

char temp[50];
char temp01[50];
char temp02[50];

char path[1024];
char wdir[1024];

strcpy(temp, find);
strcpy(wdir, cwddir);
char target = '/';
int n = 0;


DIR *dp;
struct dirent *entry;
struct stat statbuf;

if( (dp = opendir(".")) == NULL ) {
fprintf(stderr,"\tCannot open directory: %s\n", wdir);
return 0;
}

/* checking if 'find' contains '/' */
while ( strchr(find, target) != NULL ) {

/* Splitting 'find' in two parts dependent on '/'.
Second string modified by removing first char with +1. */
finddir = strtok(find, "/");
find = (strrchr(temp, '/') + 1);

while( (entry = readdir(dp)) != NULL ) {
lstat(entry->d_name,&statbuf);

if (( S_ISDIR(statbuf.st_mode) ) && ( match(finddir, entry->d_name) == true )) {
strcpy(temp01, entry->d_name);
n++;
}
}
closedir(dp);

if ( n == 1 ) {
sprintf(wdir,"%s%s%s",wdir,"/",temp01);
chdir(wdir);
n = 0;
}
else {
printf("\n\tERROR: %i directories found with the string '%s'\n\n",n,finddir);
}

}


if( (dp = opendir(".")) == NULL ) {
fprintf(stderr,"\tCannot open directory: %s\n", wdir);
return 0;
}

while( (entry = readdir(dp)) != NULL ) {
lstat(entry->d_name,&statbuf);
sprintf(path,"%s%s%s",wdir,"/",entry->d_name);

/* If regular file and contains "find" string */
if(( S_ISREG(statbuf.st_mode) ) && ( match(find, entry->d_name) == true )) {
printf("\tAvailable file: %s\n",path);
strcpy(temp02, entry->d_name);
n++;
}
}
closedir(dp);
chdir(cwddir);

if ( n == 1 ) {
sprintf(path,"%s%s%s",wdir,"/",temp02);
pathout = path;
return pathout;
}
else {
printf("\n\tERROR: %i files found with the string '%s'\n\n",n,find);
return NULL;
}

}


/* The main function that checks if two given strings match.
The first string may contain wildcard characters */
bool match(char *first, char *second)
{
// If we reach at the end of both strings, we are done
if (*first == '\0' && *second == '\0')
return true;

// Make sure that the characters after '*' are present in second string.
// This function assumes that the first string will not contain two
// consecutive '*'
if (*first == '*' && *(first+1) != '\0' && *second == '\0')
return false;

// If the first string contains '?', or current characters of both
// strings match
if (*first == '?' || *first == *second)
return match(first+1, second+1);

// If there is *, then there are two possibilities
// a) We consider current character of second string
// b) We ignore current character of second string.
if (*first == '*')
return match(first+1, second) || match(first, second+1);
return false;
}



/* Wikipedia - Chi-squared Distribution -> Cumulative distribution function
http://en.wikipedia.org/wiki/Chi-squared_distribution#Cumulative_distribution_function */
static double chisq_prob(int DOF, double chisq)
{
if (chisq < 0 || DOF < 1)
{
return 0.0;
}

double K = ((double)DOF) * 0.5;
double X = chisq * 0.5;
if(DOF == 2) // special case for Probability function
{
return exp(-1.0 * X);
}

long double PValue;
long double ln_CDF;
ln_CDF = log_igf(K, X)- log_gamma(K);
PValue = 1.0 - expl(ln_CDF);

return (double)PValue;

}


/* Returns the Natural Logarithm of the Incomplete Gamma Function
The logaritmic expression allows more accurate calculations, but the number of iterations are larger.
Wikipedia - Incomplete Gamma Function -> Evaluation formulae -> Connection with Kummer's confluent hypergeometric function
https://en.wikipedia.org/wiki/Incomplete_gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function
*/
static long double log_igf(long double K, long double X)
{
if(X < 0.0)
{
return 0.0;
}

long double Sc;
long double C;
Sc = (K * logl(X)) - X - logl(K);
C = KM(K, X);

return logl(C) + Sc;
}


/* Kummer's confluent hypergeometric function for Incomplete Gamma function
Wikipedia - Incomplete Gamma Function -> Evaluation formulae -> Connection with Kummer's confluent hypergeometric function
https://en.wikipedia.org/wiki/Incomplete_gamma_function#Connection_with_Kummer.27s_confluent_hypergeometric_function
*/
static long double KM(long double K, long double X)
{
long double Sum = 1.0;
long double Nom = 1.0;
long double Denom = 1.0;
int i = 0;

for(i = 0; i < 1000; i++) // Loops for 1000 iterations
{
Nom *= X;
K++;
Denom *= K;
Sum += (Nom / Denom);
}

return Sum;
}


/*
Returns the Natural Logarithm of the Gamma Function, based on Spouge's approximation (modification of Stirling).
The logaritmic expression allows more accurate calculations, but the number of iterations are larger.
Wikipedia - Spouge's approximation: https://en.wikipedia.org/wiki/Spouge's_approximation
*/
long double log_gamma(long double K)
{
/* A is the level of accuracy of the calculation. Spouge's approximation: the desired level of precision can only be reached with
an extra precision available so that the desired level can be built. */
double A = 15;
int k = 1;
const long double SQRT2PI = sqrtl(atanl(1.0) * 8.0);

long double Sc;

Sc = ((K + 0.5)*logl(K + A)) - (K + A) - logl(K);

long double Sum = SQRT2PI;
long double Ck = 0.0;
long double F = 1.0;

for(k = 1; k < A; k++)
{
Ck = powl(-1,k-1) * powl(A - k, k - 0.5) * expl(A - k);
Sum += (Ck / (F * (K + k)));
F *= k;
}

return logl(Sum) + Sc;
}





As a purely C code this program somewhat lacks the dynamic nature that I find convenient in Bash. It still functions and does the job, but the programming isn’t really versatile in its nature. You couldn’t really use this code for any other purpose that similar to the current application. Anything else would have to be coded into place.

The code: Bash & R

#!/bin/bash
#
# This script calculates the chi-squared value of a theoretical file compared to experimental data
# by using R statistical software
# Author: Alex Björkholm 2016


# Manual of input output
if [ "$1" == "man" ]; then
echo -e "This script calculates the chi-squared value of a theoretical file compared to experimental data by using R statistical software.
Input options:
1: input file search string, theoretical data (single file)
2: input file search string, experimental data (single/multiple files)
3: (no)/energy
experimental energy sort (optional) active if input 'energy'\n"

else

input1_theoretical=$1
input1_experimental=$2
energy=$3

declare -a files_theoretical
declare -a files_experimental
declare -a input_theoretical
declare -a input_experimental

IFS=$'\n'

#Finding input file, theoretical
files_theoretical+=($(find . -maxdepth 1 -type f -iwholename \*$input1_theoretical\* | sort -u))

#Finding input file, experimental
files_experimental+=($(find . -maxdepth 1 -type f -iwholename \*$input1_experimental\* | sort -u))


#Declare files
echo -e "\nInput files, theoretical"
for each in ${files_theoretical[@]}; do echo -e "\t$each"; done

echo -e "\nInput files, experimental"
for each in ${files_experimental[@]}; do echo -e "\t$each"; done

if [ ! ${#files_theoretical[@]} == 1 ]; then
echo -e "\nNumber of theoretical files is not singular (no. of files = ${#files_theoretical[@]})\n"
exit 2
fi



for each in ${files_theoretical[@]}; do
# Remove irrelevant headers, leaving title
sed -e '1{/\#.*/d}' -e '2{/\#.*/d}' -e 's/\#//g' -e 's/type//g' $each > temp.tot
mv temp.tot $each
done

echo -e "\nNumber of experimental files: ${#files_experimental[@]}"
for each in ${files_experimental[@]}; do
# Check files for theoretical data and problems
if [ ! $(grep -i "empire\|talys\|tendl" $each | wc -l) == 0 ]; then
echo -e "\nWarning: file $each contains string empire\|talys\|tendl"
# Exit possibility for script
echo -n "Exit script? Yes/No? >"
read input
if [[ "${input,,}" = y* ]]; then
echo "\nExiting script\n\n"; exit 2
fi
fi

# Remove irrelevant headers, leaving title
sed -e '1{/\#.*/d}' -e '2{/\#.*/d}' -e 's/\#//g' -e 's/type//g' $each > temp.tot
mv temp.tot $each
done



# Setting input files for R
for each in ${files_theoretical[@]}; do
each=$(basename $each)
input_theoretical+=($(echo "\"$each\"",))
done
theoretical=$(echo ${input_theoretical[@]} | sed s'/.$//')

# Setting input files for R
for each in ${files_experimental[@]}; do
each=$(basename $each)
input_experimental+=($(echo "\"$each\"",))
done
experimental=$(echo ${input_experimental[@]} | sed s'/.$//')



cat <<EOF > R-script.R
writeLines("\n\nExecuting R statistical software code")


experimental <- c(${experimental[@]})
noexperimental <- ${#files_experimental[@]}
theoretical <- c(${theoretical[@]})
notheoretical <- ${#files_theoretical[@]}


# Importing theoretical data from file
Rtheoretical <- read.table(theoretical, header=TRUE)
# Setting data labels for theoretical data
colnames(Rtheoretical)[1]<-paste("Energy")
colnames(Rtheoretical)[2]<-paste("xstheor")


Rexperimental<-0
# Importing experimental data from file
for (val in experimental) {
temp <- read.table(val, header=TRUE)
temp <- temp[,c("x", "y")]
# Setting data labels for experimental data
colnames(temp)[1]<-paste("Energy")
colnames(temp)[2]<-paste("xsexper")
if( ! is.data.frame(Rexperimental) ){
#Storing first parameters in Rexperimental when no data present
Rexperimental<-temp }
else {
# Merging the two data frames of experimental data
# Note: duplicate values will dissapear. This need to be resolved
# if one wants a weighted data analysis
Rexperimental<-merge(Rexperimental, temp, all = T)}
}

# Producing data table with single energy input
# Taking means of cross section for identical energy values
Rexperimental<-aggregate( xsexper ~ Energy, Rexperimental, mean )


# Merging the matching values from Rtheoretical and Rexperimental
# Getting matching values from Rtheoretical
t1 <- Rtheoretical[ Rtheoretical\$Energy %in% Rexperimental\$Energy,]
# Getting matching values from Rexperimental
t2 <- Rexperimental[ Rexperimental\$Energy %in% Rtheoretical\$Energy,]


# Merging energies and cross sections into single data table
xs.data.table <- cbind(t1, t2)
# Writing new row numbering for data set
rownames(xs.data.table)<-NULL


# Merging cross sections into single data table for chisq.test
chisq.data.table <- cbind(t1\$xstheor, t2\$xsexper)
# Setting header names for chisq.data.table
colnames(chisq.data.table)<-c("xstheor", "xsexper")
# Writing new row numbering for data set
rownames(chisq.data.table)<-c(1:nrow(chisq.data.table))


# Printing resulting table
writeLines("\nCross section data")
print(xs.data.table)
writeLines("\nChi-squared input data")
print(chisq.data.table)

# Running chi squared test
writeLines("\n\nChi-squared data")
chisq.test(chisq.data.table)

q()
EOF

# Executing R-script
Rscript R-script.R

unset IFS
fi

echo -e "\nScript finished\n\n"

Note that the integration of Bash in the script isn’t essential and the reasoning behind this solution was mainly to enable the strong file search capabilities of find in Bash. The script could have been written entirely in R, but I chose not to at the moment of coding. The script therefore also displays the capability to execute scripting via Bash. The execution is however somewhat clumsy when it involves the manufacturing of a R script followed by the script’s execution.

The difference?

As seen above both programs do more or less the same thing. The main difference is that R comes with already inbuilt functions for handling the statistics, while the C program has the statistics coded into existence.

One can of course ask which program is better. I prefer the Bash & R solution, mainly due to the versatility of the code. Adding new features isn’t difficult, if I would need any, and there is inbuilt support in R that makes the coding simpler.

How does it work?

Yeah, so how does the code work? Mainly you need the experimental data sets as well as the theoretical data set. This can be acquired via the data plot script I have previously written about, with the data output set to the table option.  We need to remember that our focus are finding how good the theoretical calculations are compared to the experimental values. So before the program is run a theoretical calculation should be preformed, whereby we need to get theoretical values (here nuclear reaction cross sections) for all the existing experimental reference points (in this case energies). After that the scripts/programs would need to be executed in the directory of the experimental and theoretical files.

The C program, when executed from the terminal command line, will ask for the theoretical file name (actually it searches for the name based on a search string, since I’m to lazy to always copy paste the full name) and then ask for the number of experimental files to compare to, where after one needs to to provide the file search strings for those files. After that the script will make the necessary checks when it executes itself where after the result is printed out in the terminal.

The Bash and R program works in a similar way, but here we can enter all the necessary data straight into the terminal when running the program. So everything is a lot smoother in execution.

Does the code work?

Well, off course. I hope that I have commented my code to the degree that it is comprehensible for anyone, but that is usually ofter quite hard to acheive. For my needs this works perfectly. It might need some readthrough before implementation so that you understand its working principles in order to get the correct results. But it should be very straight forward! Anyway, I hope this inspires.

A final note

Usually your files aren’t that clean form extra formatting. Especially if they have been produced via the table option in gnuplot. My easy fix for this is a simple bash script – no surprise there, due to my love for Bash – which cleans the table file from all the unnecessary formatting and formats it in the correct way for the Chi squared scripts (C or Bash). This code can be found here below.

#!/bin/bash
#
# This script splits and refines output table from gnuplot for further analysis
# Author: Alex Björkholm 2016


# Manual of input output
if [ "$1" == "man" ]; then
echo -e "This script splits and refines output table from gnuplot for further analysis
Input options:
1: input file
2: output name string
3: (no)/energy
experimental energy sort (optional) active if input 'energy'\n"

else

inputfile=$1
outputname=$2
energy=$3

declare -a files
IFS=$'\n'

# Using sed to remove every 'i' or 'o' at lines without '#'
# Using sed to remove lines with 'NaN'
# Using sed to remove the first line if it is empty
sed -e '/^\#/!s/i//' -e '/^\#/!s/o//' -e '/NaN.*/d' -e '1{/^$/d}' $inputfile | cat -s > temp.out

# Counting number of data sets based on the numer of titles
#N=$(($(grep -o '\#' temp.out | wc -l)/3 - 2))
N=$(($(grep -o 'title' temp.out | wc -l) - 2))

# Splitting file at every empty row
csplit -sf $outputname -b '%d.tot' temp.out '/^$/' {$N}
rm temp.out


# Removing first empty row in every file
files+=($(find . -maxdepth 1 -type f -iwholename \*$outputname\*.tot | sort -u))
for each in ${files[@]}; do
# Removing first blank row from file
sed '1{/^$/d}' $each > temp.out
mv temp.out $each
# Removing first blank space at line
sed -i 's/ //' $each
done


if [[ ! -z ${energy// } ]]; then
if [ $energy == "energy" ]; then
for each in ${files[@]}; do
# Getting all energy data points not matching empire, talys or tendl
if [ $(grep -i "empire\|talys\|tendl" $each | wc -l) == 0 ]; then
# Extracting the first number/word
cat $each | cut -d " " -f1 >> temp.out
fi
done

# Removing comment lines and empty lines
sed -e '/\#.*/d' -e '/^$/d' temp.out > energy.tot
# Sorting data file
sort -n -u energy.tot -o energy.tot
rm temp.out
fi
fi

unset IFS
fi

Exporting 3D objects from Octave

object_with_colour

Colormap/texture applied to Wavefront object format

In simplified terms; GNU Octave is the major open license version of MATLAB. with a license via the GNU General Public License. As always with non-commercial versions of programs we can see a certain lack of fine-tuning in the software as well as some missing functions from the commercial, licensed, version. I will here show some key points in working with Octave in relation to modelling of 3D data and saving it in a 3D data format. This wasn’t at all straight forward when I set out to export 3D objects.

The basics – axis

First of all: when working with 3D export from Octave you need to know what kind of axis dimensions you are working with in your plot. The question is relevant in whether your 3D output actually will relate to the 3D surface/mesh that you have created. You can test this in Octave/MATLAB with this simple command:

axis equal;

If you, after the execution of the command above, retain a 3D figure that looks good, then you can proceed with exporting your 3D image. Otherwise you need to fix your X, Y, and Z axis so that they are represented inside a quadratic cube. I solved this issue in my 3D to 3D surface plots by automatically resizing the Z axis, depending on the values in the X and Y arrays. I was before this correction to my 3D plots only able to export a very small segment of the X-Y plane.

3D formats?

So the next step is actually knowing what/where you are going to use the 3D image. Does it need to be in colour? Do you want to display it via JavaScript or just play around with it in a 3D viewer? You can find a list of the more common 3D object formats on Wikipedia, but here are some things to remember for the novice:

  • Not all formats are all inclusive with structures and colours. If you want a colourised 3D object you will need a format that actually supports it in some way.
  • Colour support is often supplied via JPEG or PNG images. This is defined as texture mapping in 3D modelling. The process of advanced texture mapping can be cumbersome in 3D modelling, but it is often very easy to apply a colour or colormap.
  • Note that even though colour is supported in the format, it doesn’t mean that it is all saved in one file. For example; Wavefront object format supports colours and textures, but the colour/texture is saved in separate files.

Octave export scripts for 3D objects

If you don’t need colour for your 3D object, then you can export your plot more or less to any format. The easiest format is the Wavefront object file format. To do this you will only need to run the following code from your Octave terminal:


% Converting 2D matrix parameters to 1D
x = (X')(:)';
y = (Y')(:)';
z = (C')(:)'/k;

% Creating an output file
fullfilename = 'image.obj'
fid = fopen(fullfilename,'w');

% Printing the data points in the 3D mesh
for i=1:size(x,2)
  fprintf(fid,'v %.2f %.2f %.4f\n',y(1,i),x(1,i),z(1,i)/255);
end

fprintf(fid,'f 1 2 3\n');
fclose(fid);

This however does not support colours, which you can see if you open up the exported object in (for example) MeshLab and try to attach a texture to it. Similar results are also attained by scripts like surf2stl.m. Why is it not possible to attach a texture? For me this is usually due to the surface basis on the exported 3D object, since I create my 3D plots via the surf()command. So essentially it is based on a surface, and not a mesh, and we need a mesh to be able to apply a texture.

You can however apply the method above generally in 3D object exports. The basics are the same, and the main challenge really is to identifying the structure of the desired file format. Our solution with the missing mesh format? Well, we need to apply a mesh based exportation. Here the credits goes to Anders Sandberg and his excellent page with multiple scripts for MATLAB output to the object file format. From there you can download the script you need for your specific purpose, but the one I use is saveobjmesh.m. The code is also available here below.

function saveobjmesh(name,x,y,z,nx,ny,nz)
% SAVEOBJMESH Save a x,y,z mesh as a Wavefront/Alias Obj file
% SAVEOBJMESH(fname,x,y,z,nx,ny,nz)
%     Saves the mesh to the file named in the string fname
%     x,y,z are equally sized matrices with coordinates.
%     nx,ny,nz are normal directions (optional)

  normals=1;
  if (nargin&amp;amp;lt;5) normals=0; end
  l=size(x,1); h=size(x,2);  

  n=zeros(l,h);

  fid=fopen(name,'w');
  nn=1;
  for i=1:l
    for j=1:h
      n(i,j)=nn;
      fprintf(fid, 'v %f %f %f\n',x(i,j),y(i,j),z(i,j));
      fprintf(fid, 'vt %f %f\n',(i-1)/(l-1),(j-1)/(h-1));
      if (normals) fprintf(fid, 'vn %f %f %f\n', nx(i,j),ny(i,j),nz(i,j)); end
      nn=nn+1;
    end
  end
  fprintf(fid,'g mesh\n');

  for i=1:(l-1)
    for j=1:(h-1)
      if (normals)
	fprintf(fid,'f %d/%d/%d %d/%d/%d %d/%d/%d %d/%d/%d\n',n(i,j),n(i,j),n(i,j),n(i+1,j),n(i+1,j),n(i+1,j),n(i+1,j+1),n(i+1,j+1),n(i+1,j+1),n(i,j+1),n(i,j+1),n(i,j+1));
      else
  	fprintf(fid,'f %d/%d %d/%d %d/%d %d/%d\n',n(i,j),n(i,j),n(i+1,j),n(i+1,j),n(i+1,j+1),n(i+1,j+1),n(i,j+1),n(i,j+1));
      end
    end
  end
  fprintf(fid,'g\n\n');
  fclose(fid);

With the function above you can export a functioning 3D object in Octave, and one that will accept colours since we now have a mesh based 3D object. The produced image is so far in grey-scale, so we still need to apply a texture in order to bring back the colours.

↑ Before and after mesh colourisation (texture mapping)

 

Texture mapping in MeshLab

The process of texture mapping is somewhat similar to what we do in Octave/MATLAB when we apply a colormap. In 3D modelling the process is defined by applying an JPEG or PNG image (which essentially is a colormap) as a texture to the 3D object. If you use MeshLab, you can apply a texture by:

  1. Loading/importing your 3D object into MeshLab.
  2. Go to: Filters (in top menu) → Texture → Set texture
  3. Supplying the correct file name to the texture map. Note: if you get an error when doing this your 3D model probably still lacks the needed mesh structure.
  4. Go to the top toolbar and activate ‘Texture’, via the texture button to the left of the ‘Light on/off’ button. Notice that the true colours will show when you turn off the light.

You can now export your mesh/3D model as a new mesh format, and thereby preserving your texture mapping. Note that even though colour may be included when you export your object, it doesn’t imply that all textures are included in your object file. You can include texture mapping in Wavefront object files, but the information of the texture is retained in your texture file(s).

3D representations of 2D images in MATLAB

The idea of an image is quite different when you try to perceive it as a matrix, especially if you twist some things around. So what does an image look like if you perceive it as a matrix with the RGB value as a z-axis while keeping the x- and y-axes ? How can you produce a 3D representation of a 2D image? Something similar to the below?

The idea of producing a 3D image from a 2D image steamed from a need of analysing and validating the roughness in microscope 2D images. My reasoning was that it had to be easier to analyse an image in three dimensions, compared to two dimensions. In hindsight I should have realised that for my purpose I was utterly wrong…

After some programming in MATLAB/Octave I had achieved my script that converts 2D images to 3D, and I had my output. However, the data you get when inserting a grey-scale image and converting it to 3D is quite overwhelming. You actually don’t get any more information out of it due to the large variation in a small grey-coloured interval. And a simple filtering didn’t really clear up the image either. The result was clearly not helpful.

2D and 3D representations of microscope images

So, what do you do with some code that is more or less useless? Well, you have some fun with it. It also became apparent that whilst visualisation in grey-scale wasn’t optimal, visualisations in colour was quite interesting. Here you can see some images, with and without filtering

Without filtering

With filtering

The big downside of the script, as well as with MATLAB/Octave in general, is the exportation to a suitable format. Is it possible to export the result to a Wavefront object file? Yes, primarily due to the important checking that the axes are of similar dimensions/values. The representation is however in black and white. Since this hasn’t any real purpose, other than some fun scripting, I will leave it here (for now).

% This script transforms a 2D RGB image to a 3D image, where the Z axis is given by the RGB number.
% For the beginner:
% 1. Save the script as 'script_2Dimage_to_3D.m'
% 2. Run the script via run 'script_2Dimage_to_3D.m'
% 3. Insert path to image you want to analyse.
% Alex Björkholm 2016
pkg load image  % Note that this needs to be commented out in MATLAB

% Load image from input
prompt = 'Enter path to image to be processed: ';
str = input(prompt,'s');
rgb_image = imread(str);
fprintf('Analysing image "%s"\n',str);

% Filtering yes/no
prompt = 'Apply filtering (yes/no): ';
yes_no = input(prompt,'s');

% Filtering
if strcmp(yes_no, 'yes')
   printf('\nApplying filter\n');
   f = ones(5) * 1/25;
   rgb_filter = imfilter(rgb_image, f, "symmetric";);
   rgbim = rgb_filter;

   figure;
   subplot ( 1, 2, 1), imshow ( rgb_image ), title('unfiltered');
   subplot ( 1, 2, 2), imshow ( rgb_filter ), title('filtered');

else
   fprintf('\nNo filter applied\n');
   rgbim = rgb_image;

end

% Converting image format and getting colormap
[C,map] = rgb2ind(rgbim);

% Setting X,Y dimensions, where X and Y has the same dimensions.
[X,Y] = meshgrid(-floor(size(rgbim, 2)/2):floor(size(rgbim, 2)/2)-1, -floor(size(rgbim, 1)/2):floor(size(rgbim, 1)/2)-1);

% Dimensions X and Y must be the same as C. This can be accomplished by adjusting the meshgrid. Automatically done below.
i = 0;
j = 0;
while ( size(X)(1) ~= size(C)(1) ) % While dimensions of x different for X and C, Note X == Y.
   i++;
   [X,Y] = meshgrid(-floor(size(rgbim, 2)/2):floor(size(rgbim, 2)/2) + 10 - j, -floor(size(rgbim, 1)/2):floor(size(rgbim, 1)/2) + 10 - i);
end

while ( size(X)(2) ~= size(C)(2) ) % While dimensions of x different for X and C, Note X == Y.
   j++;
   [X,Y] = meshgrid(-floor(size(rgbim, 2)/2):floor(size(rgbim, 2)/2) + 10 - j, -floor(size(rgbim, 1)/2):floor(size(rgbim, 1)/2) + 10 - i);
end

% If error persists, print array dimensions of X, Y, and C
XC = isequal(size(X),size(C));
YC = isequal(size(Y),size(C));
if XC ~= 1 || YC ~= 1
   fprintf('size(X):   %i   %i\n', size(X)(1), size(X)(2));
   fprintf('size(Y):   %i   %i\n', size(Y)(1), size(Y)(2));
   fprintf('size(C):   %i   %i\n', size(C)(1), size(C)(2));
end

% Setting a reasonable restriction on the Z (C) axis, which in turn enables the exportation of a reasonable 3D object
k = max(C(:)) / max(X(:));
[Z] = double(C/k);

% Setting 3D figure output
figure;
surf(X, Y, Z, double(C)), shading flat;
colormap(map); % This command returns the original image colour

3D enabled PDF?

column

2D image of a FreeCAD object

The big question when it comes to scientific 3D data; how do you most efficiently present it? A 2D image isn’t really that optimal for it. This is my quest for the optimal solution.

I would prefer a universal solution for the presentation of 3D data and figures, which means that it works across different platforms and operating systems. My quest started with 3D enabled PDFs, since PDF stands for Portable Document Format and I thought this must be the optimal solution. I quickly realised it wasn’t really that optimal, due to two primal factors:

  1. The 3D enhanced PDFs could only be interpreted by Desktop versions of Adobe Reader. This leaves out all UNIX, Linux, Mac, and Android systems. Even the iPad version of Adobe Reader is currently lacking 3D functions.
  2. 3D enhanced PDFs are highly dependable on the u3d image format and it really isn’t that universal in Linux. I tried finding a converter/program for conversion of any 3D image format to u3d. The only programs that should have this capability, native to Linux, are MeshLab and Blender, which generally works excellent but they still lack a functioning plugin for u3d conversion.

The solution for a missing u3d-exporter in Linux? Crazy as it seems, the easiest solution here is to download a Windows-copy of MeshLab and Wine the whole thing. A pretty sad solution, however, which all to often needs to be applied in Linux to get some things to work properly.

So when we now have a u3d-image; how is a 3D enabled PDF created using LaTeX? Credits goes here to Nicola Rainiero and his excellent tutorial, with all the necessary files included for trying the whole thing out (check out the zip-file at the bottom of the page).

I will here briefly go through the different steps included and provide the relevant data.

The LaTeX code

\documentclass[a4paper]{article}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{setspace}
\usepackage{caption}
\usepackage{fancyref}
\usepackage[3D]{movie15} %the style package that enables 3D rendering
\usepackage[english]{babel}
\usepackage[pdftex]{hyperref}&amp;amp;nbsp; %load after mathtools, caption, cite for links
\usepackage{cite}&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp;&amp;amp;nbsp; %load after hyperref for links

\input glyphtounicode
\pdfgentounicode=1
\setlength{\parindent}{0pt}
\onehalfspace
\pagenumbering{roman}

\numberwithin{equation}{section}
%\numberwithin{figure}{section}
%\numberwithin{table}{section}

\captionsetup{font=small, labelfont=bf, width=0.95\textwidth} %define after textwidth
%\setlength{\abovecaptionskip}{0pt}
%\setlength{\belowcaptionskip}{7.5pt}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}

\begin{document}
\includemovie[poster,
   toolbar,
   label=pt,
   text={\includegraphics[scale=1.0]{image2D.png}},
   3Droo=6.896200246789337,
   3Daac=60.000001669652114,
   3Dcoo=0.6134188175201416 0.6502023935317993 -0.8552163243293762,
   3Dc2c=-0.8354106545448303 0.3235208988189697 -0.44432342052459717,
   3Droll=-75.5946486014902,
   3Dlights=Hard,
   3Drender=SolidOutline]
   {\linewidth}{\linewidth}{image3D.u3d}\\
A 3D, schematic model. If the figure does not appear interactive, please enable this function and click on it or use a recent version of \href{http://get.adobe.com/reader}{Adobe Reader}.
\end{document}

The files that you will need are the following:

  • The 3D image “image3D.u3d”
  • The 2D image “image2D.png”
  • The “movie15” package for LaTeX. Available from CTAN.

When you have your LaTeX file and the relevant files in the same directory; the only thing left is for you to run your LaTeX-command. I prefer using pdflatex since it skips the unnecessary steps of creating DVI and PostScript files with the original latex command.

pdflatex name_of_your_pdf_file

The result? A PDF file with a 3D enhanced image, if you have access to Adobe Reader. latex-3d-pdf. Is this the optimal solution? Not really, but it is a very neat solution if you have access to Adobe Reader. My quest however, for the optimal solution, continues…

Plot script for nuclear cross section data

When I some years ago began working with nuclear cross section calculations, via the reaction codes EMPIRE and TALYS, my initial thought was; “No, I can’t start plotting all these cross section plots individually. There has to be a better way”. I think this was when I really started to code in Bash.

There was an easy solution. It took some scripting, but the solution was there. It didn’t take long either before I have my first working version of the cross section plotting program. It was a simple Bash script tied in with the plotting capabilities of gnuplot. Working in Linux gives some perks, especially through Bash with its easy syntax, wide availability to different sets of tools (for example find, grep, and sed), and the ability to run program scripts inherently. The requirements for my program grew, and I continued developing. Revision after revision the kinks and problems have ebbed away and a solid script has emerged over the years.

Why this solution?

As usual in coding the existing program really didn’t cut it, and I needed functions that wasn’t supplied by the existing solutions. One such solution is EMPIRE’s native ZVD plotting program. The advantages with this script are the following:

  • The ability to plot both EMPIRE, TALYS, EXFOR, TENDL, and your own experimental and theoretical cross section data side by side.
  • Plotting from a database based on a calculation search strings, which enables multiple plotting in one go. Note that plotting a massive numer of data sets might result in the plot legend (key in gnuplot) overflowing your gnuplot window.
  • The option to plot product cross sections for metastable and/or ground states.
  • The ability to plot data on natural targets from standard EMPIRE and TALYS calculations. I.e. the script will automatically correct the isotope cross section based on the  natural occurrence of said element.
  • The option to plot yields as well as all cross sections from a certain reaction.
  • An inbuilt ability to plot a certain energy range.
  • Excellent options for data output, and the possibility to continue analysing the data via gnuplot or other scripts.

Working notes on the code

There are still some issues with the script, mainly concerning the speed of the program. This can partly be traced down to a incomparability between the EXFOR data format (the EXFOR database is the EMPIRE copy of the latest X4toc4 master version) and gnuplot’s input formats.  The problem stems from the different notations of scientific number formats, and my workaround isn’t extremly fast. Another speed-related problem of the program will occur when/if you have amassed large libraries of nuclear data. It boils down to the simple fact that it takes time to search through massive data libraries. The solution to this is optimising the current processes as well as developing parallell processes. Sadly this isn’t really the strenght of Bash. The development in these areas will however continue.

There are also  some database related issues and functionalities in development. One is the implementation of an automatic TENDL database. Right now you will need to locate the necessary data files from the TENDL database and save them to your own file system in order to search and plot TENDL cross section data.

License / tips?

Feel free to use the program and contact me if you have any questions regarding it. If you find it helpful, please let me know and I will contine to update the code also here.

Requirements

What are needed in order to run the program?

  • Bash. If you want to run it in some other shell it should be possible after some minor modifications to the code.
  • gnuplot (preferably 5.0+). I highly recommend using the script with gnuplot v5.0+ in order to have an responsive based legend, whereby you can display and hide plots by clicking the corresponding legend key. The script also works on older versions of gnuplot, but you might need to change the output terminal from wxt to qt in code line 1,221.
  • data-plot-bash. The placement of the script on your hard drive doesn’t matter as long as you have set up your bash aliases correctly in your bashrc file.
  • data-plot-ref. This should be saved in the same directory as the script.

The code

#!/bin/bash
#
# This script runs Gnuplot for several EMPIRE and/or TALYS data files
# Gnuplot plots data for total cross section comparison with EXFOR, TENDL
# and experimental data files.
#
# Script name: data-plot
# Author: Alex Björkholm 2013-2016
# Date: 2016.08.16

################################ WORKING NOTES ################################

# NOTE: Always check with the EXFOR database for missing data. https://www-nds.iaea.org/exfor/exfor.htm
# NOTE: TENDL database not implemented

################################# DECLARATIONS ################################

echo -e "\n
     ---------------------------\n
          CROSS SECTION PLOT    \n
     ---------------------------\n"
IFS=$'\n'

declare -a data_plot
declare -a input_data_bases
declare -a data_bases
declare -a input_search

declare -a dir_empire
declare -a dir_talys
declare -a dir_tendl
declare -a dir_exper
declare -a dir_exper
declare -a dir_gnuplot

declare -a file_paths_empire
declare -a file_paths_talys
declare -a file_paths_tendl
declare -a file_paths_exfor
declare -a file_paths_exper

declare -a input_reaction
declare -a input_products
declare -a natural_isotopes

declare -a targets
declare -a targetsAA
declare -a targetsZZ
declare -a targetsNN
declare -a targetsXX
declare -a targetsPP
declare -a targetsMM
declare -a temp_targets
declare -a temp_targets_MM
declare -a temp_targets_PP

declare -a beams
declare -a beamsAA
declare -a beamsZZ
declare -a beamsNN
declare -a beamsXX
declare -a temp_beams

declare -a products
declare -a productsAA
declare -a productsZZ
declare -a productsNN
declare -a productsXX
declare -a productsMM
declare -a input_products
declare -a temp_products_MM
declare -a temp_products

declare -a reactions
declare -a temp_reactions
declare -a emissions

declare -a delete_array

cd $HOME

################################## VARIABLES ##################################

# Paths $RESEARCHDIR, $RESEARCHPROGRAMDIR, $EMPIREDIR, $TALYSDIR must be defined in bash

# Directory for this program and its reference data
programdir="$RESEARCHPROGRAMDIR/Cross section - Data plot scripts - Bash Gnuplot"
reference="$programdir/data-plot-ref.dat"

# Directories where EMPIRE or TALYS calculations are stored. Mandatory format for subdirectories is directories starting with "target,beam" with subsequent tags separated with "-". For example "75Br,h", "208Pb,p-run", "209Bi,a-omp-rot", and "202Hg,10B-omp".
dir_empire+=($(find $RESEARCHDIR -maxdepth 2 -type d -iwholename $RESEARCHDIR\*Calculations\ Empire\*))
dir_talys+=($(find $RESEARCHDIR -maxdepth 2 -type d -iwholename $RESEARCHDIR\*Calculations\ Talys\*))

# Directories where EXFOR and TENDL data are stored.
# Exfor data (download from http://www-nds.iaea.org/x4toc4-master/) should be sorted with the parseC4.py script in Empire to achieve the correct structure.
# Tendl data (download from ftp://ftp.nrg.eu/pub/www/talys/tendl2014/tendl2014.html) should be stored in the structure provided from the Tendle database.
dir_tendl+=($(find $RESEARCHDIR -maxdepth 3 -type d -iwholename $RESEARCHDIR\*Experimental\ and\ theoretical\ data\*tendl\*))
dir_exfor="$EMPIREDIR/EXFOR"
#dir_tendl="$TALYSDIR/TENDL"

# Directories where additional data in tab separated format is stored. Mandatory format for subdirectories is directories starting with "target,beam" with subsequent tags separated with "-". For example "75Br,h", "209Bi,a-omp-rot", and "202Hg,10B-omp".
dir_exper+=($(find $RESEARCHDIR -maxdepth 3 -type d -iwholename $RESEARCHDIR\*Experimental\ and\ theoretical\ data\*exper\* -not -iwholename \*chisq\*))
dir_gnuplot+=($(find $RESEARCHDIR -maxdepth 3 -type d -iwholename $RESEARCHDIR\*Experimental\ and\ theoretical\ data\*gnuplot\*))

# Directories where output from plots are saved.
mkdir -p $HOME/gnuplot-plots
gnuplot_plot="$HOME/gnuplot-plots"

# Directory for temp files used by GNUPLOT.
mkdir -p $HOME/gnuplot-temp
gnuplot_temp="$HOME/gnuplot-temp"

################################## FUNCTIONS ##################################

manual () {
   echo -e "This script runs Gnuplot for several EMPIRE and/or TALYS data files.
   Program data-plot input string:

   1.   Cross section data to plot
        product options: xs, direct, preequilibrium, compound, all-xs
        production options: activity, \#isotopes, yield, production
        emission options: xs, yield, all-xs
   2.   System(s) to plot
        options: exfor, exper, gnuplot, tendl, talys, empire
        multiple options can be chosen by writing option1_option2_option3
   3.   Reaction input
        options: '75As(h,x)76Br-75Br', 'natCr(p,x)52Mn', '52Cr(p,x)52mMn', '52Cr(p,x)52gMn' etc.
        multiple options can be chosen by writing option1_option2_option3
        product note: set outgoing particles to 'x' (calculated automatically)
   4.   Search string for EMPIRE and TALYS calculations
        multiple options can be chosen by writing option1_option2_option3
        wildcard notation * can be used to get the optimal selection
   5.   Graphical terminal type
        options: edit/no-edit/eps/epslatex/svg/table
   6.   Optional x-axis range in the form '5-33'
        if no range is defined the x-axis will autoscale to data

   Example: data-plot xs talys_empire 'natCr(p,x)52Mn' 'optimal' no-edit 5-40
            data-plot xs exfor_exper '75As(h,x)76Br-75Br' 'calculation' edit 5-50
            data-plot xs exfor_exper '75As(h,x)all' 'omp*HFB' edit 5-50\n"
exit 2 ; }

natural_conversion () {
local input=$1
unset AAA; unset NAA; unset NAP

if [[ $input != *[0-9]* ]] && [[ $input == nat[A-Z]* ]]  ; then
NZZ=$(echo $input | sed 's/nat//g')
zZZ=$(awk -v src=$NZZ '$2 == src { print $1 }' $reference)
ZZZ=($(echo ${zZZ} | sed 's/^0*//'))

while IFS=' ' read -r -a line; do
# Extracting isotope atomic numbers for natural target
AAA+=($(echo ${line[1]} | sed 's/^0*//'))
# Extracting isotope percentage for natural target
PP+=($( bc -l <<< "scale=6;${line[2]} / 100" | sed 's/^\./0./'))
done < $TALYSDIR/structure/abundance/z${zZZ}

fi

if [[ $AAA != *[0-9]* ]] || [[ $ZZZ != *[0-9]* ]]; then
echo -e "ERROR: natural_conversion function for input '$input'\n\n"; exit 2;
fi ; }


nuclide_conversion () {
local input=$1

# If format of the beam is n, p, h, a etc., resolving ZZZ and AAA
if [[ $input != *[0-9]* ]]; then
AAA=$(awk -v src=$input '$1 == src { print $5 }' $reference)
NZZ=${input}
zZZ=$(awk -v src=$input '$1 == src { print $3 }' $reference)
if [[ $input != n ]]; then ZZZ=($(echo ${zZZ} | sed 's/^0*//')); fi
if [[ $input == n ]]; then ZZZ=($(echo ${zZZ})); fi
fi

# If format of target is ZZZAAA (083209), resolving ZZZ and AAA
if [[ $input != *[a-z]* ]]; then
AAA=$(echo ${input:3:7} | sed 's/^0*//')
ZZZ=$(echo ${input:0:3} | sed 's/^0*//')
zZZ=$(echo ${input:0:3})
NZZ=$(awk -v src=$zZZ '$1 == src { print $2 }' $reference)
fi

# If format of the nuclide  is AAAXX (209Bi), resolving ZZZ and AAA
if [[ $input == *[0-9][A-Z]* ]]; then
AAA=$(echo $input | sed 's/[^0-9]*//g')
NZZ=$(echo $input | sed 's/^[0-9]*//g')
zZZ=$(awk -v src=$(echo $NZZ | sed 's/^m//g') '$2 == src { print $1 }' $reference)
ZZZ=($(echo ${zZZ} | sed 's/^0*//'))
fi

if [[ $AAA != *[0-9]* ]] || [[ $ZZZ != *[0-9]* ]]; then
echo -e "ERROR: nuclide_conversion function for input '$input'\n\n"; exit 2
fi
NNN=$(($AAA - $ZZZ)); }


p-n-emission () {
local particle_emission=$1
emission_number=0 ; NOUT=0 ; POUT=0

if [ $particle_emission == gamma ]; then particle_emission=g; else
particle_emission=$(echo $particle_emission | grep -o '[0-9]*[a-z]')

for each_emission in ${particle_emission[@]}; do
emission_particle=$(echo $each_emission | grep -o '[a-z]')
emission_number=$(echo $each_emission | grep -o '[0-9]*')

if [ -z $emission_number ]; then emission_number=1; fi

n=$(awk -v src=$emission_particle '$1 == src { print $4 }' $reference)
p=$(awk -v src=$emission_particle '$1 == src { print $3 }' $reference)
NOUT=$(($NOUT + $emission_number * $n))
POUT=$(($POUT + $emission_number * $p))

done

fi; }


particle-emission () {
NOUT=$(($TNN + $BNN - $PNN))
POUT=$(($TZZ + $BZZ - $PZZ))

particle_emission=${NOUT}n${POUT}p
if [ $particle_emission == 1n1p ]; then particle_emission=d; fi
if [ $particle_emission == 2n1p ]; then particle_emission=t; fi
if [ $particle_emission == 1n0p ]; then particle_emission=n; fi
if [ $particle_emission == 0n0p ]; then particle_emission=gamma; fi
if [[ $particle_emission == *n0p ]]; then particle_emission=$(echo $particle_emission | sed 's/0p//g'); fi
if [[ $particle_emission == *-* ]]; then particle_emission=none; fi; }


MT_number () {
local MT=()
declare -a MT

# See endf-6 manual p. 306-316. Here simplified since we're only interested in the produced nuclei
if [ $POUT == 0 -a $NOUT == 0 ];then MT+=(2); fi
if [ $POUT == 0 -a $NOUT == 1 ];then MT+=(4); fi
if [ $POUT == 1 -a $NOUT == 3 ];then MT+=(11); fi
if [ $POUT == 0 -a $NOUT == 2 ];then MT+=(16); fi
if [ $POUT == 0 -a $NOUT == 3 ];then MT+=(17); fi
if [ $POUT == 2 -a $NOUT == 3 ];then MT+=(22); fi
if [ $POUT == 6 -a $NOUT == 7 ];then MT+=(23); fi
if [ $POUT == 2 -a $NOUT == 4 ];then MT+=(24); fi
if [ $POUT == 2 -a $NOUT == 5 ];then MT+=(25); fi
if [ $POUT == 1 -a $NOUT == 1 ];then MT+=(28); fi
if [ $POUT == 4 -a $NOUT == 5 ];then MT+=(29); fi
if [ $POUT == 4 -a $NOUT == 6 ];then MT+=(30); fi
if [ $POUT == 1 -a $NOUT == 2 ];then MT+=(32); fi
if [ $POUT == 1 -a $NOUT == 3 ];then MT+=(33); fi
if [ $POUT == 2 -a $NOUT == 2 ];then MT+=(34); fi
if [ $POUT == 3 -a $NOUT == 4 ];then MT+=(35); fi
if [ $POUT == 5 -a $NOUT == 7 ];then MT+=(36); fi
if [ $POUT == 0 -a $NOUT == 4 ];then MT+=(37); fi
if [ $POUT == 1 -a $NOUT == 2 ];then MT+=(41); fi
if [ $POUT == 1 -a $NOUT == 3 ];then MT+=(42); fi
if [ $POUT == 2 -a $NOUT == 1 ];then MT+=(44); fi
if [ $POUT == 3 -a $NOUT == 3 ];then MT+=(45); fi
if [ $POUT == 0 -a $NOUT == 0 ];then MT+=(103); fi
if [ $POUT == 1 -a $NOUT == 0 ];then MT+=(103); fi
if [ $POUT == 1 -a $NOUT == 1 ];then MT+=(104); fi
if [ $POUT == 1 -a $NOUT == 2 ];then MT+=(105); fi
if [ $POUT == 2 -a $NOUT == 1 ];then MT+=(106); fi
if [ $POUT == 2 -a $NOUT == 2 ];then MT+=(107); fi
if [ $POUT == 4 -a $NOUT == 4 ];then MT+=(108); fi
if [ $POUT == 6 -a $NOUT == 6 ];then MT+=(109); fi
if [ $POUT == 2 -a $NOUT == 0 ];then MT+=(111); fi
if [ $POUT == 3 -a $NOUT == 2 ];then MT+=(112); fi
if [ $POUT == 5 -a $NOUT == 6 ];then MT+=(113); fi
if [ $POUT == 5 -a $NOUT == 5 ];then MT+=(114); fi
if [ $POUT == 2 -a $NOUT == 1 ];then MT+=(115); fi
if [ $POUT == 2 -a $NOUT == 2 ];then MT+=(116); fi
if [ $POUT == 3 -a $NOUT == 3 ];then MT+=(117); fi
if [[ -z ${MT} ]] && [[ $POUT -ge 0 ]] && [[ $NOUT -ge 0 ]];then MT+=(5); fi
MTnum=(${MT[@]}); }


sort_function_underscore () {
if [ -z $1 ]; then echo -e "ERROR: missing input in 'sort_function_underscore'"; else
for each in $(echo $1 | sed 's/_/ /g' | tr ' ' '\n' | sort -u); do echo $each; done; fi ; }

sort_function_bar () {
if [ -z $1 ]; then echo -e "ERROR: missing input in 'sort_function_bar'"; else
for each in $(echo $1 | sed 's/-/ /g' | tr ' ' '\n' | sort -u); do echo $each; done; fi ; }





############################# DEFINITION OF INPUT #############################

# CONTROLLING AND RESOLVING PARAMETERS
# Setting output and checking input for product/emission

data_plot=($(sort_function_underscore $1))
input_data_bases=($(sort_function_underscore $2))
input_reaction=($(sort_function_underscore $3))
input_search=($(sort_function_underscore $4))


# Setting output and checking input for plot axis
for each in ${data_plot[@]}; do
if [[ ${input_reaction} == *x\)* ]]; then
case $each in
man ) manual; ;;
xs | direct | preequilibrium | compound | all-xs ) echo -e "Plotting product cross section $each\n"; ;;
activity | \#isotopes | yield | production) echo -e "Plotting medical nuclide production $each\n"; ;;
* ) echo -e "Input cross section '$each' incorrect for product plot\n"; exit 2; ;;
esac; fi
if [[ ${input_reaction} != *x\)* ]]; then
case $each in
man ) manual; ;;
xs | yield | all-xs ) echo -e "Plotting emission cross section $each\n"; ;;
* ) echo -e "Input cross section '$each' incorrect for emission plot\n"; exit 2; ;;
esac; fi
done


# Setting output and checking input for systems
for each_system in ${input_data_bases[@]}; do
case $each_system in
empire | exfor | exper | talys | tendl | gnuplot ) data_bases+=(${each_system}); echo -e "Plotting ${each_system^^} data"; ;;
* ) echo -e "Input '${each_system^^}' for data system to plot incorrect\n"; exit 2; ;;
esac
done


case $5 in
edit | no-edit | eps | epslatex | svg | table ) gnuplot_output=$5; ;;
man ) manual; ;;
* ) echo -e "Input '$5' incorrect. Choose edit/no-edit/eps/epslatex/svg or 'man' for manual.\n"; exit 2; ;;
esac


case $6 in
*-* ) gnuplot_x_range=$6; echo; ;;
*) echo -e "Input '$6' for x-axis range missing. Setting auto scale\n"; ;;
esac






######################## RESOLVING REACTION PARAMETERS ########################

for each_input_reaction in ${input_reaction[@]}; do

# Setting 'reactions' array for input target and beam
input_target=($(echo $each_input_reaction | awk -F"(" '{print $1}'))
input_beam=($(echo $each_input_reaction | awk -F"(" '{print $2}' | awk -F"," '{print $1}'))

# Setting natural target isotopes
if [[ ${input_target} == *nat* ]]; then
natural_conversion ${input_target}
for (( N = 0; N < ${#PP[*]}; N++ )); do
temp_targets_PP+=($(echo ${PP[N]}))
temp_targets_MM+=($(echo nat))
temp_targets+=($(echo ${AAA[N]}${NZZ}))
done

else
# Non natural target isotopes and
if [[ ${input_target} == *[0-9][mg][A-Z][a-z] ]]; then
if [[ ${input_target} == *[0-9][g][A-Z]* ]]; then temp_targets_MM+=($(echo g)); fi
if [[ ${input_target} == *[0-9][m][A-Z]* ]]; then temp_targets_MM+=($(echo m)); fi
input_target=$(echo ${input_target} | sed 's/^\(.\{2\}\)\(.\{1\}\)\(.*\)/\1\3/')
else
temp_targets_MM+=($(echo n))
fi
temp_targets+=(${input_target})
temp_targets_PP+=($(echo 1.000000))

fi

for (( N = 0; N < ${#temp_targets[*]}; N++ )); do
temp_beams+=(${input_beam})
done


# Setting 'reactions' array for input emission and products
input_emissions=($(sort_function_bar $(echo $each_input_reaction | awk -F"," '{print $2}' | awk -F")" '{print $1}')))
input_products=($(sort_function_bar $(echo $each_input_reaction | awk -F"," '{print $2}' | awk -F")" '{print $2}')))

# Product dependent sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ ${input_products[@]} != x ]]; then

if [[ ${input_products[@]} == *all* ]]; then
# For all products input
for (( N = 0; N < ${#temp_targets[*]}; N++ )); do
for (( n = 0; n < 10; n++ )) ; do
for (( p = 0; p < 5; p++ )) ; do
nuclide_conversion ${temp_beams[N]}
BAA=${AAA}; BZZ=${ZZZ}
nuclide_conversion ${temp_targets[N]}
TAA=${AAA}; TZZ=${ZZZ}

AAA=$(($TAA + $BAA - $p - $n)); if [ $AAA -le 1 ]; then AAA=1; fi
ZZZ=$(($TZZ + $BZZ - $p)); if [ $ZZZ -le 1 ]; then ZZZ=1; fi
if [ $ZZZ -lt 100 ]; then zZZ="0$ZZZ"; else zZZ="$ZZZ"; fi
if [ $ZZZ -lt 10 ]; then zZZ="00$ZZZ"; fi
NZZ=$(awk -v src="$zZZ" '$1 ~ src { print $2 }' $reference)

temp_products+=(${AAA}${NZZ})
temp_products_MM+=($(echo NaN))

done
done
done

else
# For specific products input
for each_product in ${input_products[@]}; do
if [[ ${each_product} == *[0-9][mg][A-Z][a-z] ]]; then
if [[ ${each_product} == *[0-9][g][A-Z]* ]]; then temp_products_MM+=($(echo g)); fi
if [[ ${each_product} == *[0-9][m][A-Z]* ]]; then temp_products_MM+=($(echo m)); fi
each_product=$(echo ${each_product} | sed 's/^\(.\{2\}\)\(.\{1\}\)\(.*\)/\1\3/')
else
temp_products_MM+=($(echo NaN))
fi
temp_products+=(${each_product})
done

fi

# Resolving reaction arrays
for (( N = 0; N < ${#temp_targets[*]}; N++ )); do
for (( M = 0; M < ${#temp_products[*]}; M++ )); do

# Resolving targets
nuclide_conversion ${temp_targets[N]}
targetsAA+=(${AAA}); targetsZZ+=(${ZZZ})
targetsNN+=(${NNN}); targetsXX+=(${NZZ})
targetsPP+=(${temp_targets_PP[N]})
targetsMM+=(${temp_targets_MM[N]})
targets+=(${temp_targets[N]})

# Resolving beams
nuclide_conversion ${temp_beams[N]}
beamsAA+=(${AAA}); beamsZZ+=(${ZZZ})
beamsNN+=(${NNN})
beams+=(${temp_beams[N]})

# Resolving products
nuclide_conversion ${temp_products[M]}
productsAA+=(${AAA}); productsZZ+=(${ZZZ})
productsNN+=(${NNN}); productsXX+=(${NZZ})
productsMM+=(${temp_products_MM[M]})
products+=(${temp_products[M]})

# Resolving reactions
reactions+=($(echo ${temp_targets[N]},${temp_beams[N]}))

done
done

else
# Emission dependent sequence
# For all emission input
if [[ $input_emissions != *all* ]] && [[ ${input_products} == x ]]; then
echo -e "This option is not yet implemented \n"; exit 2

# For specific emission input
else
echo -e "This option is not yet implemented \n"; exit 2
fi

fi

done


echo -e "Reactions to be plotted"
for (( N = 0; N < ${#reactions[*]}; N++ )); do
printf " ${targets[N]}(${beams[N]},x)${products[N]} \t\t  target ocurrence: %10s \t\t target state: ${targetsMM[N]} \t\t product state: ${productsMM[N]}\n" ${targetsPP[N]}
done






##################### EMPIRE DATA AND PLOT SEQUENCE ###########################
echo -e "\nSorting EMPIRE data"

# Excluding irrelevant EMPIRE top directories based in reaction product and target
for (( N = 0; N < ${#reactions[*]}; N++ )); do
for (( M = 0; M < ${#dir_empire[*]}; M++ )); do
if [[ ${dir_empire[M]} == *${productsXX[N]}* ]] || [[ ${dir_empire[M]} == *${targetsXX[N]}* ]]; then
delete_array+=(${dir_empire[M]}); fi
done
done

dir_empire=($(for each in ${delete_array[@]}; do echo $each; done | sort -u)); unset delete_array
#for (( N = 0; N < ${#dir_empire[*]}; N++ )); do echo ${dir_empire[N]}; done


# Excluding file duplicities of EMPIRE
if [[ ${input_emissions[@]} == *x* ]] && [[ ${data_bases[@]} == *empire* ]]; then
temp_reactions=($(for each in ${reactions[@]}; do echo $each; done | sort -u))


# Finding files EMPIRE sequence
for each_dir in ${dir_empire[@]}; do
for search in ${input_search[@]}; do
for (( N = 0; N < ${#temp_reactions[*]}; N++ )); do
# Finding files for EMPIRE (product specification inside EMPIRE files)
if [[ $data_plot == *xs* ]]; then
file_paths_empire+=($(find ${each_dir} -type f -iwholename ${each_dir}\*${temp_reactions[N]}\*${search}\*.xsc -not -iwholename \*outputdir\*))
fi

# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == yield ]] || [[ $data_plot == \#isotopes ]]; then
echo -e "EMPIRE: no support for '$data_plot' data"
fi

if [[ ${#file_paths_empire[*]} == $no_files ]] || [[ ${#file_paths_empire[*]} == 0 ]]; then
echo -e "EMPIRE: invalid reaction '${temp_reactions[N]}' and '${search}' in '$(basename ${each_dir})'"
else no_files=${#file_paths_empire[*]}; fi

done
done
done
fi

file_paths_empire=($(for file_path in ${file_paths_empire[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_empire[@]}; do echo "'$file_path'"; done



# Checking if need to run sequense and correct cross section
if [ ${#file_paths_empire[*]} != 0 ]; then
if [ $data_plot == xs -o $data_plot == all-xs ]; then
for (( N = 0; N < ${#file_paths_empire[*]}; N++ )) ; do

# For product specific reactions. All reactions; resolved in reaction sequence
# Sequence based on extraction of all reaction (emission) channels from empire data file
if [[ ${input_emissions[@]} == *x* ]] && [[ ${input_products[@]} != *all* ]]; then

column_no=`cat ${file_paths_empire[$N]} | head -2 | tail -1 | wc -w`
for (( C = 8; C < $column_no; C++ )) ; do
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
dir_name=$(basename $(dirname ${file_paths_empire[$N]}))

# Finding the correct products to plot
reaction=$(echo ${dir_name} | awk -F"-" '{print $1}')
if [[ ${reactions[M]} == ${reaction} ]]; then
particle_emission=$(cat ${file_paths_empire[$N]} | head -2 | tail -1 | tr -s ' ' | cut -d ' ' -f$((C+1)) | sed 's/)//g; s/(z,//g')

p-n-emission $particle_emission
BAA=${beamsAA[M]}; BZZ=${beamsZZ[M]}
TAA=${targetsAA[M]}; TZZ=${targetsZZ[M]}

AAA=$(($TAA + $BAA - $POUT - $NOUT))
ZZZ=$(($TZZ + $BZZ - $POUT))
if [ $ZZZ -le 100 ]; then zZZ="0$ZZZ"; else zZZ="$ZZZ"; fi
if [ $ZZZ -lt 10 ]; then zZZ="00$ZZZ"; fi
NZZ=$(awk -v src="$zZZ" '$1 ~ src { print $2 }' $reference)

# Excluding irrelevant products
if [[ ${products[M]} == ${AAA}${NZZ} ]]; then

# Setting target with correct reaction abundance
if [[ ${input_target[@]} == *nat* ]]; then
columns="(\$1):(\$${C}*${targetsPP[M]})"
else
columns="(\$1):(\$${C})"
fi

[ -n "$cmd_emp" ] && cmd_emp=${cmd_emp}, || cmd_emp=''
cmd_emp="${cmd_emp} '${file_paths_empire[$N]}' using ${columns} w lp pi 2 lc palette frac 0.$((3 * $N / 5 + 2 * $C)) title 'EMPIRE-$dir_name-${AAA}${NZZ}'"

fi

fi

done
done
fi

if [[ ${input_emissions[@]} != *x* ]]; then
for (( N = 2; N <= 5; N++ )) ; do
[ -n "$cmd_emp" ] && cmd_emp=${cmd_emp}, || cmd_emp=''
dir_name=$(basename $(dirname $(dirname ${file_paths_empire[$N]})))
xs_name=$(cat ${file_paths_empire[$N]} | head -2 | tail -1 | tr -s ' ' | cut -d ' ' -f$((N+1)))
cmd_emp="${cmd_emp} '${file_paths_empire[$N]}' using 1:$N w lp pi 2 lc palette frac 0.$((3 * $N / 5)) title 'EMPIRE-${xs_name,,}-$dir_name'"

done
fi


done
fi
fi

echo EMPIRE: sorting done





##################### TALYS DATA AND PLOT SEQUENCE ############################
echo -e "\nSorting TALYS data"

# Excluding irrelevant TALYS top directories based in reaction product and target
for (( N = 0; N < ${#reactions[*]}; N++ )); do
for (( M = 0; M < ${#dir_talys[*]}; M++ )); do
if [[ ${dir_talys[M]} == *${productsXX[N]}* ]] || [[ ${dir_talys[M]} == *${targetsXX[N]}* ]]; then
delete_array+=(${dir_talys[M]}); fi
done
done

dir_talys=($(for each in ${delete_array[@]}; do echo $each; done | sort -u)); unset delete_array
#for (( N = 0; N < ${#dir_talys[*]}; N++ )); do echo ${dir_talys[N]}; done


# Finding files TALYS sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ ${data_bases[@]} == *talys* ]]; then
for each_dir in ${dir_talys[@]}; do
for search in ${input_search[@]}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do

# Finding product specific, natural target, files
if [[ ${input_target[@]} == *nat* ]]; then
search_target=$(echo nat${targetsXX[N]}); search_beam=$(echo ${beams[N]})
else
search_target=$(echo ${targets[N]}); search_beam=$(echo ${beams[N]})
fi

# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == yield ]] || [[ $data_plot == \#isotopes ]]; then
file_paths_talys+=($(find ${each_dir} -type f -iwholename \*${search_target},${search_beam}\*${search}\*/Y\*${productsZZ[N]}\*${productsAA[N]}\*.\* -not -name \*prev\*))
fi

# Finding product specific files
if [[ $data_plot == xs ]]; then
file_paths_talys+=($(find ${each_dir} -type f -iwholename \*${search_target},${search_beam}\*${search}\*/rp\*${productsZZ[N]}\*${productsAA[N]}.\* -not -name \*prev\*))
fi

if [[ ${#file_paths_talys[*]} == $no_files ]] || [[ ${#file_paths_talys[*]} == 0 ]]; then
echo -e "TALYS: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}' and '${search}' in '$(basename ${each_dir})'"
else no_files=${#file_paths_talys[*]}; fi

done
done
done
fi

file_paths_talys=($(for file_path in ${file_paths_talys[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_talys[@]}; do echo "'$file_path'"; done


# Checking if need to run sequense
if [ ${#file_paths_talys[*]} != 0 ]; then

# Defining plotting function for TALYS
cmd_talys () {
# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_talys[$N]})))
base_name=$(basename ${file_paths_talys[$N]} | sed 's/.tot//g')
product=$(echo $base_name | sed 's/^[A-Z]*//g;s/^[a-z]*//g' | awk -F"." '{print $1}')
nuclide_conversion ${product}

# Setting target with correct title/plot for natural targets
target=$(echo $dir_name | awk -F"-" '{print $1}' | awk -F"," '{print $1}')
base_target=$(echo $base_name | sed 's/.L[0-9][0-9]//' | awk -F"." '{print $2}' | sed 's/^0*//')
if [[ ! -z $base_target ]]; then
target="${target/nat/nat(${base_target})}"      # Replaces nat with $base_target
dir_name="${dir_name/nat/nat(${base_target})}"  # Replaces nat with $base_target
fi

# Setting product with correct title/plot for ground/meta state
if [[ ${base_name} == *.L[0-9][0-9]* ]] && [[ ${productsMM[@]} != *[g-m]* ]]; then
continue
else
if [[ $base_name == *[0-9].L00* ]]; then state="g"; fi
if [[ $base_name == *[0-9].L[0-9][1-9] ]]; then state="m"; fi
fi

# Setting target with correct reaction abundance
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ "nat(${targetsAA[$M]})${targetsXX[$M]}" == ${target} ]]; then targetPP=${targetsPP[M]}; fi
if [[ "nat${targetsXX[$M]}" == ${target} ]]; then targetPP=($(echo 1.000000)); fi
if [[ ${targets[$M]} == ${target} ]]; then targetPP=($(echo 1.000000)); fi
done

if [ $data_plot != all-xs -a $data_plot != production ]; then
if [ $data_plot != xs ]; then data_plot=${data_plot^}; fi
column_no=$(awk -v src=${data_plot} '{ for(i=1;i<=NF;i++){if ($i ~ src) {print i } } }' ${file_paths_talys[$N]})
if [[ ${file_paths_talys[$N]} == *Y*.tot ]]; then column_no=$(( column_no / 2)); else column_no=$(( column_no - 1)); fi
fi

# Setting correct column to plot in files
if [[ ${file_paths_talys[$N]} == *Y*.tot ]]; then
column_name=$(cat ${file_paths_talys[$N]} | head -6 | tail -1 | tr -s ' ' | cut -d ' ' -f$(( column_no * 2 )))
column_names=$(cat ${file_paths_talys[$N]} | head -6 | tail -1 | tr -s ' ' | cut -d ' ' -f$(( column_no * 2 + 1)))
ylabel="$column_name $column_names"
else
column_name=$(cat ${file_paths_talys[$N]} | head -5 | tail -1 | tr -s ' ' | cut -d ' ' -f$(( column_no + 1)))
fi

# Setting final variables
columns="(\$1):(\$${column_no}*$targetPP)"
title=$(echo TALYS-${column_name,,}-$dir_name-${AAA}${state}${NZZ})
[ -n "$cmd_tal" ] && cmd_tal=${cmd_tal}, || cmd_tal=''
cmd_tal="${cmd_tal} '${file_paths_talys[$N]}' using ${columns} w lp pi 2 lc palette frac 0.$((3 + 3 * $N / 5)) title '${title}'"; }


for (( N = 0; N < ${#file_paths_talys[*]}; N++ )) ; do

# Plotting all cross sections/medical isotope production choises
if [ $data_plot == all-xs -o $data_plot == production ]; then
if [ $data_plot == production ]; then MAXC=5; fi
if [ $data_plot == all-xs ]; then MAXC=6; fi
for (( column_no = 2; column_no < $MAXC ; column_no++ )) ; do cmd_talys; done

# Plotting specific cross section/medical isotope production choise
else
cmd_talys
fi

done
fi

echo TALYS: sorting done





##################### TENDL DATA AND PLOT SEQUENCE ############################
echo -e "\nSorting TENDL data"

# Excluding irrelevant TALYS top directories based in reaction product and target
for (( N = 0; N < ${#reactions[*]}; N++ )); do
for (( M = 0; M < ${#dir_tendl[*]}; M++ )); do
if [[ ${dir_tendl[M]} == *${productsXX[N]}* ]] || [[ ${dir_tendl[M]} == *${targetsXX[N]}* ]]; then
delete_array+=(${dir_tendl[M]}); break; fi
done
done

dir_tendl=($(for each in ${delete_array[@]}; do echo $each; done | sort -u)); unset delete_array
#for (( N = 0; N < ${#dir_tendl[*]}; N++ )); do echo ${dir_tendl[N]}; done


# Finding files TENDL sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ ${data_bases[@]} == *tendl* ]]; then
for each_dir in ${dir_tendl[@]}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do

# Finding product specific files
if [[ $data_plot == xs ]]; then
file_paths_tendl+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/\*${products[N]}\*.tot))
fi

# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == nuclides ]]; then
file_paths_tendl+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/Y\*${products[N]}\*.tot))
fi

if [[ ${#file_paths_tendl[*]} == $no_files ]] || [[ ${#file_paths_tendl[*]} == 0 ]]; then
echo -e "TENDL: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_tendl[*]}; fi

done
done
fi

file_paths_tendl=($(for file_path in ${file_paths_tendl[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_tendl[@]}; do echo "'$file_path'"; done



# Checking if need to run sequense and correct cross section xs or all-xs
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *tendl* ]]; then
if [ ${#file_paths_tendl[*]} != 0 ]; then
for (( N = 0; N < ${#file_paths_tendl[*]}; N++ )) ; do

# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_tendl[$N]})))
base_name=$(basename ${file_paths_tendl[$N]} | sed 's/.tot//g')

# Setting target with correct reaction abundance
target=$(echo $base_name | awk -F"," '{print $2}' | awk -F"(" '{print $1}' | sed 's/ //g')
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ ${targets[$M]} == ${target} ]]; then
targetPP=${targetsPP[M]}
break
fi
done

# Setting product with correct title/plot for ground/meta state
if [[ ${base_name} == *[0-9[0-9][g-m][A-Z][a-z]* ]] && [[ ${productsMM[@]} != *[g-m]* ]]; then
continue
fi

[ -n "$cmd_ten" ] && cmd_ten=${cmd_ten}, || cmd_ten=''
columns="(\$1):(\$2*$targetPP) w lp pi 1"
cmd_ten="${cmd_ten} '${file_paths_tendl[$N]}' using ${columns} lc palette frac 0.$((6 + 3 * $N / 5)) title '$base_name'"

done
fi
fi

echo TENDL: sorting done





##################### EXFOR DATA AND PLOT SEQUENCE ############################
echo -e "\nSorting EXFOR data"

# Finding files EXFOR sequence
if [[ ${input_emissions[@]} == *x* ]] && [[ $data_plot == xs ]] && [[ ${data_bases[@]} == *exfor* ]]; then
for each_dir in ${dir_exfor}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do

# Finding product specific files
if [[ ${beams[N]} != *[0-9]* ]]; then
file_paths_exfor+=($(find $dir_exfor -type f -iwholename $dir_exfor/${beams[N]}\*/\*${targetsZZ[N]}\*${targetsXX[N]}\*${targetsAA[N]}.c4))
else
file_paths_exfor+=($(find $dir_exfor -type f -iwholename $dir_exfor/other/\*${targetsZZ[N]}\*${targetsXX[N]}\*${targetsAA[N]}.c4))
fi

if [[ ${#file_paths_exfor[*]} == $no_files ]] || [[ ${#file_paths_exfor[*]} == 0 ]]; then
echo -e "EXFOR: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_exfor[*]}; fi

done
done
fi

file_paths_exfor=($(for file_path in ${file_paths_exfor[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_exfor[@]}; do echo "'$file_path'"; done



# EXFOR SEQUENCE
# Checking if need to run sequense and correct cross section
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *exfor* ]]; then
if [[ ${#file_paths_exfor[*]} != 0 ]]; then
for (( N = 0; N < ${#file_paths_exfor[*]}; N++ )) ; do

# Copying file so that the official EXFOR file keeps intact
cp ${file_paths_exfor[$N]} $gnuplot_temp
file_paths_exfor[$N]=$gnuplot_temp/$(basename ${file_paths_exfor[$N]})

# Producing plots for all products
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
NOUT=$((${targetsNN[M]} + ${beamsNN[M]} - ${productsNN[M]}))
POUT=$((${targetsZZ[M]} + ${beamsZZ[M]} - ${productsZZ[M]}))

# Number of individual data sets in EXFOR file
# each_MT_sets is the first line for each data set
MT_set=(); MT_number
for each_MT in ${MTnum[@]}; do
MT_set+=($(grep -n " "$each_MT"[[:alpha:]]*[[:blank:]]*]" ${file_paths_exfor[$N]} | grep "C4" | awk -F":" '{print $1}'))
done

for each_MT_set in ${MT_set[@]};do
each_MT_linei=$(($each_MT_set + 4)) # First line of data to plot in each data set
# Last line of data to plot in each data set. Taken from listen number of data sets in EXFOR file
each_MT_linef=$(($(sed -n $(($each_MT_linei - 3))p ${file_paths_exfor[$N]}| awk '{print $2}') + $each_MT_linei))
# Lines omitted by gnuplot and corrected counting from zero
unused_lines=$(($(head -n ${each_MT_linei} ${file_paths_exfor[$N]} | grep "\#" | wc -l) + 1))

# Omitting irrelevant data sets
not_xs_data_set=$(sed -n "$each_MT_linei"p ${file_paths_exfor[$N]} | sed 's/^\(.\{59\}\)\(.\{25\}\)\(.*\)/\2/')

if [[ -z "${not_xs_data_set// }" ]]; then
for (( P = $each_MT_linei; P < ${each_MT_linef}; P++ )); do
dataline=$(sed -n ${P}p ${file_paths_exfor[$N]})

# Setting correct scientific format for data column 8
column08pm=$(echo ${dataline} | sed 's/^\(.\{56\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column08pm// } == + ]] || [[ ${column08pm// } == - ]]; then
sed "${P} s/^\(.\{56\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{58\}\)\(.\{1\}\)\(.*\)/\1  \3/" -i ${file_paths_exfor[$N]}
fi

# Setting empty value in column 7 to zero
column07=$(echo ${dataline} | sed 's/^\(.\{50\}\)\(.\{8\}\)\(.*\)/\2/')
if [[ -z "${column07// }" ]]; then
sed "${P} s/^\(.\{50\}\)\(.\{1\}\)\(.*\)/\10\3/" -i ${file_paths_exfor[$N]}; fi

# setting correct scientific format for data column 7
column07pm=$(echo ${dataline} | sed 's/^\(.\{47\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column07pm// } == + ]] || [[ ${column07pm// } == - ]]; then
sed "${P} s/^\(.\{47\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{49\}\)\(.\{1\}\)\(.*\)/\1\2\2\3/" -i ${file_paths_exfor[$N]}
fi

# Setting correct scientific format for data column 6
column06pm=$(echo ${dataline} | sed 's/^\(.\{38\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column06pm// } == + ]] || [[ ${column06pm// } == - ]]; then
sed "${P} s/^\(.\{38\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{40\}\)\(.\{1\}\)\(.*\)/\1\2\2\3/" -i ${file_paths_exfor[$N]}
fi

# Setting empty value in column 5 to zero
column05=$(echo ${dataline} | sed 's/^\(.\{32\}\)\(.\{8\}\)\(.*\)/\2/')
if [[ -z "${column05// }" ]]; then
sed "${P} s/^\(.\{32\}\)\(.\{1\}\)\(.*\)/\10\3/" -i ${file_paths_exfor[$N]}; fi

# setting correct scientific format for data column 5
column05pm=$(echo ${dataline} | sed 's/^\(.\{29\}\)\(.\{1\}\)\(.*\)/\2/')
if [[ ${column05pm// } == + ]] || [[ ${column05pm// } == - ]]; then
sed "${P} s/^\(.\{29\}\)\(.\{1\}\)\(.*\)/\1E\2\3/" -i ${file_paths_exfor[$N]}
else
sed "${P} s/^\(.\{31\}\)\(.\{1\}\)\(.*\)/\1\2\2\3/" -i ${file_paths_exfor[$N]}
fi

done

# Energy plotted to Cross section; with energy error lower and upper limit ; with cross section error lower and upper limit
# Setting target with correct reaction abundance
columns="(\$5*1E-6):(\$7*1E+3*${targetsPP[M]}):(\$6*1E-6):(\$8*1E+3*${targetsPP[M]}) w xyerrorbars"
every="every ::$(($each_MT_linei - $unused_lines))::$(($each_MT_linef - $unused_lines))"

AUTHOR=$(sed -n "$each_MT_linei"p ${file_paths_exfor[$N]} | sed 's/^\(.\{98\}\)\(.*\)/\2/' | awk  -F"(" '{print $1}')
AUTHOR=$(echo ${AUTHOR[@],,} | sed -r 's/\<./\U&/g' | sed 's/([0-9]*)//g;s/[[:space:]]//g;s/\./\. /g;s/,Et. Al./ et al./g' )
YEAR=$(sed -n "$each_MT_linei"p ${file_paths_exfor[$N]} | awk  -F")" '{print $1}' | awk  -F"(" '{print $2}')
if [[ $YEAR > 20 ]]; then YEAR=$(echo 19${YEAR}); else YEAR=$(echo 20${YEAR}); fi
TITLE=$(echo ${AUTHOR} \(${YEAR}\), ${targets[M]}\(${beams[M]},${NOUT}n${POUT}p\)${products[M]})
TITLE=$(echo $TITLE | sed 's/  / /g;s/0p//g;s/1p/p/g;s/0n//g;s/1n/n/g;s/protons/p/g;s/alphas/a/g;s/deutrons/d/g;s/helions/h/g;')

# Plotting
if [ -n "${columns}" ]; then
[ -n "$cmd_exf" ] && cmd_exf=${cmd_exf},  || cmd_exf=''
cmd_exf="${cmd_exf} '${file_paths_exfor[$N]}' ${every} using ${columns} lc palette frac 0.$((0 + 4 * $P / 2)) title '$TITLE'"
fi

fi
done

done
done

fi
fi

echo EXFOR: sorting done





######################## EXPER DATA AND PLOT SEQUENCE #########################
echo -e "\nSorting EXPER data"

# Finding files EXPER/EXFOR sequence
if [[ ${data_bases[@]} == *exper* ]]; then
for each_dir in ${dir_exper}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do

# Finding product specific files
if [[ $data_plot == xs ]]; then
if [[ ${targetsMM[N]} == nat ]]; then
file_paths_exper+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*nat${targetsXX[N]}\*${beams[N]}/\*${products[N]}\*))
else
file_paths_exper+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/\*${products[N]}\*))
fi
fi

# Finding data files for medical isotope production
if [[ $data_plot == production ]] || [[ $data_plot == activity ]] || [[ $data_plot == nuclides ]]; then
file_paths_exper+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/Y\*${products[N]}\*))
fi

if [[ ${#file_paths_exper[*]} == $no_files ]] || [[ ${#file_paths_exper[*]} == 0 ]]; then
echo -e "EXPER: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_exper[*]}; fi

done
done
fi

file_paths_exper=($(for file_path in ${file_paths_exper[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_exper[@]}; do echo "'$file_path'"; done



# EXPER/EXFOR SEQUENCE
# Checking if need to run sequense and correct cross section xs or all-xs
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *exper* ]]; then
if [ ${#file_paths_exper[*]} != 0 ]; then
for (( N = 0; N < ${#file_paths_exper[*]}; N++ )) ; do

# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_exper[$N]})))
base_name=$(basename ${file_paths_exper[$N]} | sed 's/.tot//g;s/.x4//g')

# Setting target with correct reaction abundance
target=$(echo $base_name | awk -F"," '{print $2}' | awk -F"(" '{print $1}' | sed 's/ //g')
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ ${targets[$M]} == ${target} ]]; then
targetPP=${targetsPP[M]}; break
else
targetPP=$((1))
fi
done

[ -n "$cmd_exp" ] && cmd_exp=${cmd_exp}, || cmd_exp=''
columns="(\$1*1E-6):(\$4*1E+3*$targetPP):((\$1*1E-6)-(\$2*1E-6)):((\$1*1E-6)+(\$3*1E-6)):(((\$4*1E+3)-(\$5*1E+3))*$targetPP):(((\$4*1E+3)+(\$6*1E+3))*$targetPP) w xyerrorbars"
cmd_exp="${cmd_exp} '${file_paths_exper[$N]}' using ${columns} lc palette frac 0.$((8 + 3 * $N / 5)) title '$base_name'"

done
fi
fi

echo EXPER: sorting done





####################### GNUPLOT DATA AND PLOT SEQUENCE ########################
echo -e "\nSorting GNUPLOT data"

# Finding files GNUPLOT sequence
if [[ ${data_bases[@]} == *gnuplot* ]]; then
for each_dir in ${dir_gnuplot[@]}; do
for (( N = 0; N < ${#reactions[*]}; N++ )); do

# Finding product specific files
if [[ $data_plot == xs ]] && [[ ${input_emissions[@]} == *x* ]]; then
file_paths_gnuplot+=($(find ${each_dir} -type f -iwholename ${each_dir}/\*${reactions[N]}/\*${products[N]}\*))
fi

if [[ ${#file_paths_gnuplot[*]} == $no_files ]] || [[ ${#file_paths_gnuplot[*]} == 0 ]]; then
echo -e "GNUPLOT: invalid $data_plot product '${products[N]}' for reaction '${reactions[N]}'"
else no_files=${#file_paths_gnuplot[*]}; fi

done
done
fi

file_paths_gnuplot=($(for file_path in ${file_paths_gnuplot[@]}; do echo $file_path; done | sort -u))
#for file_path in ${file_paths_gnuplot[@]}; do echo "'$file_path'"; done



# Checking if need to run sequense and correct cross section xs or all-xs
if [[ $data_plot == *xs* ]] && [[ ${data_bases[@]} == *gnuplot* ]]; then
if [ ${#file_paths_gnuplot[*]} != 0 ]; then
for (( N = 0; N < ${#file_paths_gnuplot[*]}; N++ )) ; do

# Resolving reaction parameters from file to get correct plot titles
dir_name=$(basename $(dirname $(dirname ${file_paths_gnuplot[$N]})))
base_name=$(basename ${file_paths_gnuplot[$N]} | sed 's/.tot//g')

# Setting target with correct reaction abundance
target=$(echo $base_name | awk -F"," '{print $2}' | awk -F"(" '{print $1}' | sed 's/ //g')
for (( M = 0; M < ${#reactions[*]}; M++ )) ; do
if [[ ${targets[$M]} == ${target} ]]; then
targetPP=${targetsPP[M]}; break
else
targetPP=$((1))
fi
done

[ -n "$cmd_gnu" ] && cmd_gnu=${cmd_gnu}, || cmd_gnu=''
columns="(\$1):(\$2*$targetPP) w lp pi 1"
cmd_gnu="${cmd_gnu} '${file_paths_gnuplot[$N]}' using ${columns} lc palette frac 0.$((9 + 3 * $N / 5)) title '$base_name'"

done
fi
fi

echo GNUPLOT: sorting done





######################### GNUPLOT PLOTLINE PARAMETERS #########################

# Checking plottable files, printing statistics
echo
echo -e "Found $(echo $cmd_emp | grep -o "title" | wc -l) EMPIRE data series to plot"
echo -e "Found $(echo $cmd_tal | grep -o "title" | wc -l) TALYS data series to plot"
echo -e "Found $(echo $cmd_ten | grep -o "title" | wc -l) TENDL data series to plot"
echo -e "Found $(echo $cmd_exf | grep -o "title" | wc -l) EXFOR data series to plot"
echo -e "Found $(echo $cmd_exp | grep -o "title" | wc -l) EXPER/EXFOR data series to plot"
echo -e "Found $(echo $cmd_gnu | grep -o "title" | wc -l) GNUPLOT data series to plot\n"


# MANAGING TITLES
# Managing plot titles
TIME=`date '+%Y.%m.%d-%H.%M'`
plot="Cross section data | ${data_plot[@]} | ${input_reaction[@]} | ${input_data_bases[@]^^} | $TIME"
plot_title="Cross section data: ${input_reaction[@]}| ${input_data_bases[@]^^} $data_plot cross section(s)"

# Managing plot labels and titels
if [[ ${input_reaction} == *x* ]]; then
data_plot=${data_plot,}
if [ $data_plot == production -o $data_plot == activity -o $data_plot == yield -o $data_plot == \#isotopes ]; then
xlabel="Time (h)"
plot=$(echo $plot | sed 's/Cross section data/Nuclide production data/g')
plot_title=$(echo $plot_title | sed 's/Cross section data/Nuclide production data/g')
plot_title=$(echo $plot_title | sed 's/cross section(s)//g')

else
xlabel="Energy (MeV)"; ylabel="Cross section (mb)"
fi

else
xlabel="Emission energy (MeV)"; ylabel="Cross section (mb)"
fi

# Cleaning plot titles
plot=$(echo $plot | sed 's/_/ /g;s/-/,/g;s/|/-/g' | sed -r 's/,+/,/g')
plot_title=$(echo $plot_title | sed 's/_/ /g;s/-/,/g;s/|/,/g' | sed -r 's/,+/,/g')




# MERGING CMD COMMANDS
# Merging commands for theoretical and experimental plots
if [ ! -z $cmd_emp ]; then cmd=($cmd_emp); fi
if [ ! -z $cmd_tal ]; then
if [ -z $cmd ]; then cmd=($cmd_tal)
else cmd=$(echo $cmd, $cmd_tal); fi; fi
if [ ! -z $cmd_ten ]; then
if [ -z $cmd ]; then cmd=($cmd_ten)
else cmd=$(echo $cmd, $cmd_ten); fi; fi
if [ ! -z $cmd_exf ]; then
if [ -z $cmd ]; then cmd=($cmd_exf)
else cmd=$(echo $cmd, $cmd_exf); fi; fi
if [ ! -z $cmd_exp ]; then
if [ -z $cmd ]; then cmd=($cmd_exp)
else cmd=$(echo $cmd, $cmd_exp); fi; fi
if [ ! -z $cmd_gnu ]; then
if [ -z $cmd ]; then cmd=($cmd_gnu)
else cmd=$(echo $cmd, $cmd_gnu); fi; fi




# Checking for errors in cmd
if [ -z $cmd ]; then
echo -e "Found no data to plot. Check reaction parameters. \n"; exit 2
else echo -e "Plotting..."
fi




# GENERATING GNUPLOT SCRIPT
# Setting Gnuplot script parameter for printing
if [ $gnuplot_output == eps -o $gnuplot_output == epslatex -o $gnuplot_output == svg -o $gnuplot_output == table ]; then
if [ $gnuplot_output == eps ]; then
command_print+=("set terminal postscript enhanced color font 'lmodern,10'
set output '$gnuplot_plot/$plot.eps'
replot
");fi

if [ $gnuplot_output == epslatex ]; then
command_print+=("set terminal epslatex size 8.0,4.0 color colortext
unset title
set output '$gnuplot_plot/$plot.tex'
replot
");fi

if [ $gnuplot_output == svg ]; then
command_print+=("set terminal svg size 800,400 fname 'lmodern' fsize 10
unset title
set output '$gnuplot_plot/$plot.svg'
replot
");fi

if [ $gnuplot_output == table ]; then
command_print+=("set table
set output '$gnuplot_plot/$plot.table'
replot
");fi

else
command_print="#set terminal svg size 800,400 fname 'lmodern' fsize 10
#set output '$gnuplot_plot/$plot.svg'
#replot"
fi

if [ -z $gnuplot_x_range ]; then
xrange="#set xrange [ 0 : 40 ]"
else
x_min=`echo $gnuplot_x_range | awk -F"-" '{print $1}'`
x_max=`echo $gnuplot_x_range | awk -F"-" '{print $2}'`
xrange="set xrange [ $x_min : $x_max ]"
fi

# Writing Gnuplot script and parameters
cat <<EOF > $gnuplot_temp/plotfile.gnu
set terminal wxt size 800,400 persist
set autoscale                                    # scale axes automatically
$xrange                                          # scale x-axis
unset log                                        # remove any log-scaling
unset label                                      # remove any previous labels
set xtic auto                                    # set xtics automatically
set ytic auto                                    # set ytics automatically

set title "$plot_title"
set xlabel "$xlabel"
set ylabel "$ylabel"
set key outside left bottom horizontal center Left reverse
set key noinvert samplen 0.2 spacing 1 width 0 height 2
unset colorbox
# MATLAB jet color palette
# line styles
set style line 1 lt 1 lc rgb '#000080' #
set style line 2 lt 1 lc rgb '#0000ff' #
set style line 3 lt 1 lc rgb '#0080ff' #
set style line 4 lt 1 lc rgb '#00ffff' #
set style line 5 lt 1 lc rgb '#80ff80' #
set style line 6 lt 1 lc rgb '#ffff00' #
set style line 7 lt 1 lc rgb '#ff8000' #
set style line 8 lt 1 lc rgb '#ff0000' #
set style line 9 lt 1 lc rgb '#800000' #
# palette
set palette defined (0 0.0 0.0 0.5, 1 0.0 0.0 1.0, 2 0.0 0.5 1.0, 3 0.0 1.0 1.0, \
4 0.5 1.0 0.5, 5 1.0 1.0 0.0, 6 1.0 0.5 0.0, 7 1.0 0.0 0.0, 8 0.5 0.0 0.0 )

plot $cmd

set table
set output '$gnuplot_temp/data-table.table'
replot

unset table
replot

#In case for building an output file
${command_print[@]}
EOF


# Writing Gnuplot plotline
cat <<EOF > $gnuplot_temp/plotline.gnu
plot $cmd
EOF


# Run plot command
case $gnuplot_output in
edit ) gnuplot $gnuplot_temp/plotfile.gnu - ; ;;
no-edit | eps | epslatex | svg | table ) gnuplot $gnuplot_temp/plotfile.gnu ; ;;
esac

echo -e "\nScript finished\n\n"
unset IFS

The reference

To view the code, you will need to expand the box below.

# File name: data-plot-ref.dat
# Reference file. Comparison between nuclide names and Z-number
#
#X    Name        Z     N    A
#Z    X    Name
g    gamma        0    0    0
n    neutron     0     1    1
p    proton      1     0    1
d    deutron     1     1    2
t    triton      1     2    3
h    helion      2     1    3
a    alpha       2     2    4
001     H      Hydrogen
002     He     Helium
003     Li     Lithium
004     Be     Beryllium
005     B      Boron
006     C      Carbon
007     N      Nitrogen
008     O      Oxygen
009     F      Fluorine
010     Ne     Neon
011     Na     Sodium
012     Mg     Magnesium
013     Al     Aluminium
014     Si     Silicon
015     P      Phosphorus
016     S      Sulfur
017     Cl     Chlorine
018     Ar     Argon
019     K      Potassium
020     Ca     Calcium
021     Sc     Scandium
022     Ti     Titanium
023     V      Vanadium
024     Cr     Chromium
025     Mn     Manganese
026     Fe     Iron
027     Co     Cobalt
028     Ni     Nickel
029     Cu     Copper
030     Zn     Zinc
031     Ga     Gallium
032     Ge     Germanium
033     As     Arsenic
034     Se     Selenium
035     Br     Bromine
036     Kr     Krypton
037     Rb     Rubidium
038     Sr     Strontium
039     Y      Yttrium
040     Zr     Zirconium
041     Nb     Niobium
042     Mo     Molybdenum
043     Tc     Technetium
044     Ru     Ruthenium
045     Rh     Rhodium
046     Pd     Palladium
047     Ag     Silver
048     Cd     Cadmium
049     In     Indium
050     Sn     Tin
051     Sb     Antimony
052     Te     Tellurium
053     I      Iodine
054     Xe     Xenon
055     Cs     Caesium
056     Ba     Barium
057     La     Lanthanum
058     Ce     Cerium
059     Pr     Praseodymium
060     Nd     Neodymium
061     Pm     Promethium
062     Sm     Samarium
063     Eu     Europium
064     Gd     Gadolinium
065     Tb     Terbium
066     Dy     Dysprosium
067     Ho     Holmium
068     Er     Erbium
069     Tm     Thulium
070     Yb     Ytterbium
071     Lu     Lutetium
072     Hf     Hafnium
073     Ta     Tantalum
074     W      Tungsten
075     Re     Rhenium
076     Os     Osmium
077     Ir     Iridium
078     Pt     Platinum
079     Au     Gold
080     Hg     Mercury
081     Tl     Thallium
082     Pb     Lead
083     Bi     Bismuth
084     Po     Polonium
085     At     Astatine
086     Rn     Radon
087     Fr     Francium
088     Ra     Radium
089     Ac     Actinium
090     Th     Thorium
091     Pa     Protactinium
092     U      Uranium
093     Np     Neptunium
094     Pu     Plutonium
095     Am     Americium
096     Cm     Curium
097     Bk     Berkelium
098     Cf     Californium
099     Es     Einsteinium
100     Fm     Fermium
101     Md     Mendelevium
102     No     Nobelium
103     Lr     Lawrencium
104     Rf     Rutherfordium
105     Db     Dubnium
106     Sg     Seaborgium
107     Bh     Bohrium
108     Hs     Hassium
109     Mt     Meitnerium
110     Ds     Darmstadtium
111     Rg     Roentgenium
112     Cn     Copernicium
113     Uut Ununtrium
114     Uuq Ununquadium
115     Uup Ununpentium
116     Uuh Ununhexium
118     Uuo Ununoctium