src/misc/t-corr.c

   1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
   2
   3 /*-----------------------------------------------------------------------------
   4
   5    t-corr - calculates the velocity correlation as function of time
   6             (Eulerian correlation) from a series PIV data sets
   7
   8    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008
   9    Gerber van der Graaf <gerber_graaf@users.sourceforge.net
  10
  11    This program is free software; you can redistribute it and/or modify
  12    it under the terms of the GNU General Public License as published by
  13    the Free Software Foundation; either version 2, or (at your option)
  14    any later version.
  15
  16    This program is distributed in the hope that it will be useful,
  17    but WITHOUT ANY WARRANTY; without even the implied warranty of
  18    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19    GNU General Public License for more details.
  20
  21    You should have received a copy of the GNU General Public License
  22    along with this program; if not, write to the Free Software Foundation,
  23    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  24
  25 -----------------------------------------------------------------------------*/
  26
  27 #include <stdio.h>
  28 #include <stdlib.h>
  29 #include <assert.h>
  30 #include <math.h>
  31 #include <gpiv.h>
  32 #include "config.h"
  33 #include "git-rev.h"
  34
  35 /* #define PARFILE "scale.par"    */  /* Parameter file name */
  36 #define PARFILE "gpivrc"     /* Parameter file name */
  37 #define MAX_PIVSETS 1000     /* Maximum number of PIV data sets */
  38 #define MIN_SAMPLES 20       /* Minimum number of samples used for estimation */
  39 #define USAGE "\
  40 Usage: t-corr [-h | --help] [-p | --print] [-v | --version] \n\
  41 [-f | --file_first N]  [-l | --file_last N] [-x | --point_x F] \n\
  42 [-y | --point_y F] file_basename\n\
  43 \n\
  44 keys: \n\
  45 -h | --help:           this on-line help \n\
  46 -p | --print:          print parameters to stdout \n\
  47 -v | --version:        version number \n\
  48 -f | --file_first N    First file number of the series (default: 0) \n\
  49 -l | --file_last N     Last file number of the series (default: 0) \n\
  50 -x | --point_x F       x-position of data set to be analysed (default: 0.0) \n\
  51 -y | --point_y F       y-position of data set to be analysed (default: 0.0) \n\
  52 file_basename          File basename of a numbered series of PIV data sets \n\
  53 "
  54
  55 #define HELP  "\
  56 t-corr - calculates the Eulerian correlation from a series of PIV data sets"
  57
  58 gboolean print_par = FALSE;
  59 typedef struct _TimeCorrPar TimeCorrPar;
  60 struct _TimeCorrPar {
  61     gint first_file;
  62     gint last_file;
  63     gboolean first_file__set;
  64     gboolean last_file__set;
  65
  66     gfloat point_x;
  67     gfloat point_y;
  68     gboolean  point_x__set;
  69     gboolean  point_y__set;
  70 };
  71
  72
  73 void
  74 command_args(int argc,
  75              char *argv[],
  76              char fname[GPIV_MAX_CHARS],
  77              TimeCorrPar *time_corr_par
  78              )
  79 /* ----------------------------------------------------------------------------
  80  * Command line argument handling
  81  */
  82 {
  83     char c = '\0';
  84     int argc_next;
  85
  86     while (--argc > 0 && (*++argv)[0] == '-') {
  87
  88 /*
  89  * argc_next is set to 1 if the next cmd line argument has to be searched for;
  90  * in case that the command line argument concerns more than one char or cmd
  91  * line argument needs a parameter
  92  */
  93
  94         argc_next = 0;
  95         while (argc_next == 0 && (c = *++argv[0]))
  96
  97             switch (c) {
  98             case 'v':
  99 #ifdef GIT_HASH
 100                 printf ("git hash: %s\n", GIT_REV);
 101 #else
 102                 printf ("version: %s\n", GPIVTOOLS_VERSION);
 103 #endif
 104                 exit(0);
 105                 break;
 106             case 'h':
 107                 printf("%s\n", argv[0]);
 108                 printf("%s\n",HELP);
 109                 printf("%s\n",USAGE);
 110                 exit(0);
 111                 break;
 112             case 'p':
 113                 print_par = TRUE;
 114                 break;
 115 /*
 116  * First file number N of the series (default: 0)
 117 */
 118             case 'f':
 119                 time_corr_par->first_file = atoi(*++argv);
 120                 time_corr_par->first_file__set = TRUE;
 121                 argc_next = 1;
 122                 --argc;
 123                 break;
 124 /*
 125  * Last file number N of the series (default: 0)
 126 */
 127             case 'l':
 128                 time_corr_par->last_file = atoi(*++argv);
 129                 time_corr_par->last_file__set = TRUE;
 130                 argc_next = 1;
 131                 --argc;
 132                 break;
 133 /*
 134  * X-position (default: 0.0)
 135 */
 136             case 'x':
 137                 time_corr_par->point_x = atof(*++argv);
 138                 time_corr_par->point_x__set = TRUE;
 139                 argc_next = 1;
 140                 --argc;
 141                 break;
 142 /*
 143  * Y-position (default: 0.0)
 144 */
 145             case 'y':
 146                 time_corr_par->point_y = atof(*++argv);
 147                 time_corr_par->point_y__set = TRUE;
 148                 argc_next = 1;
 149                 --argc;
 150                 break;
 151 /*
 152  * long option keys
 153  */
 154             case '-':
 155                 if (strcmp("-help", *argv) == 0) {
 156                     printf("\n%s", argv[0]);
 157                     printf("\n%s", HELP);
 158                     printf("\n%s", USAGE);
 159                     exit(0);
 160                 } else if (strcmp("-print", *argv) == 0) {
 161                     print_par = TRUE;
 162                 } else if (strcmp("-version", *argv) == 0) {
 163 #ifdef GIT_HASH
 164                     printf ("git hash: %s\n", GIT_REV);
 165 #else
 166                     printf ("version: %s\n", GPIVTOOLS_VERSION);
 167 #endif
 168                     exit(0);
 169                 } else if (strcmp("-file_first", *argv) == 0) {
 170                     time_corr_par->first_file = atoi(*++argv);
 171                     time_corr_par->first_file__set = TRUE;
 172                     argc_next = 1;
 173                     --argc;
 174                 } else if (strcmp("-file_last", *argv) == 0) {
 175                     time_corr_par->last_file = atoi(*++argv);
 176                     time_corr_par->last_file__set = TRUE;
 177                     argc_next = 1;
 178                     --argc;
 179                 } else if (strcmp("-point_x", *argv) == 0) {
 180                     time_corr_par->point_x = atof(*++argv);
 181                     time_corr_par->point_x__set = TRUE;
 182                     argc_next = 1;
 183                     --argc;
 184                 } else if (strcmp("-point_y", *argv) == 0) {
 185                     time_corr_par->point_y = atof(*++argv);
 186                     time_corr_par->point_y__set = TRUE;
 187                     argc_next = 1;
 188                     --argc;
 189                 } else {
 190                     gpiv_error("%s: unknown option: %s", argv[0], *argv);
 191                 }
 192                 argc_next = 1;
 193                 break;
 194
 195             default:
 196                 gpiv_error(USAGE);
 197                 break;
 198             }
 199     }
 200
 201     if(argc != 1) {
 202         gpiv_error("%s: %s", argv[0], USAGE);
 203     }
 204
 205     strcpy(fname, argv[0]);
 206 }
 207
 208
 209
 210 void
 211 make_fname_out(char *fname,
 212                char *fname_out
 213                )
 214 /* ----------------------------------------------------------------------------
 215  * define filenames
 216  */
 217 {
 218     gpiv_io_make_fname(fname, ".tco", fname_out);
 219     if (print_par) printf("# Output data file: %s\n",fname_out);
 220
 221 }
 222
 223
 224
 225 void
 226 make_fname_in(char *fname,
 227               int number,
 228               char *fname_in
 229               )
 230 /* ----------------------------------------------------------------------------
 231  * define numbered input filenames
 232  */
 233 {
 234 #define GPIV_EXT_SCALED_PIV  ".sc.piv"          /* Extension of scaled piv file name (ASCII format) */
 235     snprintf(fname_in, GPIV_MAX_CHARS, "%s%d%s", fname, number, GPIV_EXT_SCALED_PIV);
 236     if (print_par) printf("# Input file: %s\n",fname_in);
 237 }
 238
 239
 240
 241 static void
 242 time_corr_par__set (TimeCorrPar *time_corr_par,
 243                     gboolean flag
 244                     )
 245 /* ----------------------------------------------------------------------------
 246  */
 247 {
 248     time_corr_par->first_file__set = flag;
 249     time_corr_par->last_file__set = flag;
 250     time_corr_par->point_x__set = flag;
 251     time_corr_par->point_y__set = flag;
 252 }
 253
 254
 255
 256 int
 257 main(int argc,
 258      char *argv[]
 259      )
 260 /* ----------------------------------------------------------------------------
 261  */
 262 {
 263     FILE *fp = NULL;
 264     gchar *err_msg = NULL;
 265     gchar fname[GPIV_MAX_CHARS],
 266         fname_in[GPIV_MAX_CHARS],
 267         fname_out[GPIV_MAX_CHARS];
 268     guint i, j, k, l, ndata_sets = 0;
 269
 270
 271     gfloat *dx = NULL, *dy = NULL;
 272     gint *flag = NULL;
 273
 274     gfloat mean_dx = 0.0, mean_dy = 0.0, sdev_dx = 0.0, sdev_dy = 0.0;
 275     gint count = 0;
 276     gfloat dx_sum = 0.0, dy_sum = 0.0, dx_sum_q = 0.0, dy_sum_q = 0.0;
 277
 278     gfloat *corr_dx, *corr_dy;
 279     gfloat *dx_sum_cov, *dy_sum_cov;
 280     gint *count_cov;
 281
 282     GpivPivData *in_data = NULL;
 283     TimeCorrPar *time_corr_par = g_new0 (TimeCorrPar, 1);;
 284
 285
 286 /*
 287  * Default parameter values
 288  */
 289     time_corr_par__set (time_corr_par, FALSE);
 290
 291     command_args (argc, argv, fname, time_corr_par);
 292     make_fname_out (fname, fname_out);
 293
 294     if (!time_corr_par->first_file__set) time_corr_par->first_file = 0;
 295     if (!time_corr_par->last_file__set) time_corr_par->last_file = 0;
 296
 297 /*
 298  * Check parameters
 299  */
 300     if (time_corr_par->last_file < time_corr_par->first_file
 301         || (time_corr_par->last_file - time_corr_par->first_file) > MAX_PIVSETS) {
 302         gpiv_error("%s: last_file smaller than first_file or total number of data sets exceeds MAX_PIVSETS=%d",
 303                    argv[0], MAX_PIVSETS);
 304     }
 305
 306
 307     if (print_par) {
 308         g_message ("first_file = %d", time_corr_par->first_file);
 309         g_message ("last_file = %d", time_corr_par->last_file);
 310         g_message ("point_x = %f", time_corr_par->point_x);
 311         g_message ("point_y = %f", time_corr_par->point_y);
 312     }
 313
 314
 315     ndata_sets = time_corr_par->last_file - time_corr_par->first_file + 1;
 316     gpiv_warning("first_file=%d last_file=%d ndata_sets=%d",
 317                  time_corr_par->first_file, time_corr_par->last_file, ndata_sets);
 318     flag = gpiv_ivector (ndata_sets);
 319     dx = gpiv_vector (ndata_sets);
 320     dy = gpiv_vector (ndata_sets);
 321
 322 /*
 323  * Obtaining the data at the sepecified point from all sets
 324  */
 325
 326     for (i = 0; i < ndata_sets; i++) {
 327         make_fname_in (fname, i + time_corr_par->first_file, fname_in);
 328         if (i == 0) {
 329
 330             if ((fp = fopen (fname_in, "r")) == NULL) {
 331                 gpiv_error ("%s: failing opening %s", argv[0], fname_in);
 332             }
 333
 334             if ((in_data = gpiv_read_pivdata (fp)) != NULL)  {
 335                 gpiv_error ("%s: %s", argv[0], err_msg);
 336             }
 337             fclose (fp);
 338
 339 /*
 340  * Filtering specified point
 341  */
 342
 343             for (k = 0; k < in_data->ny; k++) {
 344                 for (l = 0; l < in_data->nx; l++) {
 345                     if (in_data->point_x[k][l] == time_corr_par->point_x
 346                         && in_data->point_y[k][l] == time_corr_par->point_y
 347                         ) {
 348                         flag[i] = in_data->peak_no[k][l];
 349                         dx[i] = in_data->dx[k][l];
 350                         dy[i] = in_data->dy[k][l];
 351                     }
 352                 }
 353             }
 354
 355         }
 356
 357         gpiv_free_pivdata (in_data);
 358     }
 359
 360 /* t_avg()
 361  * Calculate mean and variances
 362  */
 363
 364     for (i = 0; i < ndata_sets; i++) {
 365         if (flag[i] != -1) {
 366             count++;
 367             dx_sum = dx_sum + dx[i];
 368             dy_sum = dy_sum + dy[i];
 369             dx_sum_q = dx_sum_q + dx[i] * dx[i];
 370             dy_sum_q = dy_sum_q + dy[i] * dy[i];
 371         }
 372     }
 373
 374
 375     if (count != 0) {
 376         mean_dx = dx_sum / count;
 377         sdev_dx = sqrt(dx_sum_q / (count - 1)
 378                        - dx_sum * dx_sum / (count * (count -1)));
 379
 380         mean_dy = dy_sum / count;
 381         sdev_dy = sqrt(dy_sum_q / (count - 1)
 382                        - dy_sum * dy_sum / (count * (count -1)));
 383         gpiv_warning("t-corr: mean_dx=%f sdev_dx=%f", mean_dx, sdev_dx);
 384         gpiv_warning("t-corr: mean_dy=%f sdev_dy=%f", mean_dy, sdev_dy);
 385     } else {
 386         gpiv_warning("t-corr: count equal to zero");
 387     }
 388
 389 /* t_corr()
 390  * Calculate correlation
 391  */
 392
 393
 394     if (count != 0) {
 395         count_cov = gpiv_ivector(ndata_sets);
 396         dx_sum_cov = gpiv_vector(ndata_sets);
 397         dy_sum_cov = gpiv_vector(ndata_sets);
 398         corr_dx = gpiv_vector(ndata_sets);
 399         corr_dy = gpiv_vector(ndata_sets);
 400
 401         for (i = 0; i < ndata_sets; i++) {
 402             for (j = 0; j <= i; j++) {
 403                 if (flag[i] != -1 && flag[j] != -1) {
 404 /*                     g_message("i=%d j=%d cor[%d] =>dx[%d]=%f * dx[%d]=%f",  */
 405 /*                               i, j, */
 406 /*                               i - j, */
 407 /*                               i, dx[i] - mean_dx, */
 408 /*                               j, dx[j] - mean_dx */
 409 /*                               ); */
 410
 411                     count_cov[i-j]++;
 412                     dx_sum_cov[i-j] = dx_sum_cov[i-j]
 413                         + (dx[i] - mean_dx) * (dx[j] - mean_dx);
 414                     dy_sum_cov[i-j] = dy_sum_cov[i-j]
 415                         + (dy[i] - mean_dy) * (dy[j] - mean_dy);
 416                 }
 417             }
 418         }
 419
 420 /*
 421  * corr_dx[0] = (2 * sdev_dx * sdev_dx)
 422  * corr_dy[0] = (2 * sdev_dy * sdev_dy)
 423  */
 424
 425         printf("# i   corr_dx\n");
 426         for (i = 0; i < ndata_sets - MIN_SAMPLES; i++) {
 427             corr_dx[i] = dx_sum_cov[i] / dx_sum_cov[0]
 428 /*                 * (float) count_cov[i] / (float) ndata_sets */
 429                 ;
 430             printf("%d  %f\n", i, corr_dx[i]);
 431         }
 432
 433         printf("\n\n# i   corr_dy\n");
 434         for (i = 0; i < ndata_sets - MIN_SAMPLES; i++) {
 435             corr_dy[i] = dy_sum_cov[i] / dy_sum_cov[0]
 436 /*                * (float) count_cov[i] / (float) ndata_sets  */
 437                 ;
 438             printf("%d  %f\n", i, corr_dy[i]);
 439         }
 440
 441         gpiv_free_ivector (count_cov);
 442         gpiv_free_vector  (dx_sum_cov);
 443         gpiv_free_vector  (dy_sum_cov);
 444         gpiv_free_vector  (corr_dx);
 445         gpiv_free_vector  (corr_dy);
 446     }
 447
 448     gpiv_free_vector (dy);
 449     gpiv_free_vector (dx);
 450     gpiv_free_ivector (flag);
 451     exit (0);
 452 }
 453
 454
 455