#! /bin/bash
# invoke gnuplot to compare the soil submodel transitions out of the soillay.out file
# from WEPS generated runs. (or a file with the same format)
# files are checked for matching years and then plotted for each year

# invocation: singlelayersoiltrans [laynum file-1 file-2 ... file-#]

# set up column and title correspondence
mincol=21
maxcol=28
#daysim idoy yr layer depth bszlyt bsdblk bseags bseagmn bseagm bseagmx bslagn bslmin bslagm bslmax bslagx bs0ags bsdagd rel_ag_stab rel_geo_mean freeze freeze_thaw thaw frozen wetting drying warm_puddling wet_bulk_den #daysim idoy yr layer depth bszlyt bsdblk bseags bseagmn bseagm bseagmx bslagn bslmin bslagm bslmax bslagx bs0ags bsdagd rel_ag_stab rel_geo_mean

coltitle[21]=`echo "21) Freeze Transition"`
coltitle[22]=`echo "22) Freeze-Thaw Trans"`
coltitle[23]=`echo "23) Thaw Transition"`
coltitle[24]=`echo "24) Soil Remains Frozen"`
coltitle[25]=`echo "25) Wetting Transition"`
coltitle[26]=`echo "26) Drying Transition"`
coltitle[27]=`echo "27) Warm Puddling Trans"`
coltitle[28]=`echo "28) Bulk Density Wetted"`

# grab layer number and shift command arguments
laynum=$1
shift

# test if it is a number
echo ${laynum} | grep [^0-9] > /dev/null 2>&1
if [ "$?" -eq "0" ]
then
    # laynum is not a number
    echo '['${laynum}'] is not a number'
    echo 'invocation: singlelayersoiltrans [laynum file-1 file-2 ... file-#]'
    echo 'Please enter a Layer Number as the first value.'
    exit
elif [[ ${laynum} -lt 1 ]]
then
    echo 'invocation: singlelayersoiltrans [laynum file-1 file-2 ... file-#]'
    echo 'Valid Layer Numbers are 1 to the maximum layers in file.'
    exit
fi

# loop through all file names
# extract specified layer from each soil layer file and create temporary files
let nfile=0
for nf in $*
do
    # create filename array
    let nfile=nfile+1
    filename[nfile]="$nf"
    selname[nfile]="selfile${nfile}"
    #echo ${filename[nfile]} ${selname[nfile]}

    # check for existence of input file
    if [ -f ${filename[nfile]} ]
    then
        # file exists
        echo "Extracting layer info for file: ${filename[nfile]}"

        # extract specified layer from file and put results in temporary file
        dm "if x4=${laynum} then INPUT else SKIP" < ${filename[nfile]} > ${selname[nfile]}

    else
        # file does not exist
        echo "File does not exist. Check the name: ${filename[nfile]}"
        echo "Continuing without this file."
        let nfile=nfile-1
    fi
done

# check file count
if [ $nfile -eq 0 ]
then
    echo "No valid file names. Check your command line and try again."
    exit
fi

echo ${coltitle[${colnum}]}

# zero out temp.dat file
echo '' > temp.dat

# loop through all files
for num in $(seq $nfile)
do
    # loop through all columns
    for col in $(seq 21 28)
    do
        # extract column of interest and save in temp file
        colex ${col} < ${selname[num]} > temp

        # find number of lines with values in file
        numval[col]=`stats n < temp`
        sumval[col]=`stats sum < temp`

        # find fraction of days with transitions
        fracval[col]=`echo "scale=5; ${sumval[col]}/${numval[col]}" | bc -l`
    done

    # print output line
    echo ${fracval[21]} ${fracval[22]} ${fracval[23]} ${fracval[24]} ${fracval[25]} ${fracval[26]} ${fracval[27]} ${fracval[28]} >> temp.dat

done

# remove the temp file used above
rm temp

# turn rows into columns
transpose < temp.dat > temp1.dat
mv temp1.dat temp.dat

# create gnuplot script to plot data
echo 'reset' > temp.plt

#--- start pdf ---
 echo "set terminal pdf" >> temp.plt
 echo "set output 'soil_trans.pdf'" >> temp.plt
 view_pause=0
#--- end pdf ---

#--- start screen ---
# view_pause=-1
#--- end screen ---

# create graph title line
echo "set title 'Fraction of days with soil state transitions, layer ${laynum}'" >> temp.plt

echo "set xtics rotate by -90 ( \\" >> temp.plt
# loop through all columns
for col in $(seq 21 27)
do
    adjcol=$(($col - 21))
    echo "'${coltitle[col]}' ${adjcol}, \\" >> temp.plt
done
adjcol=$((28-21))
echo "'${coltitle[28]}' ${adjcol} )" >> temp.plt

echo "set style fill pattern 4" >> temp.plt

# calculate box width and offset values
totalwidth="0.7"
boxwidth=`echo "scale=5; ${totalwidth}/${nfile}" | bc -l`
baseoffset=`echo "scale=5; (${boxwidth}/2)-(${totalwidth}/2)" | bc -l`

#echo "totalwidth $totalwidth"
#echo "boxwidth $boxwidth"
#echo "baseoffset $baseoffset"

# create set boxwidth line
echo "set boxwidth ${boxwidth}" >> temp.plt

# create plot line
echo "plot \\" >> temp.plt
# loop through all files
for num in $( seq $(($nfile - 1)) )
do
    # calculate box offset
    boxoffset=`echo "scale=5; ${baseoffset}+${boxwidth}*(${num}-1)" | bc -l`
    # output plot line
    echo "'temp.dat' using (\$0+(${boxoffset})):${num} with boxes title '${filename[num]}', \\" >> temp.plt
done
# calculate box offset
boxoffset=`echo "scale=5; ${baseoffset}+${boxwidth}*(${nfile}-1)" | bc -l`
# last line
echo "'temp.dat' using (\$0+(${boxoffset})):${nfile} with boxes title '${filename[nfile]}'" >> temp.plt

echo "pause ${view_pause}" >> temp.plt

gnuplot temp.plt

#rm temp.dat
#rm temp.plt