#! /bin/bash

# invoke gnuplot to plot comparison between two plot.out files
# on days that erosion occurs in either of the files.
# invocation: plot_plot file1 file2

# set up columns titles for plot.out file
mincol=6
maxcol=84
coltitle[6]=`echo "6) Total erosion loss (kg/m^2)"`
coltitle[7]=`echo "7) Suspension erosion loss (kg/m^2)"`
coltitle[8]=`echo "8) pm10 erosion loss (kg/m^2)"`
coltitle[9]=`echo "9) pm2.5 erosion loss (kg/m^2)"`
coltitle[10]=`echo "10) Maximum wind speed (m/s)"`
coltitle[11]=`echo "11) Wind direction (degrees from N)"`
coltitle[12]=`echo "12) Precipitation (mm)"`
coltitle[13]=`echo "13) Snow depth (mm)"`
coltitle[14]=`echo "14) Soil surface water content (kg h20/kg soil)"`
coltitle[15]=`echo "15) Crust fraction"`
coltitle[16]=`echo "16) Crust thickness (mm)"`
coltitle[17]=`echo "17) Loose mass on crust (Mg/m^2)"`
coltitle[18]=`echo "18) Loose mass crust cover fraction"`
coltitle[19]=`echo "19) Ridge height (mm)"`
coltitle[20]=`echo "20) Ridge width (mm)"`
coltitle[21]=`echo "21) Ridge spacing (mm)"`
coltitle[22]=`echo "22) Ridge orientation (degrees from N)"`
coltitle[23]=`echo "23) Ridge spacing in the wind direction (mm)"`
coltitle[24]=`echo "24) Dike height (mm)"`
coltitle[25]=`echo "25) Dike spacing (mm)"`
coltitle[26]=`echo "26) Random roughness (mm)"`
coltitle[27]=`echo "27) Aggregate size geometric mean diameter (mm)"`
coltitle[28]=`echo "28) Aggregate size geometric standard deviation (mm)"`
coltitle[29]=`echo "29) Aggregate size minimum size limit (mm)"`
coltitle[30]=`echo "30) Aggregate size maximum size limit (mm)"`
coltitle[31]=`echo "31) Aggregate stability (J/cm^2)"`
coltitle[32]=`echo "32) Soil bulk density (Mg/m^3)"`
coltitle[33]=`echo "33) Total Flat cover fraction"`
coltitle[34]=`echo "34) Total mass weighted average height (m)"`
coltitle[35]=`echo "35) Total stem area index"`
coltitle[36]=`echo "36) Total leaf area index"`
coltitle[37]=`echo "37) Crop height (m)"`
coltitle[38]=`echo "38) Crop Stem area index"`
coltitle[39]=`echo "39) Crop Leaf area index"`
coltitle[40]=`echo "40) Crop biodrag"`
coltitle[41]=`echo "41) Residue mass weighted average height (m)"`
coltitle[42]=`echo "42) Residue stem area index"`
coltitle[43]=`echo "43) Residue leaf area index"`
coltitle[44]=`echo "44) Residue biodrag"`
coltitle[45]=`echo "45) Aggregate coefficient of abrasion (1/m)"`
coltitle[46]=`echo "46) Crust coefficient of abrasion (1/m)"`
coltitle[47]=`echo "47) PM10 abrasion fraction"`
coltitle[48]=`echo "48) PM10 emission fraction"`
coltitle[49]=`echo "49) PM10 breakage fraction"`
coltitle[50]=`echo "50) PM2.5 abrasion fraction"`
coltitle[51]=`echo "51) PM2.5 emission fraction"`
coltitle[52]=`echo "52) PM2.5 breakage fraction"`
coltitle[53]=`echo "53) Surface soil fraction less than 10 microns initial condition"`
coltitle[54]=`echo "54) Surface soil fraction less than 100 microns initial condition"`
coltitle[55]=`echo "55) Surface soil fraction less than 0.84 mm initial condition"`
coltitle[56]=`echo "56) Surface soil fraction less than 2.0 mm initial condition"`
coltitle[57]=`echo "57) Surface soil fraction less than 10 microns"`
coltitle[58]=`echo "58) Surface soil fraction less than 100 microns"`
coltitle[59]=`echo "59) Surface soil fraction less than 0.84 mm"`
coltitle[60]=`echo "60) Surface soil fraction less than 2.0 mm"`
coltitle[61]=`echo "61) 1 - erosion was entered"`
coltitle[62]=`echo "62) 1 - snow depth prevented erosion"`
coltitle[63]=`echo "63) anemometer located friction velocity"`
coltitle[64]=`echo "64) site surface random roughness adjusted friction velocity"`
coltitle[65]=`echo "65) site surface oriented roughness adjusted friction velocity"`
coltitle[66]=`echo "66) site biodrag adjusted friction velocity"`
coltitle[67]=`echo "67) site friction velocity"`
coltitle[68]=`echo "68) bare friction velocity threshold"`
coltitle[69]=`echo "69) flat cover increases threshold"`
coltitle[70]=`echo "70) surface wetness increases threshold"`
coltitle[71]=`echo "71) ag density increases threshold"`
coltitle[72]=`echo "72) resultant threshold friction velocity"`
coltitle[73]=`echo "73) ratio anemometer located friction velocity"`
coltitle[74]=`echo "74) ratio site surface random roughness adjusted friction velocity"`
coltitle[75]=`echo "75) ratio site surface oriented roughness adjusted friction velocity"`
coltitle[76]=`echo "76) ratio site biodrag adjusted friction velocity"`
coltitle[77]=`echo "77) ratio bare friction veolocity threshold"`
coltitle[78]=`echo "78) ratio flat cover increases threshold"`
coltitle[79]=`echo "79) ratio surface wetness increases threshold"`
coltitle[80]=`echo "80) ratio ag density increases threshold"`
coltitle[81]=`echo "81) Surface soil fraction less than 0.84 mm"`
coltitle[82]=`echo "82) Surface soil fraction greater than 2 mm"`
coltitle[83]=`echo "83) Surface soil aerodynamic roughness (mm)"`
coltitle[84]=`echo "84) Soil surface fraction which does not emit"`

# check for presence of command arguments
if [ ! -n "$1" ]
then
    echo 'invoke gnuplot to plot comparison between two plot.out files'
    echo 'on days that erosion occurs in either of the files.'
    echo 'invocation: plot_plot file1 file2'
    exit
fi

# select lines that have a erosion in either file
dm "if x$mincol < 0.0 then INLINE else SKIP" < $1 > lines1.out
dm "if x$mincol < 0.0 then INLINE else SKIP" < $2 > lines2.out

# get list of line numbers to be extracted
cat lines1.out lines2.out | sort -g | uniq > lines.out
#selected_lines=`cat lines.out`
linecnt=$(wc -w < lines.out)
if [ "${linecnt}" -gt "1000" ]
then
    selected_lines=$(cat lines.out | linex 1-1000 )
else
    selected_lines=$(cat lines.out)
fi
echo $(wc -w < lines.out) erosion events > plot.rep
echo $selected_lines

# extract same lines from both files and remove blank lines
cat $1 | linex $selected_lines | dm "if N > 0 then INPUT else SKIP" > plot1.out
cat $2 | linex $selected_lines | dm "if N > 0 then INPUT else SKIP" > plot2.out

# indicate daysim and day of year for extracted lines
echo "daysim  day-of-year year erosion1 swc1 gmd1 gsd1 erosion2 swc2 gmd2 gsd2" >> plot.rep
dm x1 x2 x5 x6 x12 x16 x17 < plot1.out > rep1.out
dm x6 x12 x16 x17 < plot2.out > rep2.out
abut rep1.out rep2.out >> plot.rep
cat plot.rep

# limit output to days when erosion differences are greatest
plotcnt=100
daycount=$(dm "abs(x4-x8)" < plot.rep | stats n)
if [ "${daycount}" -gt "${plotcnt}" ]
then
    # select ten highest values
    limitvalue=$(dm "abs(x4-x8)" < plot.rep | sort -r | linex ${plotcnt} )
    dm "if abs(x4-x8) >= ${limitvalue} then x1 else SKIP" < plot.rep > daysim.sel
    # redo line extraction
    rm -f plot1.sel plot2.sel
    while read record
    do
        dm "if x1=${record} then INPUT else SKIP" < plot1.out >> plot1.sel
        dm "if x1=${record} then INPUT else SKIP" < plot2.out >> plot2.sel
    done < daysim.sel
    # move selected values into plot files used below
    mv plot1.sel plot1.out
    mv plot2.sel plot2.out
fi

# plot data file columns in comparison plots
plotfile=plot.tmp

echo reset > $plotfile

echo "set terminal pdf" >> $plotfile
echo "set output 'plot.pdf'" >> $plotfile

# label the x and y axes with input file name
echo "set xlabel '$1'" >> $plotfile
echo "set ylabel '$2'" >> $plotfile

# extract daysim for point labels
colex 1 < plot1.out > column3.out
for colnum in $(seq $mincol $maxcol)
do

  # extract the numbered column from the two files and place them in the plot file
  colex $colnum < plot1.out > column1.out
  colex $colnum < plot2.out > column2.out
  abut column1.out column2.out column3.out > ${colnum}.out
  
  # create lines in plot file to plot ${colnum}.out
  echo "set title '${coltitle[$colnum]}'" >> $plotfile
  echo "plot x notitle lt 3, '${colnum}.out' using 1:2 with points notitle, '' using 1:2:3 with labels font 'Helvetica,2' tc lt 1 notitle" >> $plotfile
#  echo "pause -1" >> $plotfile

done

gnuplot $plotfile

# clean up all the files created

rm lines1.out
rm lines2.out
rm lines.out
rm plot1.out
rm plot2.out
rm rep1.out
rm rep2.out
rm column1.out
rm column2.out
rm column3.out

rm $plotfile
for colnum in $(seq $mincol $maxcol)
do
  rm ${colnum}.out
done