#!/bin/tcsh
#
# tab6 [conc6_output [element [start_concentration [added_concentration]]]]
#
# This C-Shell script processes output files from "conc6",
# resulting from an EQ6 run.
#
# The final output generated by this script is a table with:
#  - the 1st column giving the X-Axis stepping variable (Xi, pH, or Eh),
#  - all following columns giving the element speciation (in %).
# It is necessary to restrict the species to only those containing a
# given (user-defined) element. The stoichiometric coefficents of this
# element in each species are not automatically determined,
# so the user has to provide them.
#
# The name of the output file is the same as its corresponding input file,
# but instead of the extension ".6c" it has ".6d".
#
# The output can be post-processed to create graphs on the screen with
# the program 'nplot' (default) or 'mplot'.
#
# Author: Dr.Vinzenz Brendler, Forschungszentrum Dresden-Rossendorf e.V.
#                              Institute of Radiochemistry
#
# Version: 0.96                Date: 26-jan-2011

onintr FINISH
stty icanon

if ($#argv < 1) then
	echo " "
	ls *.6c >& /dev/null && ls *.6c
	echo " "
	echo "Specify the name of the concentration file: "
	echo " "
	echo -n "     < "
	set argv = (`gets`)
endif
if ($#argv > 4) then
	echo "Error: No more than four arguments are allowed: "
	echo "   conc6_output [ element [ start_concentration [ added_concentration ] ] ]"
	exit
endif
if (! -r $argv[1]) then
	echo "Error: Input file $argv[1] is not existent or not readable! "
	exit
endif
if (-z $argv[1]) then
	echo "Error: Input file $argv[1] is empty! "
	exit
endif

set CONC = $argv[1]			# File with Concentration Table
set DIST = $CONC:r.6d		# Table of Species Distribution (Percentages)

if (-e $DIST) then
	echo -n "Output file $DIST already exists, overwrite y/(n) ? "
	set rr = (`gets`)
	if ($rr != y && $rr != Y) exit
endif

# Specify the element for which the species concentration is looked for

set SYMBOL = U		# Uranium as Default Species Selected for Distribution

if ($#argv > 1) then
	set ELEMENT = $argv[2]
else
	echo " "
	echo -n "Element wanted for distribution ( <CR> = $SYMBOL ) : "
	set ELEMENT = (`gets`)
	if ( ! $#ELEMENT) set ELEMENT = $SYMBOL
endif

# Specify the total start analytical concentration of the above element

if ($#argv > 2) then
	set TOTAL = $argv[3]
else
	echo -n "Total analytical start concentration of $ELEMENT : "
	set TOTAL = (`gets`)
	if ( ! $#TOTAL) exit
	echo " "
endif

# Specify the added analytical concentration of the above element

if ($#argv > 3) then
	set ADDED = $argv[4]
else
	set ADDED = 0.0
endif

# Dimension the vector of the stoichiometric coefficents and ask for them

set COEFFS = `head -1 $CONC`
echo "Stoichiometric coeffcients of $ELEMENT ( <CR> = 1 ) in"
@ COUNTER = 0
foreach entry ($COEFFS)
	@ COUNTER++
	if ($COUNTER < 5) continue
	echo -n "   $entry : "
	set rr = (`gets`)
	if ( ! $#rr) set rr = 1
	if ($rr < 0) set rr = 0
	set COEFFS[$COUNTER] = $rr
end
echo " "

# Ask for the variable to be used as X-axis (default: pH)

echo "Please specify the variable to be used as X-axis (default: pH)"
echo " 1  ....... reaction progress Xi"
echo "(2) ....... pH"
echo " 3  ....... redox potential Eh"
echo -n " > "
set XCOL = (`gets`)
if ( (! $#XCOL) || ($XCOL < 1) || ($XCOL > 3)) set XCOL = 2
echo " "
switch ($XCOL)
case "1":
		set XAXIS = "Xi"; set DELTA = 0.001; set XOUT = "%10.3e"; breaksw
case "2":
		set XAXIS = "pH"; set DELTA = 0.01;  set XOUT = "%8.2f"; breaksw
case "3":
		set XAXIS = "Eh"; set DELTA = 0.001; set XOUT = "%8.3f"; breaksw
endsw

# Ask for the grid size for the X axis (to reduce table size)

echo -n "Smallest grid size for $XAXIS ( <CR> = $DELTA ) : "
set GRID = (`gets`)
if ( ! $#GRID) set GRID = $DELTA
echo " "

# Normalize to the total analytical concentration of the wanted element
# and multiply with stoichiometric coefficents.
# Table entries differing from the previous line in the x axis values
# less than the grid size are not printed.
# Change all space sequences into tabs

echo $COEFFS | cat - $CONC | \
awk 'BEGIN  {	Cstart = '$TOTAL' \
				Cadded = '$ADDED' \
				LAST = 0.0 \
				LIMIT = '$GRID' * '$GRID' }\
	NR == 1 {	for (i = 5; i <= NF; i++) coeff[i] = $i * 100.0 ; next }\
	NR == 2 {	printf "#6E#%s", "'$XAXIS'" \
				for (i = 5; i <= NF; i++) { \
					solid[i] = index($i, "(s)") \
					printf " %s", $i \
				} \
				printf "\n"; next } \
		    {	DELTA = ($'$XCOL' - LAST) ; if (DELTA * DELTA < LIMIT) next }\
		    {	Xi = $1 ; LAST = $'$XCOL' ; printf "'$XOUT'", LAST \
		    	Vsol = 1.0 + Xi \
		    	Ctotal = (Cstart + Xi * Cadded) / Vsol \
				for (i = 5; i <= NF; i++) { \
					if (solid[i] > 0) $i = $i / Vsol \
					printf " %8.2f", $i / Ctotal * coeff[i] \
				} \
				printf "\n" } \
' | sed 's/  */	/g' > $DIST && echo "        Successfully created: $DIST"

# Ask the user if further data treatment, namely graph preparation
# (gnuplot & postscript) is wanted

echo " "
echo -n "Continue with the plot routine y/(n) ? "
set rr = (`gets`)
echo " "
if ($rr == y || $rr == Y) nplot $DIST $ELEMENT

# Delete all temporary files

FINISH:

exit

# Future enhancements:
#