src/post/vorstra.c

   1 /* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */
   2
   3 /*---------------------------------------------------------------------------
   4
   5   vorstra - calculates differential quantities (like vorticity and
   6   strain) from PIV data
   7
   8
   9    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008
  10    Gerber van der Graaf <gerber_graaf@users.sourceforge.net
  11
  12    This program is free software; you can redistribute it and/or modify
  13    it under the terms of the GNU General Public License as published by
  14    the Free Software Foundation; either version 2, or (at your option)
  15    any later version.
  16
  17    This program is distributed in the hope that it will be useful,
  18    but WITHOUT ANY WARRANTY; without even the implied warranty of
  19    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20    GNU General Public License for more details.
  21
  22    You should have received a copy of the GNU General Public License
  23    along with this program; if not, write to the Free Software Foundation,
  24    Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  25
  26 ----------------------------------------------------------------------------
  27
  28   Some formulea:
  29      vorty_z = dV/dx - dU/dy
  30      sh_strain = dU/dy + dV/dx
  31       n_strain = dU/dx + dV/dy
  32
  33
  34 -----------------------------------------------------------------------------*/
  35
  36 #include <stdio.h>
  37 #include <stdlib.h>
  38 #include <math.h>
  39 #include <glib.h>
  40 #include <gpiv.h>
  41
  42 /* #define PARFILE "vorstra.par" */     /* Parameter file name */
  43 #define PARFILE "gpivrc"        /* Parameter file name */
  44
  45 #define GNUPLOT_DISPLAY_COLOR "DarkBlue"
  46 #define GNUPLOT_DISPLAY_SIZE 250
  47
  48 #define USAGE "\
  49 Usage: vorstra | vorty | nstrain | sstrain \n\
  50                [-d int] [-g] [--no_g] [-h | --help] [-o] [-s] [-n] \n\
  51                [-p | --print] [-v | --version] [filename] < stdin > stdout \n\
  52 \n\
  53 keys: \n\
  54 -d N:                  differential type; central (0), least squares \n\
  55                        (1), richardson (2), circulation method (3)  \n\
  56 -g:                    graphical visualization with gnuplot (needs -f) \n\
  57 --no_g:                suppresses graphical visualization \n\
  58 -h | --help:           this on-line help \n\
  59 -n:                    normal strain analyses \n\
  60 -o:                    vorticity analyses \n\
  61 -p | --print:          print parameters to stdout  \n\
  62 -s:                    shear strain analyses \n\
  63 -v | --version;        version number\n\
  64 filename:              input PIV data file,. Substitutes stdin and stdout \n\
  65 "
  66
  67 #ifdef DEBUG
  68 #define USAGE_DEBUG "\
  69 Developers version also contains: \n\
  70                [-p_main][-p_differential] \n\
  71 \n\
  72 keys: \n\
  73 -p_'function' N: prints data to be generated in the function; the \n\
  74                  higher N, the more detailed the output. \n\
  75                  For N = 10, err_vec will exit at the end of the function"
  76 #endif
  77
  78
  79 #define HELP  "\
  80 vorstra calculates the differential quantities vorticity, shear strain and \n\
  81 normal strain from PIV data."
  82
  83 #define RCSID "$Id: vorstra.c,v 2.16 2008-10-09 14:45:01 gerber Exp $"
  84
  85
  86 gboolean use_stdin_stdout = FALSE;
  87 gboolean verbose = FALSE;
  88
  89 int gnuplot = 0;
  90 gboolean gnuplot__set = FALSE;
  91
  92 #ifdef DEBUG
  93 /*
  94  * Parameters for development version
  95  */
  96 int print_main = 0, print_differential = 0;
  97 #endif
  98
  99
 100
 101 static void
 102 command_args(int argc,
 103              char *argv[],
 104              char fname[GPIV_MAX_CHARS],
 105              GpivPostPar * piv_post_par
 106              )
 107 /* ----------------------------------------------------------------------------
 108  * Command line argument handling
 109  */
 110 {
 111     char c;
 112     int argc_next;
 113
 114
 115     while (--argc > 0 && (*++argv)[0] == '-') {
 116         argc_next = 0;
 117
 118 /*
 119  * argc_next is set to 1 if the next cmd line argument has to be
 120  * searched for; in case that the command line argument concerns more
 121  * than one char or cmd line argument needs a parameter
 122  */
 123         while (argc_next == 0 && (c = *++argv[0])) {
 124             switch (c) {
 125 /*
 126  * Use Revision Control System (RCS) for version
 127  */
 128             case 'v':
 129                 printf("\n%s\n", RCSID);
 130                 exit(0);
 131
 132 /*
 133  * differential operator type
 134  */
 135             case 'd':
 136                 piv_post_par->diff_type__set = TRUE;
 137                 piv_post_par->diff_type = atoi(*++argv);
 138                 argc_next = 1;
 139                 --argc;
 140                 break;
 141
 142 /*
 143  * graphic output with gnuplot
 144  */
 145             case 'g':
 146               gnuplot = 1;
 147               gnuplot__set = TRUE;
 148               break;
 149             case 'h':
 150                 printf("\n%s", argv[0]);
 151                 printf("\n%s", HELP);
 152                 printf("\n%s\n", USAGE);
 153                 exit(0);
 154
 155 /*
 156  * normal strain analysis
 157  */
 158             case 'n':
 159                 piv_post_par->operator_vorstra__set = TRUE;
 160                 piv_post_par->operator_vorstra = GPIV_N_STRAIN;
 161                 break;
 162
 163 /*
 164  * vorticity analysis
 165  */
 166             case 'o':
 167                 piv_post_par->operator_vorstra__set = TRUE;
 168                 piv_post_par->operator_vorstra = GPIV_VORTICITY;
 169                 break;
 170             case 'p':
 171 #ifdef DEBUG
 172                 if (strcmp(*argv, "p_main") != 0) {
 173 #endif
 174                     verbose = TRUE;
 175 #ifdef DEBUG
 176                 } else if (strcmp("p_main", *argv) == 0) {
 177                     print_main = atoi(*++argv);
 178                     --argc;
 179                     argc_next = 1;
 180                 }
 181 #endif
 182                 break;
 183 /*
 184  * shear strain analysis
 185  */
 186             case 's':
 187                 piv_post_par->operator_vorstra__set = TRUE;
 188                 piv_post_par->operator_vorstra = GPIV_S_STRAIN;
 189                 break;
 190
 191 /*
 192  * long option keys
 193  */
 194             case '-':
 195                 if (strcmp("-help", *argv) == 0) {
 196                     printf("\n%s", argv[0]);
 197                     printf("\n%s", HELP);
 198                     printf("\n%s", USAGE);
 199                     exit(0);
 200                 } else if (strcmp("-print", *argv) == 0) {
 201                     verbose = TRUE;
 202                 } else if (strcmp("-version", *argv) == 0) {
 203                     printf("%s\n", RCSID);
 204                     exit(0);
 205
 206 /*
 207  * do not plot with gnuplot
 208  */
 209                 } else if (strcmp("-no_g", *argv) == 0) {
 210                     gnuplot = 0;
 211                     gnuplot__set = TRUE;
 212                     argc_next = 1;
 213
 214                 } else {
 215                     gpiv_error("%s: unknown option: %s", argv[0], *argv);
 216                 }
 217                 argc_next = 1;
 218                 break;
 219
 220             default:
 221                 gpiv_error("\n%s error : unknown option: %s\n",
 222                         argv[0], *argv);
 223                 break;
 224             }
 225         }
 226     }
 227
 228     /*
 229      * Check if filename or stdin /stdout is used
 230      */
 231     if (argc == 1) {
 232         use_stdin_stdout = FALSE;
 233         strcpy(fname, argv[argc - 1]);
 234     } else if (argc == 0) {
 235         use_stdin_stdout = TRUE;
 236         verbose = FALSE;
 237     } else {
 238         gpiv_error("\n%s error: unknown argument: %s\n", argv[0], *argv);
 239     }
 240
 241 }
 242
 243
 244
 245 static gchar *
 246 make_fname(char *fname_in,
 247            char *fname_parameter,
 248            char *fname_out,
 249            GpivPostPar *piv_post_par
 250            )
 251 /* ----------------------------------------------------------------------------
 252  * generate filenames
 253  */
 254 {
 255     gchar *err_msg = NULL;
 256     gchar *fname_base = NULL;
 257
 258     if (fname_in == NULL ) {
 259         err_msg = "make_fname: \"fname_in == NULL\"";
 260         return (err_msg);
 261     }
 262
 263 /*
 264  * Stripping filename
 265  */
 266     fname_base = g_strdup(fname_in);
 267     strtok(fname_base, ".");
 268
 269 /*
 270  * filenames for output
 271  */
 272     gpiv_io_make_fname(fname_base, GPIV_EXT_PAR, fname_parameter);
 273     if (verbose)
 274         printf("# Parameter file: %s\n", fname_parameter);
 275
 276     if (piv_post_par->operator_vorstra == GPIV_VORTICITY) {
 277         gpiv_io_make_fname (fname_base, GPIV_EXT_VOR, fname_out);
 278     } else if (piv_post_par->operator_vorstra == GPIV_S_STRAIN) {
 279         gpiv_io_make_fname (fname_base, GPIV_EXT_SSTR, fname_out);
 280     } else if (piv_post_par->operator_vorstra == GPIV_N_STRAIN) {
 281         gpiv_io_make_fname (fname_base, GPIV_EXT_NSTR, fname_out);
 282     } else {
 283         gpiv_error("%s: non valid operation", RCSID);
 284     }
 285     if (verbose)
 286         printf("# output data file: %s\n", fname_out);
 287
 288     g_free (fname_base);
 289     return (err_msg);
 290 }
 291
 292
 293
 294 int
 295 main(int argc,
 296      char *argv[]
 297      )
 298 /*-----------------------------------------------------------------------------
 299 */
 300 {
 301     char *err_msg = NULL, *c = NULL;
 302     FILE *fp_par_dat, *fp;
 303     gchar fname_in[GPIV_MAX_CHARS],
 304         fname_out[GPIV_MAX_CHARS],
 305         fname_parameter[GPIV_MAX_CHARS];
 306
 307     GpivPostPar *piv_post_par = g_new (GpivPostPar, 1);;
 308     GpivPivData *in_data = NULL;
 309     GpivScalarData *out_data = NULL;
 310
 311
 312     verbose = FALSE;
 313     use_stdin_stdout = FALSE;
 314     piv_post_par->operator_vorstra__set = FALSE;
 315     piv_post_par->diff_type__set = FALSE;
 316
 317 /*
 318  * Define GpivOperation type from program name, which is a symbolic link to
 319  * gpiv_vorstra
 320  */
 321     if ((c = strstr(argv[0], "vorstra")) != NULL) {
 322
 323     } else if ((c = strstr(argv[0], "vorty")) != NULL) {
 324         piv_post_par->operator_vorstra = GPIV_VORTICITY;
 325         piv_post_par->operator_vorstra__set = TRUE;
 326
 327    } else if ((c = strstr(argv[0], "nstrain")) != NULL) {
 328         piv_post_par->operator_vorstra = GPIV_N_STRAIN;
 329         piv_post_par->operator_vorstra__set = TRUE;
 330
 331     } else if ((c = strstr(argv[0], "sstrain")) != NULL) {
 332         piv_post_par->operator_vorstra = GPIV_S_STRAIN;
 333         piv_post_par->operator_vorstra__set = TRUE;
 334
 335     } else {
 336         gpiv_error("vorstra: unvalid program name or symlink");
 337     }
 338
 339
 340     command_args (argc, argv, fname_in, piv_post_par);
 341     if (verbose) {
 342       printf("# %s\n# Command line options:\n", RCSID);
 343       gpiv_post_print_parameters (NULL, piv_post_par);
 344     }
 345
 346
 347     if (use_stdin_stdout == FALSE) {
 348         make_fname(fname_in, fname_parameter, fname_out, piv_post_par);
 349
 350 /*
 351  * Prints command line parameters to par-file
 352  */
 353         if ((fp_par_dat = fopen(fname_parameter, "a")) == NULL)
 354              gpiv_error("\n%s: failure opening %s for input",
 355                         argv[0], fname_parameter);
 356       fprintf(fp_par_dat, "# %s\n# Command line options:\n", argv[0]);
 357       gpiv_post_print_parameters (fp_par_dat, piv_post_par);
 358
 359 /*
 360  * Reading parametes from PARFILE (and writing to data par-file)
 361  */
 362         gpiv_scan_parameter (GPIV_POSTPAR_KEY, PARFILE, piv_post_par, verbose);
 363         if ((err_msg =
 364              gpiv_scan_resourcefiles (GPIV_POSTPAR_KEY, piv_post_par, verbose))
 365             != NULL)  gpiv_error ("%s: %s", argv[0], err_msg);
 366         gpiv_post_print_parameters (fp_par_dat, piv_post_par);
 367         fclose (fp_par_dat);
 368
 369
 370     } else {
 371         gpiv_scan_parameter (GPIV_POSTPAR_KEY, PARFILE, piv_post_par, verbose);
 372         if ((err_msg =
 373              gpiv_scan_resourcefiles (GPIV_POSTPAR_KEY, piv_post_par, verbose))
 374             != NULL)  gpiv_error ("%s: %s", argv[0], err_msg);
 375     }
 376
 377     gpiv_post_check_parameters_read (piv_post_par, NULL);
 378
 379
 380 /*
 381  * Check parameters on correct values and adjust belonging variables
 382  */
 383     if (piv_post_par->diff_type__set) {
 384          if (piv_post_par->diff_type != 0
 385              && piv_post_par->diff_type != 1
 386              && piv_post_par->diff_type != 2
 387              && piv_post_par->diff_type != 3)
 388               gpiv_error("%s error: no value differential differentiator type",
 389                          argv[0]);
 390     }
 391
 392     if (use_stdin_stdout == TRUE && gnuplot == 1) {
 393         gpiv_error ("%s: filename has to be used in combination with 'gnuplot'\n",
 394                     argv[0]);
 395     }
 396
 397     if ((piv_post_par->operator_vorstra == GPIV_S_STRAIN
 398          || piv_post_par->operator_vorstra == GPIV_N_STRAIN)
 399         && piv_post_par->diff_type == GPIV_CIRCULATION) {
 400         gpiv_error ("\n%s error: strain can not be calculated by circulation method",
 401                     argv[0]);
 402     }
 403
 404
 405 /*
 406  * Reading input PIV data
 407  */
 408     if (use_stdin_stdout == TRUE) {
 409         fp = stdin;
 410     } else {
 411         if ((fp = fopen (fname_in, "rb")) == NULL) {
 412             gpiv_error ("%s: Failure opening %s for input", argv[0], fname_in);
 413         }
 414     }
 415
 416     if ((in_data = gpiv_read_pivdata (fp)) == NULL) {
 417         gpiv_error ("%s: %s", argv[0], err_msg);
 418     }
 419     if (use_stdin_stdout == FALSE) fclose (fp);
 420
 421
 422 /*
 423  * here the function calls of post-processing; calculating vorticity, strain
 424  */
 425
 426     if ((out_data = gpiv_post_vorstra (in_data, piv_post_par)) == NULL) {
 427         gpiv_error ("%s: %s", argv[0], err_msg);
 428     }
 429
 430 /*
 431  * Adding comment to the data
 432  * And writing to output
 433  */
 434     g_free (out_data->comment);
 435     out_data->comment = g_strdup_printf ("# Software: %s\n", RCSID);
 436     out_data->comment = gpiv_add_datetime_to_comment (out_data->comment);
 437     if (piv_post_par->operator_vorstra == GPIV_VORTICITY) {
 438         out_data->comment = g_strconcat (out_data->comment,
 439                                          "\n# Data type: vorticity [1/s]\n",
 440                                          NULL);
 441     } else  if (piv_post_par->operator_vorstra == GPIV_N_STRAIN) {
 442         out_data->comment = g_strconcat (out_data->comment,
 443                                          "\n# Data type: normal strain [1/s]\n",
 444                                          NULL);
 445     } else  if (piv_post_par->operator_vorstra == GPIV_S_STRAIN) {
 446         out_data->comment = g_strconcat (out_data->comment,
 447                                          "\n# Data type: shear strain [1/s]\n",
 448                                          NULL);
 449     }
 450
 451
 452     if (use_stdin_stdout == TRUE) {
 453         fp = stdout;
 454     } else {
 455         if ((fp = fopen (fname_out, "wb")) == NULL) {
 456             gpiv_error ("%s: Failure opening %s for output", argv[0], fname_out);
 457         }
 458     }
 459
 460     if (gnuplot == 1) {
 461         if ((err_msg = gpiv_write_sc_griddata (fp, out_data, TRUE)) != NULL) {
 462             gpiv_error ("%s: %s", argv[0], err_msg);
 463         }
 464     } else {
 465         if ((err_msg = gpiv_write_scdata (fp, out_data, TRUE)) != NULL) {
 466             gpiv_error ("%s: %s", argv[0], err_msg);
 467         }
 468     }
 469     if (use_stdin_stdout == FALSE) fclose (fp);
 470
 471
 472 /*
 473  * Freeing allocated memory of matrices
 474  */
 475     gpiv_free_pivdata(in_data);
 476
 477
 478 /*
 479  * Graphical output with gnuplot
 480  */
 481     if (gnuplot) {
 482         gchar title[GPIV_MAX_CHARS];
 483
 484         snprintf (title, GPIV_MAX_CHARS, "vorticity of %s", fname_in);
 485         gpiv_scalar_gnuplot (fname_out, title, GNUPLOT_DISPLAY_COLOR,
 486                             GNUPLOT_DISPLAY_SIZE);
 487     }
 488
 489     exit (0);
 490 }
 491
 492
 493