#!/bin/tcsh
#
# scan_conc [[EQ3/6 database] EQ3NR input file]
# scan_log  [[EQ3/6 database] EQ3NR input file]
# scan_gas  [[EQ3/6 database] EQ3NR input file]
# scan_ph   [[EQ3/6 database] EQ3NR input file]
# scan_eh   [[EQ3/6 database] EQ3NR input file]
# scan_t    [[EQ3/6 database] EQ3NR input file]
#
# This script performs several "single point" calculations with EQ3NR to
# simulate titration curves. It is an alias to scan_conc, scan_log,
# scan_gas, scan_eh, scan_ph and scan_t.
# According to the command line name, either the concentration of
# one component (linear or logarithmic), the gas partial pressure,
# the redox potential, the pH or the temperature is varied.
#
# Possible precipitates are announced, but not taken into account for
# the mass balances.
#
# The final output is a long list only containing
# the actual value of the varied variable, the pH, the true and the
# stoichiometric ionic strength, the concentration of all dissolved
# species and the total concentration of all strict basis species.
# Solids are included, if they are listed in the mineral saturation
# output block. Also all gases are included, provided they end with (g).
#
# Attention! All automatic input substitutions only work for input files
# using the same key words as in the examples from the EQ3/6 manual !
# And moreover, the input file needs to be in the W-format and is not
# allowed to contain several concatenated problem descriptions !
#
# The name of the output file is similar to the input file,
# but with the extension ".3l" instead of ".3i".
# The output can be post-processed for various element distributions with
# another csh script: "conc3"
#
# Author: Dr.Vinzenz Brendler, FZ Dresden-Rossendorf e.V.
#                              Institute of Radiochemistry
#
# Version: 2.00                Date: 31-Aug-2010

onintr FINISH
stty icanon

setenv PATH $PATH\:$EQ36CO		# include EQ3/6 binaries
								# (and also 'gets' on some systems)
set TDDB = fzr					# Default thermodynamic database

# Check for calling arguments

while ($#argv < 1 || $#argv > 2)
	echo ""
	echo "Error: Only ONE or TWO arguments are allowed:"
	echo "       an EQ3/6 thermodynamic database (optional, default is $TDDB)"
	echo "       and an EQ3NR input description file *.3i"
	echo " "
	ls *.3i >& /dev/null && ls *.3i
	echo " "
	echo    "Specify the name of the EQ3/6 database (optional)"
	echo -n "and select an EQ3NR input file from the above list: "
	set argv = (`gets`)
	echo " "
end

# Set internal variables depending on the number of arguments

set CALL = $0
set SCAN = $CALL:t				# Variable to be scanned
if ($#argv == 2) then
	set DATA = $1				# EQ3/6 data base to be used
	set FILE = $2				# EQ3NR Input File to be processed
else
	set DATA = $TDDB	 		# Default EQ3/6 data base
	set FILE = $1				# EQ3NR Input File to be processed
	echo " "
	echo "The scan will be performed using the default data base: $DATA ! "
endif

# Check for input file existence / read permission

if (! -r $FILE) then
	echo ""
	echo "Error: Input file $FILE does not exist or is not readable! "
	exit
endif

# Set temporary files and output

set NOCOM = scan$$.3n			# EQ3NR Input File without Comments
set OUT   = scan$$.3o			# Temporary EQ3NR-Output-File
set INPUT = scan$$.3i			# Temporary EQ3NR-Input-File
set PICK  = scan$$.3p			# Temporary EQ3NR-Pickup-File
set LIST  = $FILE:r.3l			# Output List of Concentrations
set DUMMY = /dev/null			# Dummy input for awk

# Check if output file already exists

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

# Identify the kind of titration requested by the user
# and set variables appropriately

switch ($SCAN)
case "scan_log":
		# scan a component concentration (logarithmic grid)
	set LOG
case "scan_conc":
		# scan a component concentration
	echo " "
	echo -n "Component with variable concentration: "
	set CVC   = (`gets`)
	if ($#CVC < 1) then
		echo "Invalid component formula"
		exit
	endif
	set RANGE = "concentration range in mol/L"
	set VAR   = ($CVC)
	breaksw
case "scan_gas":
		# scan a gas component partial pressure
	echo " "
	echo -n "Component with variable partial pressure: "
	set CVC   = (`gets`)
	if ($#CVC < 1) then
		echo "Invalid component formula"
		exit
	endif
	set RANGE = "pressure range in log(atm)"
	set VAR   = ($CVC)
	breaksw
case "scan_ph":
		# scan the pH
	set CVC   = "H+"
	set RANGE = "pH range"
	set VAR   = "pH"
	breaksw
case "scan_t":
		# scan the temperature
	set RANGE = "T range in degree Celsius"
	set VAR   = "T"
	breaksw
case "scan_eh":
		# scan the redox potential
	set REDOX = `awk '/iopt1-10=/ {print $2}' $FILE`
	if ($REDOX != -1) then
		echo " "
		echo "Input file can not be used for redox scan"
		echo "Set iopt1 to -1 (Now: $REDOX)"
		exit
	endif
	set RANGE = "redox potential range in V"
	set VAR   = "Eh"
	breaksw
endsw

# Check for presence of selected basic species

if ($?CVC) then
	grep -q -e "^data file master species= $CVC" $FILE
	if($status != 0) then
		echo "$CVC is not present in the input file"
		exit
	endif
endif

# Input of the range of the scan variable:
# start value, end value, step number or step size / scan increment

echo " "
echo "Input of $RANGE (start, end, step)"
echo -n "Start value : "
set X1 = (`gets`)
if ($X1 == "") then
	echo "Invalid start value"
	exit
endif
echo -n "  End value : "
set X2 = (`gets`)
if ($X2 == "") then
	echo "Invalid end value"
	exit
endif
if ($X2 == $X1) then
	echo "Identical start and end value"
	exit
endif
if ($?LOG) then
	set SPC = "L"		# logarithmic spacing of scan variable
	echo -n "  Number of steps : "
	set MAXLOOP = (`gets`)
	if ($MAXLOOP == "" || $MAXLOOP < 2) then
		echo "Invalid number of steps"
		exit
	endif
	set DELTA = `awk 'BEGIN { printf "%11.5e\n" , exp((log('$X2') - log('$X1'))/'$MAXLOOP') }' $DUMMY`
else
	set SPC = "E"		# equal spacing of scan variable
	echo -n "  Step size : "
	set DELTA = (`gets`)
	if ($DELTA == "" || $DELTA == 0) then
		echo "Invalid step size"
		exit
	endif
	set MAXLOOP = `awk 'BEGIN { printf "%2.0f\n" , ('$X2' - '$X1')/'$DELTA' }' $DUMMY`
endif
echo " "

# Remove all comment lines from the input file to make later parsing
# and substituting unambigious

sed '/^\*/d' $FILE > $NOCOM

# Reformat start value of scan variable

set XNOW = `awk 'BEGIN { printf "%11.5e\n" , '$X1' }' $DUMMY`

# Define scan type identifier for output file
# This identifier helps further processing of the output and follows a strict
# syntax: #[3|6|M][E|L]#[pH|Eh|T|$CVC]
# The first character specifies the code ran: EQ3NR, EQ6 or MINTEQA2
# The second character specifies the scan spacing: equal or logarithmic.
# The last part gives the scanned variable: pH, redox potential, temperature or
# a species ($CVC) with variable concentration or partial pressure.

set TOKEN = "#3"$SPC"#"$VAR[1]

# Initialize scan loop

@ MAXLOOP++
@ COUNTER = 0
cat $DUMMY > $LIST

# Start loop over the scanned variable

while ($COUNTER < $MAXLOOP)

	# Build a modified copy of the original input file and run EQ3NR

	switch ($SCAN)
	case "scan_conc":
	case "scan_gas":
	case "scan_log":
		sed "/master .*= $CVC/,/csp=/s/csp= .*/csp= $XNOW/" $NOCOM > $INPUT
		breaksw
	case "scan_ph":
		sed "/master .*= $CVC/,/csp=/s/csp= .*/csp= -$XNOW/" $NOCOM > $INPUT
		breaksw
	case "scan_t":
		sed "s/tempc= .*/tempc=  $XNOW/" $NOCOM > $INPUT
		breaksw
	case "scan_eh":
		# Eh is the only scan that can produce negative variables
		# Then the minus sign will increase the string length
		if ($XNOW =~ -[0-9].* ) then
			sed "s/ fep= ............/ fep= $XNOW/" $NOCOM > $INPUT
		else
			sed "s/ fep= ............/ fep=  $XNOW/" $NOCOM > $INPUT
		endif
		breaksw
	endsw

	echo " "
	echo "++++++++++++++ $VAR now at $XNOW ++++++++++++"
	echo " "

	# vi $INPUT
	runeq3 $DATA $INPUT

	# Add the actual value of the scan variable to the collection list
	# preceeded by a tag indicating the type of the scanned variable

	echo $TOKEN $XNOW >> $LIST

	# Extract the following information:
	# - pH,
	# - true and stoichiometric ionic strength,
	# - concentrations of all aqueous species,
	# - total concentration of all basic components,
	# - all supersaturated solids and
	# - all gas fugacities

	awk '\
		/^ *modified NBS pH scale/		{print $3 , $5} \
		/^ *True ionic strength=/		{print "TI" , $4} \
		/^ *Stoi.* ionic strength=/		{print "SI" , $4} \
		/^ *Species *Molality/,/ --- M/	{print $1 , $2} \
		/^ *Mineral  *Log/,/ approx/	{if ($2 > 0.0) print $1 "(s)" , $2} \
		/--- Current/,/Single ion/		{if ($2 + 0.0 > 0.0) print $1 "(t)" , $2} \
		/^ *Gas  *Fugacity/,/ --- Reading/	{if ($2 + 0.0 > 0.0) print $1 , $2} \
	' $OUT | sed -n '/approx/d \
		/[0-9]/p' >> $LIST
	#' | sed '/approx/d' |\
	#sed -n '/[0-9]/p' >> $LIST

	# Increment scan variable and loop counter

	if ($?LOG) then
		set XNOW = `awk 'BEGIN { printf "%11.5e\n" , ('$XNOW' * '$DELTA') }' $DUMMY`
	else
		set XNOW = `awk 'BEGIN { printf "%11.5e\n" , ('$XNOW' + '$DELTA') }' $DUMMY`
	endif
	@ COUNTER++

end
echo " "
echo "        Successfully created: $LIST"

# Ask the user if further data treatment, namely
# creating a concentration table is wanted.

echo " "
echo -n "Continue with creating a concentration table y/(n) ? "
set rr = (`gets`)
echo " "
if ($rr == y || $rr == Y) conc3 $LIST

# Delete all intermediate EQ3NR files

FINISH:
/bin/rm -f $OUT $INPUT $PICK $NOCOM

# Future enhancements planned:
#
# - Restructuring, putting the scanned variable as 1st calling argument,
#   dropping the cmd-aliases.
#
# - Mark the value to be varied already in the input template with a
#   special token, makes later exchange easier
#
# - 2D Scans
#
# - Incorporation of the runeq36 scripts to avoid a lot of overhead by
#   invoking these scripts
#
# - Better parsing of the input files
#