#!/bin/tcsh
#
# conc3 [speciation_list [element [start_concentration]]]
#
# This C-Shell script processes output files from
# "scan_conc", "scan_log", "scan_ph", "scan_t", or "scan_eh"
# EQ3NR run series or from MINTEQA2 runs processed with "minsweep".
#
# The final output generated by this script is a table with:
#  - the 1st column giving the stepping value,
#  - the 2nd column giving the pH,
#  - the 3rd column giving the true ionic strength,
#  - the 4th column giving the stoichiometric ionic strength,
#  - all following columns giving the real concentration of the species
#                          or the mineral saturation indices.
# 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 processing of solids is triggered by the element specification.
# If the element is followed by a "+", also the minerals and their saturation
# indices log (IP / Ksp) are extracted. If log (IP / Ksp) becomes positive, the
# mineral is oversaturated and should precipitate if not kinetically
# hindered. The list of minerals containing the requested element is
# built from a global list created by the C-Shell script "minerals".
#
# If the element is followed by a "-", also the total concentration of the
# respective strict basis species is extracted, useful e.g. in studies of
# the pH dependence of mineral solubilities.
#
# If the element is followed by a "=", also the partial pressures of all
# gases are extracted.
#
# If "Q" is specified as (fake) element, all mineral indices (and only them) are
# listed.
#
# The name of the output file is the same as the corresponding input file,
# but with the extension ".[3M]c" instead of ".[3M]l".
# 
# The output from "conc3" can be post-processed with the program "tab3"
# to create species distribution tables.
#
# The executable 'list2table' must be in the $PATH.
#
# Author: Dr.Vinzenz Brendler, FZ Dresden-Rossendorf e.V.
#                              Institute of Radiochemistry
#
# Version: 2.70                Date: 19-Aug-2010

onintr FINISH
stty icanon

if ($#argv < 1) then
	echo " "
	ls *.[3M]l >& /dev/null && ls *.[3M]l
	echo " "
	echo "Specify the name of the speciation list to be processed: "
	echo " "
	echo -n "     < "
	set argv = (`gets`)
endif
if ($#argv > 3) then
	echo "Error: No more than three arguments are allowed: "
	echo "   speciation_list [ element [ start_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 FILE = $argv[1]			# Speciation List to be Processed
set TEMP = $FILE:r.$$		# Temporary Output File

# Determine which speciation code was used
# and if the output file is already existent

set CODE = `sed -n '1s/^#\([36M]\).*/\1/p' $FILE`
set CONC = $FILE:r.${CODE}c	# Table of Concentrations

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/files/team/eq36/data
set SOLIDS = $EQ36DA/minerals.list
set TMIN   = $FILE:r.min		# Temporary Mineral List

# Check for request for mineral saturation indices.

# Case 1: Build up a list of relevant minerals for a specifc element

set SATURATION = 0
if ( $ELEMENT =~ *+ ) then
	set SATURATION = 1
	set ELEMENT = (`echo $ELEMENT | sed 's/+//' `)
	awk '/	'$ELEMENT'([^abcdefghijklmnopqrstuvwxyz]|$)/ { print $1 "(s)" }' $SOLIDS > $TMIN
	set NMIN = `cat $TMIN | wc -l`
	# vi $TMIN
endif

# Case 2: Build up a list of all relevant minerals

if ( $ELEMENT =~ Q ) then
	set SATURATION = 1
	awk '{ print $1 "(s)" }' $SOLIDS > $TMIN
	set NMIN = `cat $TMIN | wc -l`
	set TOTAL = 1e-5
	# vi $TMIN
endif

# Check for request for solubility calculations

set SOLUBILITY = 0
if ( $ELEMENT =~ *- ) then
	set SOLUBILITY = 1
	set ELEMENT = (`echo $ELEMENT | sed 's/-//' `)
endif

# Check for request for gas partial pressures

set GASES = 0
if ( $ELEMENT =~ *= ) then
	set GASES = 1
	set ELEMENT = (`echo $ELEMENT | sed 's/=//' `)
endif

# Specify the total analytical concentration of the above element,
# reject any input other than positive numbers.

if ($#argv > 2) then
	set TOTAL = $argv[3]
else
	if ( $ELEMENT =~ Q ) then
		set TOTAL = 1e-5
	else
		echo " "
		echo -n "Total analytical start concentration of $ELEMENT : "
		set TOTAL = (`gets`)
		if ( ! $#TOTAL) exit
		echo " "
	endif
endif
set TOTAL = `awk 'BEGIN { print (0.0 + '$TOTAL') }'`
if ( $TOTAL == 0 || $TOTAL =~ -* ) exit

# Extract only the really interesting lines from collected data:
# stepping value, pH, ionic strengths, saturation index (if requested),
# gas fugacities, and those aqueous species containing $ELEMENT
# and with the concentration exceeding 0.01 % of the analytical total.

awk 'BEGIN { LIMIT = (0.0001 * '$TOTAL') }\
	/^#/				{ print ; next } \
	/^pH /				{ print ; next } \
	/\(s\)/				{ if ('$SATURATION' == 1) print ; next } \
	/\(g\)/				{ if ('$GASES' == 1) print ; next } \
	/^[TS]I /			{ print ; next } \
	/'$ELEMENT'.*\(t\)/	{ if ('$SOLUBILITY' == 1) print ; next } \
	/'$ELEMENT'[^abcdefghijklmnopqrstuvwxyz]/	{ if ($2 + 0 > LIMIT) print } \
' $FILE > $TEMP

# In case of request for mineral saturation indices,
# skip all irrelevant minerals

if ($SATURATION == 1) then
	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
endif

# Transform the 1-column list into a table (with a seperate C-program).
# A maximum of 1000 iterations and 100 species will 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 | sed 's/ /	/g' > $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) tab3 $CONC $ELEMENT $TOTAL

# Delete all temporary files.

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

exit

# Future enhancements:
#
# - Enable correct handling of MINTEQA2 mineral names
# - What happens if both "+" and "-" are used after an element's specification ?