#!/usr/bin/perl -w

#   t-avg - Calculates statistics of a **series** of piv data.
#
#   Copyright (C) 2002 Gerber van der Graaf

#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2, or (at your option)
#   any later version.

#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.

#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software Foundation,
#   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
#------------------------------------------------------------------

$VERSION = q$Id: t-avg.pl,v 1.6 2006/03/04 12:37:08 gerber Exp $;
$HELP = "Calculates time-avaraged mean and rms from a **series** of piv data";
$USAGE = "gpiv_t-avg [-h|-help] [-n|-none] [-p|-print] [-relrms] [-uv] [-v|-version] 
[-dn|dir_name] [-fd|-first_dir N] [-ld|-last_dir N] [-dx|dir_prefix]
[-fb|-file_basename] [-fext|-file_extension] [-ff|-first_file N] 
[-fl|-last_file N] [-fx|file_prefix]


keys:
-h:      on-line help
-n:      suppresses real execution
-p:      prints process parameters/variables to stdout
-relrms: prints relative rms, related to the absolute mean value
-uv:     writes u and v values to separate outputs
-v:      prints version
-dn:     directory name to be processed
-df:     starting at first directory number N
-dl:     ending at last directory number N
-dx:     prefix numbering to directory base name
-fb:     file base name to be processed
-fext:   add an extension after the file basename + number (without 
         leading \".\")
-ff:     starting at first file number N
-fl:     ending at last file number N
-fx:     prefix numbering to file base name
";

#----------------- Command arguments handling ----------
$opt_h = 0;
$opt_n = 0;
$opt_p = 0;
$opt_uv = 0;
$opt_relrms = 0;
$opt_v = 0;
$opt_dx = 0;
$opt_fx = 0;

use Getopt::Long;
GetOptions('h|help', 'n|none', 'p|print', 'v|version',
'relrms',
'uv',
'dx|dir_prefix', 'fx|file_prefix',
'fb|file_basename=s' => \$file_base_name, 
'fext|file_extension=s' => \$file_ext,
'ff|first_file=i' => \$file_first_nr,
'fl|last_file=i' => \$file_last_nr,
'dn|dir_name=s' => \$dir_base_name, 
'df|dir_first=i' => \$dir_first_nr,
'dl|dir_last=i' => \$dir_last_nr);



if ($opt_h) {
  print ("$HELP\n");
  print ("$USAGE\n");
  exit;
}
if ($opt_v) {
  print ("$VERSION\n");
  exit;
}

if ($#ARGV != -1) {
  print ("Usage: $USAGE\n");
  exit;
}

#---------------------------------------------------- Initializing variables
$first = 1;
$k = 0;

if (!$file_base_name) {$file_base_name = "";}
if (!$file_first_nr) {$file_first_nr = 0;}   # First data set to be scanned
if (!$file_last_nr) {$file_last_nr = 2; }    # Last data set to be scanned


if (!$dir_base_name) {
  $dir_base_name_defined = 0; $dir_base_name = "./"; 
} else {
  $dir_base_name_defined = 1;
}
if (!$dir_first_nr) {$dir_first_nr = 0;}
if (!$dir_last_nr) {$dir_last_nr = 0;}


if ($opt_uv) {
    if ($file_ext) {
        $outfile_A = $dir_base_name.$file_base_name.".ta_U".".$file_ext".".piv";
        $outfile_B = $dir_base_name.$file_base_name.".ta_V".".$file_ext".".piv";
    } else {
        $outfile_A = $dir_base_name.$file_base_name.".ta_U.piv";
        $outfile_B = $dir_base_name.$file_base_name.".ta_V.piv";
    }
} else {
    if ($file_ext) {
        $outfile_A = $dir_base_name.$file_base_name.".ta_mean".".$file_ext".".piv";
        $outfile_B = $dir_base_name.$file_base_name.".ta_sdev".".$file_ext".".piv";
    } else {
        $outfile_A = $dir_base_name.$file_base_name.".ta_mean.piv";
        $outfile_B = $dir_base_name.$file_base_name.".ta_sdev.piv";
    }
}

if ($opt_n || ($file_base_name && $opt_p)) {
  print ("\nfile_base_name=$file_base_name file_first_nr=$file_first_nr file_last_nr=$file_last_nr\n");
}

if ($opt_n || ($dir_base_name && $opt_p)) {
  print ("\ndir_base_name=$dir_base_name dir_nr_first=$dir_first_nr dir_nr_last=$dir_last_nr\n");
}

#--------------------------------------------- Print filenames
if ($opt_p || $opt_n) {  
  printf ("Output: %s\n", $file_name);
  if ($opt_uv) {
    printf ("Data U-file: %s\n", $outfile_A);
    printf ("Data V-file: %s\n", $outfile_B);
  } else {
    printf ("Data mean file: %s\n", $outfile_A);
    printf ("Data svar file: %s\n", $outfile_B);
  }
}

#---------------------------------------- loops for directory and data sets:
for ($dir_nr = $dir_first_nr; $dir_nr <= $dir_last_nr; $dir_nr++) {

  if ($dir_base_name_defined) {
    if (!$opt_dx) {
      $target_dir = $dir_base_name.$dir_nr."/";
    } else {
      $target_dir = $dir_nr.$dir_base_name."/";
    }
  } else { 
    $target_dir = $dir_base_name;
  }

  if ($opt_p || $opt_n) {
    printf ("\ntarget_dir[%d] = %s\n", $dir_nr, $target_dir);
  }

  for ($file_nr = $file_first_nr; $file_nr <= $file_last_nr; $file_nr++) {

    if (!$opt_fx) {
#print ("append\n");
      $file_name=$target_dir.$file_base_name.$file_nr;
      if ($file_ext) {$file_name=$file_name.".$file_ext"};
    } else {
#print ("prefix\n");
      $file_name=$target_dir.$file_nr.$file_base_name;
      if ($file_ext) {$file_name=$file_name.".$file_ext"};
    }

    $file_name_piv = $file_name.".piv";
    $infile = $file_name_piv;


#--------------------------------------------- Print input filename
    if ($opt_p || $opt_n) {  
      printf ("Data piv file: %s\n", $file_name_piv);
    }

    if (!$opt_n) { 
    open (IN,"$infile") || die 'can\'t open $infile';


#----------------------------------------- Starts reading the lines with data
    while (<IN>) {
      chomp;
#---------------------------------- skip blank lines or ; sign at first col
      if ($_ ne "" &&  (!/^(\#|;)/)) { 
#------------------------------------ substitutes eventually komma's by dot
#------------------------------ for example at LaVision's PIV program DaVis
	s/,/./g;
#---------------------------------------- removes eventually leading spaces
	s/^ *//;
#------------------- splitting line $_ up in words; #same as @words = split
	  @words = split(/\s+/,$_);


#------------ As the total number of data is unknown in forward, the variables
# are set to zero when the first data set is scanned. Provided all data sets 
# will have the same number (and positions!!) 
	if ($first) {
	  $total_dx[$k] = 0;
	  $total_dy[$k] = 0;
	  $total_dx_2[$k] = 0;
	  $total_dy_2[$k] = 0;
	  $count[$k] = 0;
	}

	$x_pos[$k] = $words[0];
	$y_pos[$k] = $words[1];
	$dx[$k] = $words[2];
	$dy[$k] = $words[3];
	$snr[$k] = $words[4];
	$peak_nr[$k] = $words[5];

# to avoid warning message $snr is used twice
	$snr[$k] = 0.0;

#------------------------------------------------ Take sum and sum of squares
	if ($peak_nr[$k] != -1) {
	  $total_dx[$k] += $dx[$k];
	  $total_dy[$k] += $dy[$k];
	  $total_dx_2[$k] += $dx[$k]**2;
	  $total_dy_2[$k] += $dy[$k]**2;
	  $count[$k]++;
	}
	$k++;
      }
    }

    if ($first) {
      $ndata = $k;
      if ($opt_p ) {
        printf ("\ndata = %d\n", $ndata);
      }
      $first = 0;
    }
    $k = 0;
    close (IN) || die 'can\'t close $infile';
  }
}
}

#-------------------------------------- Calculating mean and rms or variances
if (!$opt_n) {
for ($k=0; $k < $ndata; $k++) {
  if ($count[$k] > 1) {
    $mean_x[$k]= $total_dx[$k]/$count[$k];
    $svar_x[$k]= sqrt(($total_dx_2[$k]/($count[$k]-1) - 
		    $total_dx[$k]**2 / ($count[$k]*($count[$k]-1))));
		      
    $mean_y[$k]= $total_dy[$k]/$count[$k];
    $svar_y[$k]= sqrt(($total_dy_2[$k]/($count[$k]-1) - 
		    $total_dy[$k]**2 / ($count[$k]*($count[$k]-1))));
		      
  } else {
    $mean_x[$k] = $k;
    $svar_x[$k] = 0.;

    $mean_y[$k] = $k;
    $svar_y[$k] = 0.;
    $count[$k] = -1;
  }
  $mean_xy[$k] = sqrt(($mean_x[$k]**2 + $mean_y[$k]**2));
  $svar_xy[$k] = sqrt(($svar_x[$k]**2 + $svar_y[$k]**2));

}



#----------------------------------------------------------- Writing results
open (OUT_A,">$outfile_A") || die 'can\'t open $outfile_A';
open (OUT_B,">$outfile_B") || die 'can\'t open $outfile_B';
if ($opt_uv) {
    printf(OUT_A "#$VERSION");
    printf(OUT_A "\n#   x           y     x_mean    x_sdev    count");
    printf(OUT_B "#$VERSION");
    printf(OUT_B "\n#   x           y     y_mean    y_sdev    count");

    for ($k=0; $k < $ndata; $k++) {
      printf(OUT_A "\n%-f %-f %-f %-f %-f  %-2d",
	     $x_pos[$k], $y_pos[$k], $mean_x[$k], $svar_x[$k], $count[$k]);

      printf(OUT_B "\n%-f %-f %-f %-f %-f  %-2d",
	     $x_pos[$k], $y_pos[$k], $mean_y[$k], $svar_y[$k], $count[$k]);
    }

} elsif ($opt_relrms) {
    printf(OUT_A "#$VERSION");
    printf(OUT_A "\n#   x           y     x_mean    y_mean sdev/mean count");
    printf(OUT_B "#$VERSION");
    printf(OUT_B "\n#   x           y     x_sdev    y_sdev sdev sdev/mean count");

    for ($k=0; $k < $ndata; $k++) {
      printf(OUT_A "\n%f %f %f %f %f %d",
	     $x_pos[$k], $y_pos[$k], $mean_x[$k], $mean_y[$k],
#	     $mean_xy[$k], 
$svar_xy[$k] / $mean_xy[$k],
$count[$k]);

      printf(OUT_B "\n%f %f %f %f %f %f %d",
	     $x_pos[$k], $y_pos[$k], $svar_x[$k], $svar_y[$k],
	     $svar_xy[$k], $svar_xy[$k] / $mean_xy[$k],  $count[$k]);
    }

} else {
    printf(OUT_A "#$VERSION");
    printf(OUT_A "\n#   x           y     x_mean    y_mean    sdev  count");
    printf(OUT_B "#$VERSION");
    printf(OUT_B "\n#   x           y     x_sdev    y_sdev    sdev  count");

    for ($k=0; $k < $ndata; $k++) {
      printf(OUT_A "\n%-f %-f %-f %-f %-f  %-2d",
	     $x_pos[$k], $y_pos[$k], $mean_x[$k], $mean_y[$k], 
	     $svar_xy[$k], $count[$k]);

      printf(OUT_B "\n%-f %-f %-f %-f %-f  %-2d",
	     $x_pos[$k], $y_pos[$k], $svar_x[$k], $svar_y[$k], 
	     $svar_xy[$k], $count[$k]);
    }
}

close (OUT_A) || die 'can\'t close $outfile';
close (OUT_B) || die 'can\'t close $outfile';
}

#
# Thats all folks.
#
