#! /bin/bash
# Last edited on 2009-12-12 00:19:09 by stolfi

pha="$1"; shift
logscale="$1"; shift
datafile="$1"; shift

echo "${0##*/} pha = ${pha} logscale = ${logscale} datafile = ${datafile}" 1>&2

seppref="/tmp/$$-sep"
epsfile="/tmp/$$.eps"

if [[ ${pha} -eq 0 ]]; then
  xrange="set xrange [2000.8:2010.8]"
elif [[ ${pha} -eq 1 ]]; then
  xrange="set xrange [2000.8:2004.8]"
elif [[ ${pha} -eq 2 ]]; then
  xrange="set xrange [2006.3:2010.7]"
else
  echo "invalid pha = [${pha}]" 1>&2 ; exit 1
fi

errmax="25000"; # Max absolute error.
ratmax="2.50";  # Max data/model ratio.

if [[ ${logscale} -gt 0 ]]; then 
  yfunc="(100*(column(i)+0.0)/column(j))";
  ymin="(100.0/ratmax/1.05)"
  ymax="(100.0*ratmax*1.05)"
  ymed="100.0"
  ytitle="mod/obs (%)"
  yscale="set logscale y; set ytics \\
    ( \\
      '%8.0f'  50 0, \\
      '%8.0f'  60 0, \\
      '%8.0f'  70 0, \\
      '%8.0f'  80 0, \\
      '%8.0f'  90 0, \\
      '%8.0f' 100 0, \\
      ''      110 1, \\
      '%8.0f' 120 0, \\
      ''      130 1, \\
      ''      140 1, \\
      '%8.0f' 150 0, \\
      ''      160 1, \\
      ''      170 1, \\
      ''      180 1, \\
      ''      190 1, \\
      '%8.0f' 200 0, \\
      ''      210 1, \\
      ''      220 1, \\
      ''      230 1, \\
      ''      240 1, \\
      '%8.0f' 250 0  \\
    )"
else
  yfunc="((column(i)+0.0) - column(j))";
  ymin="(-1.05*errmax)"
  ymax="(+1.05*errmax)"
  ymed="0.0"
  ytitle="mod-obs (art)"
  yscale="unset logscale y; set ytics 5000 format '%8.0f'; set mytics 5"
fi

# In this plot we assume that {su} is 0 or 1:
for su in 0 1; do
  cat ${datafile} \
    | gawk -v su=${su} \
        ' /^ *([\#]|$)/{ next; } 
          ($6 == su){ print; next; } 
          // { print ""; }
        ' \
    > ${seppref}-sz-${su}.txt
done

gnuplot <<EOF
set term postscript eps color solid "TimesRoman" 24
set output "${epsfile}"
set size 1.5, 0.75

errmax=${errmax}
ratmax=${ratmax}
yearsz(i)=(column(i)/365.25 + 2001)
yfunc(i,j)=${yfunc}

set samples 200

# set key reverse left Left samplen 1.0
set key samplen 1.0

set grid xtics
set xtics 1.0 format "        %4.0f"
set mxtics 12
${xrange}

${yscale}
set yrange[${ymin}:${ymax}]
set grid mytics ytics

plot \
  "${seppref}-sz-1.txt" using (yearsz(1)):(yfunc(7,5)) \
      title "${ytitle}" with linespoints lt 1 lw 1.0 pt 7 ps 0.75 lc rgb '#3366dd', \
  "${seppref}-sz-0.txt" using (yearsz(1)):(yfunc(7,5)) notitle with points pt 6 lc rgb '#6666aa', \
  (${ymed}) notitle with lines lt 1 lc '#778899'

quit
EOF

cat ${epsfile}
gv ${epsfile}
rm -f ${epsfile} ${seppref}-*.txt