src/misc/t-avg.pl

   1 #!/usr/bin/perl -w
   2
   3 #   t-avg - Calculates statistics of a **series** of piv data.
   4 #
   5 #   Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008
   6 #   Gerber van der Graaf <gerber_graaf@users.sourceforge.net
   7
   8 #   This program is free software; you can redistribute it and/or modify
   9 #   it under the terms of the GNU General Public License as published by
  10 #   the Free Software Foundation; either version 2, or (at your option)
  11 #   any later version.
  12
  13 #   This program is distributed in the hope that it will be useful,
  14 #   but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 #   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 #   GNU General Public License for more details.
  17
  18 #   You should have received a copy of the GNU General Public License
  19 #   along with this program; if not, write to the Free Software Foundation,
  20 #   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  21 #------------------------------------------------------------------
  22
  23 $VERSION = q$Id: t-avg.pl,v 1.6 2006-03-04 12:37:08 gerber Exp $;
  24 $HELP = "Calculates time-avaraged mean and rms from a **series** of piv data";
  25 $USAGE = "gpiv_t-avg [-h|-help] [-n|-none] [-p|-print] [-relrms] [-uv] [-v|-version]
  26 [-dn|dir_name] [-fd|-first_dir N] [-ld|-last_dir N] [-dx|dir_prefix]
  27 [-fb|-file_basename] [-fext|-file_extension] [-ff|-first_file N]
  28 [-fl|-last_file N] [-fx|file_prefix]
  29
  30
  31 keys:
  32 -h:      on-line help
  33 -n:      suppresses real execution
  34 -p:      prints process parameters/variables to stdout
  35 -relrms: prints relative rms, related to the absolute mean value
  36 -uv:     writes u and v values to separate outputs
  37 -v:      prints version
  38 -dn:     directory name to be processed
  39 -df:     starting at first directory number N
  40 -dl:     ending at last directory number N
  41 -dx:     prefix numbering to directory base name
  42 -fb:     file base name to be processed
  43 -fext:   add an extension after the file basename + number (without
  44          leading \".\")
  45 -ff:     starting at first file number N
  46 -fl:     ending at last file number N
  47 -fx:     prefix numbering to file base name
  48 ";
  49
  50 #----------------- Command arguments handling ----------
  51 $opt_h = 0;
  52 $opt_n = 0;
  53 $opt_p = 0;
  54 $opt_uv = 0;
  55 $opt_relrms = 0;
  56 $opt_v = 0;
  57 $opt_dx = 0;
  58 $opt_fx = 0;
  59
  60 use Getopt::Long;
  61 GetOptions('h|help', 'n|none', 'p|print', 'v|version',
  62 'relrms',
  63 'uv',
  64 'dx|dir_prefix', 'fx|file_prefix',
  65 'fb|file_basename=s' => \$file_base_name,
  66 'fext|file_extension=s' => \$file_ext,
  67 'ff|first_file=i' => \$file_first_nr,
  68 'fl|last_file=i' => \$file_last_nr,
  69 'dn|dir_name=s' => \$dir_base_name,
  70 'df|dir_first=i' => \$dir_first_nr,
  71 'dl|dir_last=i' => \$dir_last_nr);
  72
  73
  74
  75 if ($opt_h) {
  76   print ("$HELP\n");
  77   print ("$USAGE\n");
  78   exit;
  79 }
  80 if ($opt_v) {
  81   print ("$VERSION\n");
  82   exit;
  83 }
  84
  85 if ($#ARGV != -1) {
  86   print ("Usage: $USAGE\n");
  87   exit;
  88 }
  89
  90 #---------------------------------------------------- Initializing variables
  91 $first = 1;
  92 $k = 0;
  93
  94 if (!$file_base_name) {$file_base_name = "";}
  95 if (!$file_first_nr) {$file_first_nr = 0;}   # First data set to be scanned
  96 if (!$file_last_nr) {$file_last_nr = 2; }    # Last data set to be scanned
  97
  98
  99 if (!$dir_base_name) {
 100   $dir_base_name_defined = 0; $dir_base_name = "./";
 101 } else {
 102   $dir_base_name_defined = 1;
 103 }
 104 if (!$dir_first_nr) {$dir_first_nr = 0;}
 105 if (!$dir_last_nr) {$dir_last_nr = 0;}
 106
 107
 108 if ($opt_uv) {
 109     if ($file_ext) {
 110         $outfile_A = $dir_base_name.$file_base_name.".ta_U".".$file_ext".".piv";
 111         $outfile_B = $dir_base_name.$file_base_name.".ta_V".".$file_ext".".piv";
 112     } else {
 113         $outfile_A = $dir_base_name.$file_base_name.".ta_U.piv";
 114         $outfile_B = $dir_base_name.$file_base_name.".ta_V.piv";
 115     }
 116 } else {
 117     if ($file_ext) {
 118         $outfile_A = $dir_base_name.$file_base_name.".ta_mean".".$file_ext".".piv";
 119         $outfile_B = $dir_base_name.$file_base_name.".ta_sdev".".$file_ext".".piv";
 120     } else {
 121         $outfile_A = $dir_base_name.$file_base_name.".ta_mean.piv";
 122         $outfile_B = $dir_base_name.$file_base_name.".ta_sdev.piv";
 123     }
 124 }
 125
 126 if ($opt_n || ($file_base_name && $opt_p)) {
 127   print ("\nfile_base_name=$file_base_name file_first_nr=$file_first_nr file_last_nr=$file_last_nr\n");
 128 }
 129
 130 if ($opt_n || ($dir_base_name && $opt_p)) {
 131   print ("\ndir_base_name=$dir_base_name dir_nr_first=$dir_first_nr dir_nr_last=$dir_last_nr\n");
 132 }
 133
 134 #--------------------------------------------- Print filenames
 135 if ($opt_p || $opt_n) {
 136   printf ("Output: %s\n", $file_name);
 137   if ($opt_uv) {
 138     printf ("Data U-file: %s\n", $outfile_A);
 139     printf ("Data V-file: %s\n", $outfile_B);
 140   } else {
 141     printf ("Data mean file: %s\n", $outfile_A);
 142     printf ("Data svar file: %s\n", $outfile_B);
 143   }
 144 }
 145
 146 #---------------------------------------- loops for directory and data sets:
 147 for ($dir_nr = $dir_first_nr; $dir_nr <= $dir_last_nr; $dir_nr++) {
 148
 149   if ($dir_base_name_defined) {
 150     if (!$opt_dx) {
 151       $target_dir = $dir_base_name.$dir_nr."/";
 152     } else {
 153       $target_dir = $dir_nr.$dir_base_name."/";
 154     }
 155   } else {
 156     $target_dir = $dir_base_name;
 157   }
 158
 159   if ($opt_p || $opt_n) {
 160     printf ("\ntarget_dir[%d] = %s\n", $dir_nr, $target_dir);
 161   }
 162
 163   for ($file_nr = $file_first_nr; $file_nr <= $file_last_nr; $file_nr++) {
 164
 165     if (!$opt_fx) {
 166 #print ("append\n");
 167       $file_name=$target_dir.$file_base_name.$file_nr;
 168       if ($file_ext) {$file_name=$file_name.".$file_ext"};
 169     } else {
 170 #print ("prefix\n");
 171       $file_name=$target_dir.$file_nr.$file_base_name;
 172       if ($file_ext) {$file_name=$file_name.".$file_ext"};
 173     }
 174
 175     $file_name_piv = $file_name.".piv";
 176     $infile = $file_name_piv;
 177
 178
 179 #--------------------------------------------- Print input filename
 180     if ($opt_p || $opt_n) {
 181       printf ("Data piv file: %s\n", $file_name_piv);
 182     }
 183
 184     if (!$opt_n) {
 185     open (IN,"$infile") || die 'can\'t open $infile';
 186
 187
 188 #----------------------------------------- Starts reading the lines with data
 189     while (<IN>) {
 190       chomp;
 191 #---------------------------------- skip blank lines or ; sign at first col
 192       if ($_ ne "" &&  (!/^(\#|;)/)) {
 193 #------------------------------------ substitutes eventually komma's by dot
 194 #------------------------------ for example at LaVision's PIV program DaVis
 195         s/,/./g;
 196 #---------------------------------------- removes eventually leading spaces
 197         s/^ *//;
 198 #------------------- splitting line $_ up in words; #same as @words = split
 199           @words = split(/\s+/,$_);
 200
 201
 202 #------------ As the total number of data is unknown in forward, the variables
 203 # are set to zero when the first data set is scanned. Provided all data sets
 204 # will have the same number (and positions!!)
 205         if ($first) {
 206           $total_dx[$k] = 0;
 207           $total_dy[$k] = 0;
 208           $total_dx_2[$k] = 0;
 209           $total_dy_2[$k] = 0;
 210           $count[$k] = 0;
 211         }
 212
 213         $x_pos[$k] = $words[0];
 214         $y_pos[$k] = $words[1];
 215         $dx[$k] = $words[2];
 216         $dy[$k] = $words[3];
 217         $snr[$k] = $words[4];
 218         $peak_nr[$k] = $words[5];
 219
 220 # to avoid warning message $snr is used twice
 221         $snr[$k] = 0.0;
 222
 223 #------------------------------------------------ Take sum and sum of squares
 224         if ($peak_nr[$k] != -1) {
 225           $total_dx[$k] += $dx[$k];
 226           $total_dy[$k] += $dy[$k];
 227           $total_dx_2[$k] += $dx[$k]**2;
 228           $total_dy_2[$k] += $dy[$k]**2;
 229           $count[$k]++;
 230         }
 231         $k++;
 232       }
 233     }
 234
 235     if ($first) {
 236       $ndata = $k;
 237       if ($opt_p ) {
 238         printf ("\ndata = %d\n", $ndata);
 239       }
 240       $first = 0;
 241     }
 242     $k = 0;
 243     close (IN) || die 'can\'t close $infile';
 244   }
 245 }
 246 }
 247
 248 #-------------------------------------- Calculating mean and rms or variances
 249 if (!$opt_n) {
 250 for ($k=0; $k < $ndata; $k++) {
 251   if ($count[$k] > 1) {
 252     $mean_x[$k]= $total_dx[$k]/$count[$k];
 253     $svar_x[$k]= sqrt(($total_dx_2[$k]/($count[$k]-1) -
 254                     $total_dx[$k]**2 / ($count[$k]*($count[$k]-1))));
 255
 256     $mean_y[$k]= $total_dy[$k]/$count[$k];
 257     $svar_y[$k]= sqrt(($total_dy_2[$k]/($count[$k]-1) -
 258                     $total_dy[$k]**2 / ($count[$k]*($count[$k]-1))));
 259
 260   } else {
 261     $mean_x[$k] = $k;
 262     $svar_x[$k] = 0.;
 263
 264     $mean_y[$k] = $k;
 265     $svar_y[$k] = 0.;
 266     $count[$k] = -1;
 267   }
 268   $mean_xy[$k] = sqrt(($mean_x[$k]**2 + $mean_y[$k]**2));
 269   $svar_xy[$k] = sqrt(($svar_x[$k]**2 + $svar_y[$k]**2));
 270
 271 }
 272
 273
 274
 275 #----------------------------------------------------------- Writing results
 276 open (OUT_A,">$outfile_A") || die 'can\'t open $outfile_A';
 277 open (OUT_B,">$outfile_B") || die 'can\'t open $outfile_B';
 278 if ($opt_uv) {
 279     printf(OUT_A "#$VERSION");
 280     printf(OUT_A "\n#   x           y     x_mean    x_sdev    count");
 281     printf(OUT_B "#$VERSION");
 282     printf(OUT_B "\n#   x           y     y_mean    y_sdev    count");
 283
 284     for ($k=0; $k < $ndata; $k++) {
 285       printf(OUT_A "\n%-f %-f %-f %-f %-f  %-2d",
 286              $x_pos[$k], $y_pos[$k], $mean_x[$k], $svar_x[$k], $count[$k]);
 287
 288       printf(OUT_B "\n%-f %-f %-f %-f %-f  %-2d",
 289              $x_pos[$k], $y_pos[$k], $mean_y[$k], $svar_y[$k], $count[$k]);
 290     }
 291
 292 } elsif ($opt_relrms) {
 293     printf(OUT_A "#$VERSION");
 294     printf(OUT_A "\n#   x           y     x_mean    y_mean sdev/mean count");
 295     printf(OUT_B "#$VERSION");
 296     printf(OUT_B "\n#   x           y     x_sdev    y_sdev sdev sdev/mean count");
 297
 298     for ($k=0; $k < $ndata; $k++) {
 299       printf(OUT_A "\n%f %f %f %f %f %d",
 300              $x_pos[$k], $y_pos[$k], $mean_x[$k], $mean_y[$k],
 301 #            $mean_xy[$k],
 302 $svar_xy[$k] / $mean_xy[$k],
 303 $count[$k]);
 304
 305       printf(OUT_B "\n%f %f %f %f %f %f %d",
 306              $x_pos[$k], $y_pos[$k], $svar_x[$k], $svar_y[$k],
 307              $svar_xy[$k], $svar_xy[$k] / $mean_xy[$k],  $count[$k]);
 308     }
 309
 310 } else {
 311     printf(OUT_A "#$VERSION");
 312     printf(OUT_A "\n#   x           y     x_mean    y_mean    sdev  count");
 313     printf(OUT_B "#$VERSION");
 314     printf(OUT_B "\n#   x           y     x_sdev    y_sdev    sdev  count");
 315
 316     for ($k=0; $k < $ndata; $k++) {
 317       printf(OUT_A "\n%-f %-f %-f %-f %-f  %-2d",
 318              $x_pos[$k], $y_pos[$k], $mean_x[$k], $mean_y[$k],
 319              $svar_xy[$k], $count[$k]);
 320
 321       printf(OUT_B "\n%-f %-f %-f %-f %-f  %-2d",
 322              $x_pos[$k], $y_pos[$k], $svar_x[$k], $svar_y[$k],
 323              $svar_xy[$k], $count[$k]);
 324     }
 325 }
 326
 327 close (OUT_A) || die 'can\'t close $outfile';
 328 close (OUT_B) || die 'can\'t close $outfile';
 329 }
 330
 331 #
 332 # Thats all folks.
 333 #