#!/bin/bash

export LANG=en_US

print_help() {
    echo "anasum version 1.0" && echo && \
    echo "options: -i <input>  set input path(s) to <input>" && \
    echo "         -o <output> set output path to <output>" && \
    echo "         -l <depth>  process <depth> subdirectories" && \
    echo "         -p          replace old output path w/o prompting" && \
    echo "         -h          display this help and exit" && echo && \
    echo "all path names must be given with a trailing slash" && echo
}

DEPTH=1
NOPROMPT=FALSE
while getopts :i:o:l:ph OPT
do
  case $OPT in
  i) INPUT=$INPUT" "$OPTARG ;; 
  o) OUTPUT=$OPTARG ;;
  p) NOPROMPT=TRUE ;;
  l) DEPTH=$OPTARG ;;
  h) print_help && exit 0 ;;
  \?)
    shift `expr $OPTIND - 1`
    if [ "$1" = "--clean" ]; then print_help && exit 0
    else 
      echo -n "anasum: error: unrecognized option "
      if [ $OPTARG != "-" ]; then echo "'-$OPTARG'. try '-h'"
      else echo "'$1'. try '-h'"
      fi
      print_help && exit 1
    fi
    shift 1
    OPTIND=1
  esac
done
if [ "$INPUT" = "" ]; then
  echo "anasum: no input paths"
  exit 1
fi
if [ "$OUTPUT" = "" ]; then
  echo "anasum: no output paths"
  exit 1
fi
if test -d $OUTPUT; then
  if [ "$NOPROMPT" = "TRUE" ]; then
    rm -r $OUTPUT
  else
    if ! read -t30 -p "anasum: overwrite directory '$OUTPUT' (y/n) ? " ANS; then
      echo -e "\nanasum: timeout"
      exit 1
    else 
      if [ "$ANS" != "y" ]; then
        echo "anasum: invalid output path"
        exit 1
      fi
      rm -r $OUTPUT
    fi
  fi
fi
mkdir $OUTPUT
SDIR=`echo "$INPUT" | cut -d' ' -f2`
echo -n "scanning directory '$SDIR' ... "
if ! cd $SDIR; then exit 1; fi
FILES=`find . -maxdepth $DEPTH -name \*.dat`; NFILES=0
for I in $FILES; do (( ++NFILES )); done
cd $OLDPWD
echo "$NFILES files"

for FILE in $FILES; do

  echo -n "file '$FILE': "

  CFS=""
  for DIR in $INPUT; do
    if ! test -f $DIR$FILE; then
      echo "anasum: file '$FILE' not found in '$DIR'. skipped."
    else
      CFS=$CFS" "$DIR$FILE
    fi
  done
  if [ "$CFS" = "" ]; then continue; fi

  awk 'BEGIN{ ofile="'$OUTPUT$FILE'"; init=0; write=0; }{
    if (init==0) {
      tp=$1; md=tp%10; uf=0; of=0; euf=0; eof=0; n=0; nb=$2-2;
      for (i=0;i<=nb;++i) { x[i]=$3+i*($4-$3)/nb; y[i]=0; y2[i]=0; ps[i]=0; }
      print "x-range ["x[0]","x[nb]"], "$2-2" bins, type "tp;
      printf "processing files ";
      write=1;
    }
    if (FNR==1) {
      if (FILENAME==ofile) nextfile;
      printf ".";
      if (md==0) cn=$7; else cn=$9;
      if (cn<0) cn+=2^32;
      if (md==0) { uf+=cn*$5; of+=cn*$6; }
      else { uf+=cn*$5; of+=cn*$7; 
        euf+=cn*((cn-1)*$6+$5*$5); eof+=cn*((cn-1)*$8+$7*$7); }
      l=0;
      n+=cn;
      ++init;
      next;
    }
    dev=x[l]-$1;
    if (x[l]>1.0e-12 || x[l]<-1.0e-12) dev=dev/x[l];
    if (dev>1.0e-3 || dev<-1.0e-3) {
      print "\nerror: bins "l" ("x[l]","$1") " \
        " do not coincide in "FILENAME" (rel. dev. "dev"). abort.";
      write=0;
      exit;
    }
    y[l]+=cn*$2;
    if (md>0) y2[l]+=cn*((cn-1)*$3+$2*$2);
    if (md>1) ps[l]+=cn*$4;
    ++l;
  }END{ 
    if (write>0) {
      dls=split(ofile,dirs,"/");
      cdir=dirs[1];
      for (i=2;i<dls;++i) {
        cdir=cdir"/"dirs[i];
        system("if ! test -d "cdir"; then mkdir "cdir"; fi");
      }
      print " done\ncombined "init" files into "ofile;
      if (md==0) {
        print tp" "nb+2" "x[0]" "x[nb]" "uf/n" "of/n" "n > ofile;
        for (i=0;i<=nb;++i) { print x[i]" "y[i]/n > ofile; } 
      }
      else {
        euf=sqrt((euf/n-uf*uf/(n*n))/(n-1));
        eof=sqrt((eof/n-of*of/(n*n))/(n-1));
        print tp" "nb+2" "x[0]" "x[nb]" "uf/n" "euf" "of/n" "eof" "n > ofile;
        for (i=0;i<=nb;++i) { 
          y2[i]=(y2[i]/n-y[i]*y[i]/(n*n))/(n-1);
          if (md==1) print x[i]" "y[i]/n" "sqrt(y2[i]) > ofile; 
          else print x[i]" "y[i]/n" "sqrt(y2[i])" "ps[i]/n > ofile; 
        } 
      }
    }
  }' $CFS

done