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
  33 /* #define PARFILE "scale.par"    */  /* Parameter file name */
  34 #define PARFILE "gpivrc"     /* Parameter file name */
  35 #define MAX_PIVSETS 1000     /* Maximum number of PIV data sets */
  36 #define MIN_SAMPLES 20       /* Minimum number of samples used for estimation */
  37 #define USAGE "\
  38 Usage: t-corr [-h | --help] [-p | --print] [-v | --version] \n\
  39 [-f | --file_first N]  [-l | --file_last N] [-x | --point_x F] \n\
  40 [-y | --point_y F] file_basename\n\
  41 \n\
  42 keys: \n\
  43 -h | --help:           this on-line help \n\
  44 -p | --print:          print parameters to stdout \n\
  45 -v | --version:        version number \n\
  46 -f | --file_first N    First file number of the series (default: 0) \n\
  47 -l | --file_last N     Last file number of the series (default: 0) \n\
  48 -x | --point_x F       x-position of data set to be analysed (default: 0.0) \n\
  49 -y | --point_y F       y-position of data set to be analysed (default: 0.0) \n\
  50 file_basename          File basename of a numbered series of PIV data sets \n\
  51 "
  52
  53 #define HELP  "\
  54 t-corr - calculates the Eulerian correlation from a series of PIV data sets"
  55
  56 #define RCSID "$Id: t-corr.c,v 1.6 2008-09-25 13:08:34 gerber Exp $"
  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                 printf("%s\n", RCSID);
 100                 exit(0);
 101                 break;
 102             case 'h':
 103                 printf("%s\n", argv[0]);
 104                 printf("%s\n",HELP);
 105                 printf("%s\n",USAGE);
 106                 exit(0);
 107                 break;
 108             case 'p':
 109                 print_par = TRUE;
 110                 break;
 111 /*
 112  * First file number N of the series (default: 0)
 113 */
 114             case 'f':
 115                 time_corr_par->first_file = atoi(*++argv);
 116                 time_corr_par->first_file__set = TRUE;
 117                 argc_next = 1;
 118                 --argc;
 119                 break;
 120 /*
 121  * Last file number N of the series (default: 0)
 122 */
 123             case 'l':
 124                 time_corr_par->last_file = atoi(*++argv);
 125                 time_corr_par->last_file__set = TRUE;
 126                 argc_next = 1;
 127                 --argc;
 128                 break;
 129 /*
 130  * X-position (default: 0.0)
 131 */
 132             case 'x':
 133                 time_corr_par->point_x = atof(*++argv);
 134                 time_corr_par->point_x__set = TRUE;
 135                 argc_next = 1;
 136                 --argc;
 137                 break;
 138 /*
 139  * Y-position (default: 0.0)
 140 */
 141             case 'y':
 142                 time_corr_par->point_y = atof(*++argv);
 143                 time_corr_par->point_y__set = TRUE;
 144                 argc_next = 1;
 145                 --argc;
 146                 break;
 147 /*
 148  * long option keys
 149  */
 150             case '-':
 151                 if (strcmp("-help", *argv) == 0) {
 152                     printf("\n%s", argv[0]);
 153                     printf("\n%s", HELP);
 154                     printf("\n%s", USAGE);
 155                     exit(0);
 156                 } else if (strcmp("-print", *argv) == 0) {
 157                     print_par = TRUE;
 158                 } else if (strcmp("-version", *argv) == 0) {
 159                     printf("%s\n", RCSID);
 160                     exit(0);
 161                 } else if (strcmp("-file_first", *argv) == 0) {
 162                     time_corr_par->first_file = atoi(*++argv);
 163                     time_corr_par->first_file__set = TRUE;
 164                     argc_next = 1;
 165                     --argc;
 166                 } else if (strcmp("-file_last", *argv) == 0) {
 167                     time_corr_par->last_file = atoi(*++argv);
 168                     time_corr_par->last_file__set = TRUE;
 169                     argc_next = 1;
 170                     --argc;
 171                 } else if (strcmp("-point_x", *argv) == 0) {
 172                     time_corr_par->point_x = atof(*++argv);
 173                     time_corr_par->point_x__set = TRUE;
 174                     argc_next = 1;
 175                     --argc;
 176                 } else if (strcmp("-point_y", *argv) == 0) {
 177                     time_corr_par->point_y = atof(*++argv);
 178                     time_corr_par->point_y__set = TRUE;
 179                     argc_next = 1;
 180                     --argc;
 181                 } else {
 182                     gpiv_error("%s: unknown option: %s", argv[0], *argv);
 183                 }
 184                 argc_next = 1;
 185                 break;
 186
 187             default:
 188                 gpiv_error(USAGE);
 189                 break;
 190             }
 191     }
 192
 193     if(argc != 1) {
 194         gpiv_error("%s: %s", argv[0], USAGE);
 195     }
 196
 197     strcpy(fname, argv[0]);
 198 }
 199
 200
 201
 202 void
 203 make_fname_out(char *fname,
 204                char *fname_out
 205                )
 206 /* ----------------------------------------------------------------------------
 207  * define filenames
 208  */
 209 {
 210     gpiv_io_make_fname(fname, ".tco", fname_out);
 211     if (print_par) printf("# Output data file: %s\n",fname_out);
 212
 213 }
 214
 215
 216
 217 void
 218 make_fname_in(char *fname,
 219               int number,
 220               char *fname_in
 221               )
 222 /* ----------------------------------------------------------------------------
 223  * define numbered input filenames
 224  */
 225 {
 226 #define GPIV_EXT_SCALED_PIV  ".sc.piv"          /* Extension of scaled piv file name (ASCII format) */
 227     snprintf(fname_in, GPIV_MAX_CHARS, "%s%d%s", fname, number, GPIV_EXT_SCALED_PIV);
 228     if (print_par) printf("# Input file: %s\n",fname_in);
 229 }
 230
 231
 232
 233 static void
 234 time_corr_par__set (TimeCorrPar *time_corr_par,
 235                     gboolean flag
 236                     )
 237 /* ----------------------------------------------------------------------------
 238  */
 239 {
 240     time_corr_par->first_file__set = flag;
 241     time_corr_par->last_file__set = flag;
 242     time_corr_par->point_x__set = flag;
 243     time_corr_par->point_y__set = flag;
 244 }
 245
 246
 247
 248 int
 249 main(int argc,
 250      char *argv[]
 251      )
 252 /* ----------------------------------------------------------------------------
 253  */
 254 {
 255     FILE *fp = NULL;
 256     gchar *err_msg = NULL;
 257     gchar fname[GPIV_MAX_CHARS],
 258         fname_in[GPIV_MAX_CHARS],
 259         fname_out[GPIV_MAX_CHARS];
 260     guint i, j, k, l, ndata_sets = 0;
 261
 262
 263     gfloat *dx = NULL, *dy = NULL;
 264     gint *flag = NULL;
 265
 266     gfloat mean_dx = 0.0, mean_dy = 0.0, sdev_dx = 0.0, sdev_dy = 0.0;
 267     gint count = 0;
 268     gfloat dx_sum = 0.0, dy_sum = 0.0, dx_sum_q = 0.0, dy_sum_q = 0.0;
 269
 270     gfloat *corr_dx, *corr_dy;
 271     gfloat *dx_sum_cov, *dy_sum_cov;
 272     gint *count_cov;
 273
 274     GpivPivData *in_data = NULL;
 275     TimeCorrPar *time_corr_par = g_new0 (TimeCorrPar, 1);;
 276
 277
 278 /*
 279  * Default parameter values
 280  */
 281     time_corr_par__set (time_corr_par, FALSE);
 282
 283     command_args (argc, argv, fname, time_corr_par);
 284     make_fname_out (fname, fname_out);
 285
 286     if (!time_corr_par->first_file__set) time_corr_par->first_file = 0;
 287     if (!time_corr_par->last_file__set) time_corr_par->last_file = 0;
 288
 289 /*
 290  * Check parameters
 291  */
 292     if (time_corr_par->last_file < time_corr_par->first_file
 293         || (time_corr_par->last_file - time_corr_par->first_file) > MAX_PIVSETS) {
 294         gpiv_error("%s: last_file smaller than first_file or total number of data sets exceeds MAX_PIVSETS=%d",
 295                    argv[0], MAX_PIVSETS);
 296     }
 297
 298
 299     if (print_par) {
 300         g_message ("first_file = %d", time_corr_par->first_file);
 301         g_message ("last_file = %d", time_corr_par->last_file);
 302         g_message ("point_x = %f", time_corr_par->point_x);
 303         g_message ("point_y = %f", time_corr_par->point_y);
 304     }
 305
 306
 307     ndata_sets = time_corr_par->last_file - time_corr_par->first_file + 1;
 308     gpiv_warning("first_file=%d last_file=%d ndata_sets=%d",
 309                  time_corr_par->first_file, time_corr_par->last_file, ndata_sets);
 310     flag = gpiv_ivector (ndata_sets);
 311     dx = gpiv_vector (ndata_sets);
 312     dy = gpiv_vector (ndata_sets);
 313
 314 /*
 315  * Obtaining the data at the sepecified point from all sets
 316  */
 317
 318     for (i = 0; i < ndata_sets; i++) {
 319         make_fname_in (fname, i + time_corr_par->first_file, fname_in);
 320         if (i == 0) {
 321
 322             if ((fp = fopen (fname_in, "r")) == NULL) {
 323                 gpiv_error ("%s: failing opening %s", argv[0], fname_in);
 324             }
 325
 326             if ((in_data = gpiv_read_pivdata (fp)) != NULL)  {
 327                 gpiv_error ("%s: %s", argv[0], err_msg);
 328             }
 329             fclose (fp);
 330
 331 /*
 332  * Filtering specified point
 333  */
 334
 335             for (k = 0; k < in_data->ny; k++) {
 336                 for (l = 0; l < in_data->nx; l++) {
 337                     if (in_data->point_x[k][l] == time_corr_par->point_x
 338                         && in_data->point_y[k][l] == time_corr_par->point_y
 339                         ) {
 340                         flag[i] = in_data->peak_no[k][l];
 341                         dx[i] = in_data->dx[k][l];
 342                         dy[i] = in_data->dy[k][l];
 343                     }
 344                 }
 345             }
 346
 347         }
 348
 349         gpiv_free_pivdata (in_data);
 350     }
 351
 352 /* t_avg()
 353  * Calculate mean and variances
 354  */
 355
 356     for (i = 0; i < ndata_sets; i++) {
 357         if (flag[i] != -1) {
 358             count++;
 359             dx_sum = dx_sum + dx[i];
 360             dy_sum = dy_sum + dy[i];
 361             dx_sum_q = dx_sum_q + dx[i] * dx[i];
 362             dy_sum_q = dy_sum_q + dy[i] * dy[i];
 363         }
 364     }
 365
 366
 367     if (count != 0) {
 368         mean_dx = dx_sum / count;
 369         sdev_dx = sqrt(dx_sum_q / (count - 1)
 370                        - dx_sum * dx_sum / (count * (count -1)));
 371
 372         mean_dy = dy_sum / count;
 373         sdev_dy = sqrt(dy_sum_q / (count - 1)
 374                        - dy_sum * dy_sum / (count * (count -1)));
 375         gpiv_warning("t-corr: mean_dx=%f sdev_dx=%f", mean_dx, sdev_dx);
 376         gpiv_warning("t-corr: mean_dy=%f sdev_dy=%f", mean_dy, sdev_dy);
 377     } else {
 378         gpiv_warning("t-corr: count equal to zero");
 379     }
 380
 381 /* t_corr()
 382  * Calculate correlation
 383  */
 384
 385
 386     if (count != 0) {
 387         count_cov = gpiv_ivector(ndata_sets);
 388         dx_sum_cov = gpiv_vector(ndata_sets);
 389         dy_sum_cov = gpiv_vector(ndata_sets);
 390         corr_dx = gpiv_vector(ndata_sets);
 391         corr_dy = gpiv_vector(ndata_sets);
 392
 393         for (i = 0; i < ndata_sets; i++) {
 394             for (j = 0; j <= i; j++) {
 395                 if (flag[i] != -1 && flag[j] != -1) {
 396 /*                     g_message("i=%d j=%d cor[%d] =>dx[%d]=%f * dx[%d]=%f",  */
 397 /*                               i, j, */
 398 /*                               i - j, */
 399 /*                               i, dx[i] - mean_dx, */
 400 /*                               j, dx[j] - mean_dx */
 401 /*                               ); */
 402
 403                     count_cov[i-j]++;
 404                     dx_sum_cov[i-j] = dx_sum_cov[i-j]
 405                         + (dx[i] - mean_dx) * (dx[j] - mean_dx);
 406                     dy_sum_cov[i-j] = dy_sum_cov[i-j]
 407                         + (dy[i] - mean_dy) * (dy[j] - mean_dy);
 408                 }
 409             }
 410         }
 411
 412 /*
 413  * corr_dx[0] = (2 * sdev_dx * sdev_dx)
 414  * corr_dy[0] = (2 * sdev_dy * sdev_dy)
 415  */
 416
 417         printf("# i   corr_dx\n");
 418         for (i = 0; i < ndata_sets - MIN_SAMPLES; i++) {
 419             corr_dx[i] = dx_sum_cov[i] / dx_sum_cov[0]
 420 /*                 * (float) count_cov[i] / (float) ndata_sets */
 421                 ;
 422             printf("%d  %f\n", i, corr_dx[i]);
 423         }
 424
 425         printf("\n\n# i   corr_dy\n");
 426         for (i = 0; i < ndata_sets - MIN_SAMPLES; i++) {
 427             corr_dy[i] = dy_sum_cov[i] / dy_sum_cov[0]
 428 /*                * (float) count_cov[i] / (float) ndata_sets  */
 429                 ;
 430             printf("%d  %f\n", i, corr_dy[i]);
 431         }
 432
 433         gpiv_free_ivector (count_cov);
 434         gpiv_free_vector  (dx_sum_cov);
 435         gpiv_free_vector  (dy_sum_cov);
 436         gpiv_free_vector  (corr_dx);
 437         gpiv_free_vector  (corr_dy);
 438     }
 439
 440     gpiv_free_vector (dy);
 441     gpiv_free_vector (dx);
 442     gpiv_free_ivector (flag);
 443     exit (0);
 444 }
 445
 446
 447