#!/bin/tcsh
#
# conc6 [EQ6_output [element]]
#
# This C-Shell Script extracts the concentration information from
# an EQ6 output file.
#
# The final output generated by this script is a table with:
#  - all species names in the 1st row and from then onwards
#  - the 1st column giving the reaction progress Xi,
#  - the 2nd column giving the pH,
#  - the 3rd column giving the Eh,
#  - the 4th column giving the ionic strength,
#  - all following columns giving the real concentration of the species
#    (molality for dissolved species and total moles for minerals).
# It is necessary to restrict the species to only those containing a
# given (user-defined) element. Only species with concentrations larger
# than 0.01 % of the total analytical start concentration are considered.
#
# The list of minerals containing the requested element is
# built from a global list created by the C-Shell script "minerals".
#
# The name of the output file is the same as the corresponding EQ6 file,
# but with the extension ".6c" instead of ".6o".
#
# The output from "conc6" can be post-processed with the program "tab6"
# to create species distribution tables.
#
# The executable 'list2table' must be in the PATH.
#
# Author: Dr.Vinzenz Brendler, Forschungszentrum Dresden-Rossendorf e.V.
#                              Institute of Radiochemistry
#
# Version: 2.12                Date: 07-feb-2011

onintr FINISH
stty icanon

if ($#argv < 1) then
	echo " "
	ls *.6o >& /dev/null && ls *.6o
	echo " "
	echo "Specify the name of the EQ6 output file to be processed: "
	echo " "
	echo -n "     < "
	set argv = (`gets`)
endif
if ($#argv > 2) then
	echo "Error: No more than two arguments are allowed: "
	echo "   input_file [ elementaselection ]"
	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 FILE = $argv[1]			# EQ6 - Output File to be Processed
set CONC = $FILE:r.6c		# Table of Concentrations
set TEMP = $FILE:r.$$		# Temporary Output File

if (-e $CONC) then
	echo -n "Output file $CONC 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

# Global list of mineral stoichiometry created by the 'minerals' script

if (! $?EQ36DA) set EQ36DA = /net/filet/team/eq36/data
set SOLIDS = $EQ36DA/minerals.list
set TMIN   = $FILE:r.min		# Temporary Mineral List

# Build a list of relevant minerals

awk '/	'$ELEMENT'([^abcdefghijklmnopqrstuvwxyz]|$)/ { print $1 "(s)" }' $SOLIDS > $TMIN
set NMIN = `cat $TMIN | wc -l`
echo "Number of relevant minerals : " $NMIN

# Extract the total analytical concentration of the above element
# in the aqueous phase only

set TOTAL = `awk '/Component/,/input./ {if ($1 == "'$ELEMENT'") print $2}' $FILE`
if ( ! $#TOTAL) set TOTAL = 0.0
echo "Total concentration of $ELEMENT : " $TOTAL

# Check if the EQ6 run was a titration with a solution
# If so, extract concentration of wanted element in the titrant

set TITRATION = `grep nmodl1= $FILE | awk '{print $2}'`
if ($TITRATION == 1) then
	set ADD = `awk '/reactant=/,/endit./ {if ($1 == "'$ELEMENT'") print $2}' $FILE`
	echo "Added concentration of $ELEMENT : " $ADD
else
	set ADD = 0
endif

# Extraction of useful information from EQ6 output file.
# All empty lines and all lines before the 1st stepping are discarded.
# The following information is collected:
#  - Reaction progress Xi
#  - pH
#  - Eh
#  - Ionic strength
#  - all molalities of species containing $ELEMENT and exceeding
#    0.01 % of the total aqueous concentration of $ELEMENT
#  - total moles of all minerals

sed '\
/^$/d\
1,/^  *Stepping to zi/d' $FILE |\
awk 'BEGIN { LIMIT = (0.0001 * '$TOTAL') }\
	/^ *Reaction/         { print "#6#Xi" , $4 }\
	/^ *modified NBS pH/  { print "pH" , $5 ; print "Eh" , $6 }\
	/^ *Ionic strength =/ { print "IS" , $4 }\
	/^ *Species/,/ Mass/  { if (($2 + 0 > LIMIT) && ($1 ~ /'$ELEMENT'[^abcdefghijklmnopqrstuvwxyz]/))\
							print $1 , $2 }\
	/^ *Product/,/Mass, grams/  { if ($3 + 0 > 0) print $1 "(s)" , $3 }\
' > $TEMP

# echo -n "CR >> " ; gets
# vi $TEMP
# Skip all irrelevant minerals

cat $TMIN $TEMP | awk 'BEGIN { nlist = '$NMIN' } \
	NR <= nlist { mineral[NR] = $1 ; next } \
	/\(s\)/	{ for (i=0; i<nlist; i++) if ($1 == mineral[i]) print ; next } \
	{ print } \
' > $CONC
mv $CONC $TEMP
# echo -n "CR >> " ; gets
# vi $TEMP

# Transform the 1-column list into a table (with a seperate C-program).
# A maximum of 1000 iterations and 100 species can be processed,
# this restriction comes from the source code of 'list2table'

set NITER = `grep "^#" $TEMP | wc -l`
if ($NITER > 1000) then
	echo "Too many ($NITER) parameter steps, 1000 is the maximum ! "
endif
list2table $TEMP > $CONC && (echo "        Successfully created: $CONC")

# Ask the user if further data treatment, namely multiplying with
# stoichiometric coefficients and creating a distribution table is wanted

echo " "
echo -n "Continue with creating a distribution table y/(n) ? "
set rr = (`gets`)
echo " "
if ($rr == y || $rr == Y) tab6 $CONC $ELEMENT $TOTAL $ADD

# Delete all temporary files

FINISH:
/bin/rm -f $TEMP $TMIN

exit

# Future enhancements:
#