#! /bin/bash

# extract wind direction, total erosion from daily erosion events.
# requires erosion subumdel flag with second bit set and command line flag -e1.
# multiple runs will be compared. When erosion days do not match, a day with zero
#erosion values will be created. Most useful to compare runs on a single site with
#varying row direction.

# invocation: compare_daily [run_dir1 run_dir2 ... run_dirn]

# grab and count command arguments
let ncmdarg=0
for ncmd in $*
do
    let ncmdarg=ncmdarg+1
    cmdarg[ncmdarg]="$ncmd"
    # echo "Command argument $ncmdarg = ${cmdarg[ncmdarg]}"
done

# check command argument count, zero means none
if [[ ncmdarg -eq 0 ]]
then
    echo 'compare_daily [run_dir1 run_dir2 ... run_dirn]'
    exit
fi

# parse command arguments checking existence of specfied directories and required files
# initialize total number of run directories found
let totrundir=0
for ncmd in $( seq 1 $ncmdarg )
do
    # test if directory exists
    echo "${cmdarg[ncmd]}" | grep [^0-9] > /dev/null 2>&1
    if [ -d "${cmdarg[ncmd]}" ]
    then
        # directory exists
        # increment count
        let totrundir+=1

        # create comparison file name for simulation day and wind direction
        compare_file[totrundir]=$(echo $(basename ${cmdarg[ncmd]})_${totrundir}.dde)
        compare_name[totrundir]=$(echo $(basename ${cmdarg[ncmd]})_${totrundir})
        echo "File ${totrundir}: ${compare_file[totrundir]}"

        # Pull daysim and wind direction from erod.sweep
        find "${cmdarg[ncmd]}/sae_in_out_files/" -name "erod.sweep" -exec grep -E '<WEPS_SimulationDay>|<SCI_WindDirection>' {} \; | sed 's/<SCI_WindDirection>//' | sed 's/<\/SCI_WindDirection>//' | sed 's/<WEPS_SimulationDay>//' | sed 'N;s/<\/WEPS_SimulationDay>\n//g' | sort -g > temp.dwd

        # Pull daysim and total erosion from erod.egrd
        find "${cmdarg[ncmd]}/sae_in_out_files/" -name "erod.egrd" -exec grep -E '# WEPS erosion day mon yr daysim|repeat of total, salt/creep, susp, PM10, PM2.5' {} \; | sed 'N;s/\n//g' | sed -E 's/(\s+\S+){9}\s+//' | sed 's/#//' | sed -E 's/(\s+\S+){7}\s+/ /' | dm s1 s2 | sort -g > temp.dte

        # place into one file and process daysim, wind direction, total erosion
        # test for line count match between temp.dwd temp.dte
        temp1=$(cat temp.dwd | wc -l)
        temp2=$(cat temp.dte | wc -l)
        if [[ ${temp1} -eq ${temp2} ]]
        then
            abut temp.dwd temp.dte | dm x1 x2 x4 > ${compare_file[totrundir]}
        else
            echo "Different number of erod.sweep and erod.egrd daily files. exiting."
            exit
        fi
    else
        # directory does not exist
        echo "${cmdarg[ncmd]} directory does not exist. Exiting."
        exit
    fi
done

echo "Found ${totrundir} run directories"

# define the bins, month names and wind directions
#totbin=25
ulimit="1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5 25.5 30.5 35.5 40.5 45.5"
totpane=12
totdir=16
totdir1=15
direct=( 0 22.5 45 67.5 90 112.5 135 157.5 180 202.5 225 247.5 270 292.5 315 337.5 )
datvals=( nevents toterod avgerod )
# count, sumtotal erosion, and event average total erosion for each cardinal direction

# intialize scaling variables
maxnumevent=0
maxtoterode=0
maxavgerode=0

# loop through all input files
for nrundir in $( seq 1 ${totrundir} )
do
    echo "${nrundir} - ${compare_file[nrundir]}"

    numevent=( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 )
    toterode=( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 )
    avgerode=( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 )

    # read in file
    while read record
    do
        # select values for each wind direction
        fdir=$(echo "${record}" | cut -f2)

        for ndir in $( seq 0 ${totdir1} )
        do
            if [ "${fdir}" = "${direct[ndir]}" ]
            then
                let "numevent[ndir] += 1"
                if [ ${numevent[ndir]} -gt ${maxnumevent} ]
                then
                    maxnumevent=${numevent[ndir]}
                fi
                fterod=$(echo "${record}" | cut -f3)
                toterode[ndir]=$(echo "${toterode[ndir]} + ${fterod}" | bc)
                if (( $(echo "${toterode[ndir]} > ${maxtoterode}" | bc -l) ))
                then
                    maxtoterode=${toterode[ndir]}
                fi
                #echo "${direct[ndir]} :: ${numevent[ndir]} ; ${maxnumevent} : ${toterode[ndir]} ; ${maxtoterode}"
            fi
        done
    done < ${compare_file[nrundir]}

    echo "" > ${compare_file[nrundir]}.dat
    echo "" > ${compare_file[nrundir]}.dat
    for ndir in $( seq 0 ${totdir1} )
    do
        if [ ${numevent[ndir]} -gt 0 ]
        then
            avgerode[ndir]=$(echo "${toterode[ndir]} / ${numevent[ndir]}" | bc -l)
            if (( $(echo "${avgerode[ndir]} > ${maxavgerode}" | bc -l) ))
            then
                maxavgerode=${avgerode[ndir]}
            fi
        fi
        # write values to file for graphing
        echo "${direct[ndir]} ${numevent[ndir]} ${toterode[ndir]} ${avgerode[ndir]}" >> ${compare_file[nrundir]}.dat
    done

    echo "MAX :: ${maxnumevent} : ${maxtoterode} : ${maxavgerode}"

done

colmaxarr=( ${maxnumevent} ${maxtoterode} ${maxavgerode} )

# set up pane sizes for a 3 across, 4 down
tacross=3
tdown=4
wsize=$( echo "scale=4; 0.999/${tacross}" | bc -l )
hsize=$( echo "scale=4; 0.999/${tdown}" | bc -l )

# set up pane origins
let "eacross = $tacross - 1"
let "edown = $tdown - 1"
for nh in $( seq 0 $edown )
do
    for nw in $( seq 0 $eacross )
    do
        let "npane = $nh * $tacross + $nw"
        worigin[npane]=$( echo "scale=4; ${nw}*${wsize}" | bc -l )
        horigin[npane]=$( echo "scale=4; (${edown}-${nh})*${hsize}" | bc -l )

        echo "ORG: $nh : $nw : $npane : ${worigin[npane]} : ${horigin[npane]}"

    done
done

# create gnuplot scripts for wind direction day count, total erosion. event average erosion
plotfile="compare_daily.plt"
outfile="compare_daily.pdf"
echo "set terminal pdf font 'Helvetica,3'" > ${plotfile}
echo "set output '${outfile}'" >> ${plotfile}
echo "set size 1.0, 1.0" >> ${plotfile}
echo "set origin 0.0, 0.0" >> ${plotfile}

for nrundir in $( seq 1 ${totrundir} )
do
    let "endpane = totpane - 1"
    for colfile in $( seq 2 $((tacross+1)) )
    do
        colpane=$(( ((nrundir-1)*tacross + (colfile-2)) % totpane ))

        if [ ${colpane} -eq 0 ]
        then
            # this command necessary to display all plots on some terminals
            echo "set multiplot" >> ${plotfile}
        fi
        # w-width h-height
        echo "set size ${wsize}, ${hsize}" >> ${plotfile}
        echo "set origin ${worigin[colpane]}, ${horigin[colpane]}" >> ${plotfile}
        if [ ${colfile} -eq 2 ]
        then
            echo "set label 1 '${compare_name[nrundir]}' noenhanced at graph 0.5,1.0 center front font 'Helvetica,4'" >> ${plotfile}
        fi
        echo "set label 3 '${datvals[colfile-2]}' at graph 0.02,0.92 center front font 'Helvetica,4'" >> ${plotfile}
        # plot the wind rose
        echo "set size square" >> ${plotfile}
        echo "unset key" >> ${plotfile}
        echo "unset xlabel" >> ${plotfile}
        echo "unset ylabel" >> ${plotfile}
        echo "set angles degrees" >> ${plotfile}
        echo "set polar" >> ${plotfile}
        echo "set angles degrees" >> ${plotfile}
        echo "unset border" >> ${plotfile}

        echo "unset xtics" >> ${plotfile}
        echo "unset ytics" >> ${plotfile}
        echo "set rtics autofreq axis scale 0.5 font 'Helvetica,3'" >> ${plotfile}
        echo "unset raxis" >> ${plotfile}
        echo "set grid polar 22.5" >> ${plotfile}
        echo "set style line 10 dashtype (2,2,2,2) lc 0 lw 0.3" >> ${plotfile}
        echo "set grid ls 10" >> ${plotfile}
        echo "set size ratio 1" >> ${plotfile}
        echo "set rrange [0:${colmaxarr[colfile-2]}]" >> ${plotfile}
        echo "plot '${compare_file[nrundir]}.dat' index 0 using (90-\$1):${colfile} with impulse linewidth 1.25" >> ${plotfile}

        echo "INDEX: ${nrundir} : ${colfile} : ${colpane}"

        if [ ${colpane} -eq ${endpane} ]
        then
            # this command necessary to display all plots on some terminals
            echo "unset multiplot" >> ${plotfile}
        fi
    done

done

gnuplot ${plotfile}

rm temp.*
rm *.dde*

exit

# define scale for wind direction percentages
probmax=30

dirfile="alldir.out"
calmfile="allcalm.out"
thistfile="temp.out"
histfile="dirhist.out"
tmaxminfile="temp1.out"
maxminfile="allmaxminhr.out"
plotfile="allplot.plt"
outfile="windplot.pdf"

# initialize output files
echo "# Probability of wind from the sixteen cardinal directions, months are columns" > ${dirfile}
echo "# Probability of no wind, calm, by month, stations are columns" > ${calmfile}
let "dend = totdir - 1"
for ndir in $( seq 0 ${dend} )
do
    echo "# Histograms of wind velocity from the ${direct[ndir]} degree cardinal direction for each month, columns: 1-bins 2-13-months , etc, stations are separated by double blank lines" > ${ndir}${histfile}
done
echo "# month, maximum to minimum wind velocity ratio, hour of maximum" > ${maxminfile}
echo "# GNUPLOT commands to plot wind station data records" > ${plotfile}

# combine files into data files for plotting
for nsta in $( seq 1 ${totsta} )
do
    # initialize temporary histogram file
    echo "${ulimit}" > ${thistfile}
    echo "${tmonths}" > ${tmaxminfile}

    #initialize read counters
    linecnt=0
    dircnt=0
    histcnt=0

    # create double blank line to seperate indexes for gnuplot for separate stations
    echo "" >> ${dirfile}
    echo "" >> ${dirfile}
    echo "" >> ${calmfile}
    echo "" >> ${calmfile}
    echo "" >> ${maxminfile}
    echo "" >> ${maxminfile}
    echo "" >> ${ndir}${histfile}
    echo "" >> ${ndir}${histfile}

    # read in file
    while read record
    do
        let linecnt=linecnt+1
        #echo "Line: ${linecnt}"
        #echo "${record}"

        # check line number for appropriate action
        if [[ ${linecnt} -eq 1 ]]
        then
            # grab station name
            staname[nsta]=`echo ${record} | sed s/#//`
        fi
        if [[ ${linecnt} -ge 3 ]]
        then
            if [[ ${linecnt} -le 18 ]]
            then
                # these are the wind direction probabilies
                # directly output into file for gnuplot
                # prepend the direction labels
                echo "Reading wind direction probabilities for ${direct[dircnt]} degrees"
                echo "${direct[dircnt]} ${record}" >> ${dirfile}
                let "dircnt += 1" 
            fi
            if [[ ${linecnt} -eq 19 ]]
            then
                # probabilities for calm
                echo "Reading wind direction probabilities for calm"
                # create calm file for this station
                echo "-${probmax} ${record}" >> ${calmfile}
            fi
        fi
        if [[ ${linecnt} -ge 20 ]]
        then
            if [[ ${linecnt} -le 211 ]]
            then
                # these are the wind velocity histograms
                # use one dim array to simulate three dim, store rows sequentially
                # rows are nhist long (ie columns), with totdir rows, with twelve of these blocks for the months
                moncnt=`expr $histcnt / $totdir`
                dircnt=`expr $histcnt % $totdir`
                let "index = $moncnt * $totdir + $dircnt"
                echo "Reading histogram for ${months[moncnt]} for ${direct[dircnt]} degrees"
                # check for number of values in each record
                bincnt=0
                for value in ${record}
                do
                    let "bincnt += 1"
                done
                # not all records have a full bin count
                while [ $bincnt -lt 25 ]
                do
                    # pad record out to full bin count
                    # all missing values are 1.0
                    record="${record} 1.0"
                    let "bincnt += 1" 
                done
                echo "${record}" >> ${thistfile}
                let "histcnt += 1" 
            fi
        fi
        if [[ ${linecnt} -ge 212 ]]
        then
            if [[ ${linecnt} -le 213 ]]
            then
                # extract max min ratio and hour of max
                echo "${record}" >> ${tmaxminfile}
            fi
        fi

    done < ${filename[nsta]}

    # create transposed block file with histogram values down the columns
    # transpose has 100 line limit, workaround
    # abut (and maybe gnuplot has line length limit, workaround
    let "dend = totdir - 1"
    for ndir in $( seq 0 ${dend} )
    do
        # create double blank line to seperate indexes for gnuplot for separate stations
        echo "" >> ${ndir}${histfile}
        echo "" >> ${ndir}${histfile}

        echo "Transposing direction ${direct[$ndir]} degrees"
        # first line has x values
        linex 1 < ${thistfile} | transpose > temp
        let "mend = totmon - 1"
        getlines=""
        for nmon in $( seq 0 ${mend} )
        do
            let "getline = 2 + $ndir + $nmon * $totdir"
            getlines=`echo "${getlines} ${getline}"`
        done
        linex ${getlines} < ${thistfile} | transpose > temp1
        abut temp temp1 >> ${ndir}${histfile}
    done

    # append to max min hour file with columns
    cat ${tmaxminfile} | transpose >> ${maxminfile}

done

# set up pane sizes for a 4 across, 3 down
tacross=4
tdown=3
wsize=$( echo "scale=4; 0.999/${tacross}" | bc -l )
hsize=$( echo "scale=4; 0.999/${tdown}" | bc -l )

# set up pane origins
let "eacross = $tacross - 1"
let "edown = $tdown - 1"
for nh in $( seq 0 $edown )
do
    for nw in $( seq 0 $eacross )
    do
        let "nmon = $nw * $tdown + $nh"
        worigin[nmon]=$( echo "scale=4; ${nw}*${wsize}" | bc -l )
        horigin[nmon]=$( echo "scale=4; (${edown}-${nh})*${hsize}" | bc -l )
    done
done

# create gnuplot scripts for wind direction probabilities
echo "set terminal pdf font 'Helvetica,3'" >> ${plotfile}
echo "set output '${outfile}'" >> ${plotfile}
echo "set size 1.0, 1.0" >> ${plotfile}
echo "set origin 0.0, 0.0" >> ${plotfile}

for nsta in $( seq 1 ${totsta} )
do
    let "staidx = nsta - 1"
    echo "set multiplot" >> ${plotfile}
    echo "set label 1 '${staname[nsta]}' at screen 0.5,0.98 center front font 'Helvetica,4'" >> ${plotfile}
    index=0
    let "end = totmon - 1"
    for nmon in $( seq 0 ${end} )
    do
        let "col = 2 + ${index}"
        # w-width h-height
        echo "set size ${wsize}, ${hsize}" >> ${plotfile}
        echo "set origin ${worigin[nmon]}, ${horigin[nmon]}" >> ${plotfile}
        echo "set label 3 '${months[nmon]}' at graph 0.05,0.95 center front font 'Helvetica,4'" >> ${plotfile}
        #echo "set title '${months[nmon]}, ${stanum[nsta]}'" >> ${plotfile}
        # plot the wind rose
        echo "set size square" >> ${plotfile}
        echo "unset key" >> ${plotfile}
        echo "unset xlabel" >> ${plotfile}
        echo "unset ylabel" >> ${plotfile}
        echo "set angles degrees" >> ${plotfile}
        echo "set polar" >> ${plotfile}
        echo "set grid polar 22.5" >> ${plotfile}
        echo "unset border" >> ${plotfile}
        echo "set xtics 0,5,${probmax} axis scale 0.5 font 'Helvetica,3'" >> ${plotfile}
        echo "set ytics 0,${probmax},${probmax} axis scale 0.5 font 'Helvetica,3'" >> ${plotfile}
        echo "set size ratio 1" >> ${plotfile}
        echo "set xrange [-${probmax}:${probmax}]" >> ${plotfile}
        echo "set yrange [-${probmax}:${probmax}]" >> ${plotfile}
        echo "plot '${dirfile}' index ${staidx} u (90-\$1):(\$${col}) w i lw 5" >> ${plotfile}
        # plot calm on same pane
        echo "unset polar" >> ${plotfile}
        let "col = $nmon + 2"
        echo "plot '${calmfile}' index ${staidx} u 1:${col} w i lw 5" >> ${plotfile}
        let "index += 1"
    done
    # this command necessary to display all plots on some terminals
    echo "unset multiplot" >> ${plotfile}
done

# get ready for next graph
echo "reset" >> ${plotfile}

# create gnuplot scripts for histograms
let "enddir = totdir - 1"
for ndir in $( seq 0 ${enddir} )
do
    echo "set multiplot" >> ${plotfile}
    echo "set label 1 '${direct[ndir]} degrees' at screen 0.5,0.99 right front font 'Helvetica,4'" >> ${plotfile}
    echo "set label 2 ' CDF wind (m/s)' at screen 0.5,0.99 left front font 'Helvetica,4'" >> ${plotfile}
    let "endmon = totmon - 1" 
    for nmon in $( seq 0 ${endmon} )
    do
        let "col = 2 + ${nmon}"
        # w-width h-height
        echo "set size ${wsize}, ${hsize}" >> ${plotfile}
        echo "set origin ${worigin[nmon]}, ${horigin[nmon]}" >> ${plotfile}
        echo "set label 3 '${months[nmon]}' at graph 0.01,0.95 left front font 'Helvetica,4'" >> ${plotfile}

        #echo "set title '${months[nmon]}' font 'Helvetica,3'" >> ${plotfile}
        # plot probability scale
        if [ $nmon -eq 0 ]
        then
            echo "set key inside bottom right" >> ${plotfile}
        else
            echo "unset key" >> ${plotfile}
        fi
        echo "unset xlabel" >> ${plotfile}
        echo "unset ylabel" >> ${plotfile}
        echo "unset border" >> ${plotfile}
        # only put labels on left and bottom axes
        if [ $nmon -eq 2 ] || [ $nmon -eq 5 ] || [ $nmon -eq 8 ] || [ $nmon -eq 11 ]
        then
            echo "set xtics ('1.5' 1.5,'8.5' 8.5,'20.5' 20.5,'45.5' 45.5) scale 0.5 font 'Helvetica,3'" >> ${plotfile}
        else
            echo "set xtics ('' 1.5,'' 8.5,'' 20.5,'' 45.5 ) scale 0.5 font 'Helvetica,1'" >> ${plotfile}
        fi
        if [ $nmon -eq 0 ] || [ $nmon -eq 1 ] || [ $nmon -eq 2 ]
        then
            echo "set ytics ('1' invnorm(0.01),'5' invnorm(0.05),\\" >> ${plotfile}
            echo "'20' invnorm(0.2),'40' invnorm(0.4),'60' invnorm(0.6),\\" >> ${plotfile}
            echo "'80' invnorm(0.8),'90' invnorm(0.9),\\" >> ${plotfile}
            echo "'95' invnorm(0.95),'99' invnorm(0.99),'99.9' invnorm(0.999)) scale 0.5 font 'Helvetica,3'" >> ${plotfile}
        else
            echo "set ytics ('' invnorm(0.01),'' invnorm(0.05),\\" >> ${plotfile}
            echo "'' invnorm(0.2),'' invnorm(0.4),'' invnorm(0.6),\\" >> ${plotfile}
            echo "'' invnorm(0.8),'' invnorm(0.9),\\" >> ${plotfile}
            echo "'' invnorm(0.95),'' invnorm(0.99),'' invnorm(0.999)) scale 0.5 font 'Helvetica,1" >> ${plotfile}
        fi
        #echo "set size ratio 1" >> ${plotfile}
        echo "set yrange [invnorm(0.01):invnorm(0.999)]" >> ${plotfile}
        echo "set xrange [1.5:45.5]" >> ${plotfile}
        echo "set logscale x" >> ${plotfile}
        echo "set grid" >> ${plotfile}
        echo "plot \\" >> ${plotfile}
        for nsta in $( seq 1 ${totsta} )
        do
            let "staidx = nsta - 1"
            if [ ${nsta} -eq ${totsta} ]
            then
                echo "'${ndir}${histfile}' index ${staidx} u 1:(invnorm(\$${col})) with lines title '${stanum[nsta]}'" >> ${plotfile}
            else
                echo "'${ndir}${histfile}' index ${staidx} u 1:(invnorm(\$${col})) with lines title '${stanum[nsta]}', \\" >> ${plotfile}
            fi
        done

    done
    # this command necessary to display all plots on some terminals
    echo "unset multiplot" >> ${plotfile}
done

# get ready for next graphs
echo "reset" >> ${plotfile}

echo "set ylabel 'Max/Min Wind Ratio'" >> ${plotfile}
boxwidth=$( echo "scale=4; 1/(${totsta}+1)" | bc -l )
echo "set boxwidth ${boxwidth} absolute" >> ${plotfile}
echo "set style fill solid 1.00 border -1" >> ${plotfile}
echo "set grid nopolar" >> ${plotfile}
echo "set grid noxtics nomxtics ytics nomytics noztics nomztics \\" >> ${plotfile}
echo "nox2tics nomx2tics noy2tics nomy2tics nocbtics nomcbtics" >> ${plotfile}
echo "set xrange [-0.5:11.5]" >> ${plotfile}
echo "set yrange [1:]" >> ${plotfile}
echo "plot \\" >> ${plotfile}
for nsta in $( seq 1 ${totsta} )
do
    let "staidx = nsta - 1"
    boxloc=$( echo "scale=4; -0.5+${nsta}*${boxwidth}" | bc -l )
    if [ ${nsta} -eq ${totsta} ]
    then
        echo "'${maxminfile}' index ${staidx} using (\$0+${boxloc}):2:xticlabels(1) with boxes title '${stanum[nsta]}'" >> ${plotfile}
    else
        echo "'${maxminfile}' index ${staidx} using (\$0+${boxloc}):2 with boxes title '${stanum[nsta]}', \\" >> ${plotfile}
    fi
done

echo "set ylabel 'Hour of Maximum Wind'" >> ${plotfile}
echo "set yrange [0:]" >> ${plotfile}
echo "plot \\" >> ${plotfile}
for nsta in $( seq 1 ${totsta} )
do
    let "staidx = nsta - 1"
    boxloc=$( echo "scale=4; -0.5+${nsta}*${boxwidth}" | bc -l )
    if [ ${nsta} -eq ${totsta} ]
    then
        echo "'${maxminfile}' index ${staidx} using (\$0+${boxloc}):3:xticlabels(1) with boxes title '${stanum[nsta]}'" >> ${plotfile}
    else
        echo "'${maxminfile}' index ${staidx} using (\$0+${boxloc}):3 with boxes title '${stanum[nsta]}', \\" >> ${plotfile}
    fi
done

gnuplot ${plotfile}