/* -*- Mode: C; indent-tabs-mode: nil; c-basic-offset: 4 c-style: "K&R" -*- */

/*-----------------------------------------------------------------------------

   suta - subtracts time-avaraged velocities (local mean velocities) from
          the PIV estimators

   Copyright (C) 2002, 2003, 2004 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.  

-----------------------------------------------------------------------------*/

#include <stdio.h> 
#include <stdlib.h> 
#include <assert.h> 
#include <math.h>
#include <gpiv.h>

/* #define PARFILE "scale.par"    */  /* Parameter file name */
#define PARFILE "gpivrc"     /* Parameter file name */
#define MAX_PIVSETS 1000     /* Maximum number of PIV data sets */
#define MIN_SAMPLES 20       /* Minimum number of samples used for estimation */
#define USAGE "\
Usage: suta [-h | --help] [-p | --print] [-v | --version] \n\
filename_ta filename \n\
\n\
keys: \n\
-h | --help:           this on-line help \n\
-p | --print:          print parameters to stdout \n\
-v | --version:        version number \n\
filename               File name of PIV data set (without .piv) \n\
filename_ta            File name of time-averaged (local means) PIV data set \n\
                       (full name) \n\
"

#define HELP  "\
suta - subtracts time-avaraged velocities (local mean) from the PIV estimators"

#define RCSID "$Id: suta.c,v 1.4 2006/01/31 14:18:04 gerber Exp $"

gboolean print_par = FALSE;

void 
command_args(int argc, 
             char *argv[], 
             char fname_in[GPIV_MAX_CHARS],
             char fname_mean[GPIV_MAX_CHARS]
             )
/* ----------------------------------------------------------------------------
 * Command line argument handling
 */
{
    char c = '\0';
    int argc_next;

    while (--argc > 0 && (*++argv)[0] == '-') {

/*
 * argc_next is set to 1 if the next cmd line argument has to be searched for; 
 * in case that the command line argument concerns more than one char or cmd 
 * line argument needs a parameter 
 */

        argc_next = 0;
	while (argc_next == 0 && (c = *++argv[0]))
            
            switch (c) {
            case 'v':
                printf("%s\n", RCSID); 
                exit(0);
                break;
            case 'h':
                printf("%s\n", RCSID); 
                printf("%s\n",HELP);
                printf("%s\n",USAGE);
                exit(0);
                break;
            case 'p':
                print_par = TRUE;
                break;

/*
 * long option keys
 */
	    case '-':
		if (strcmp("-help", *argv) == 0) {
                    printf("\n%s", RCSID);
                    printf("\n%s", HELP);
                    printf("\n%s", USAGE);
                    exit(0);
                } else if (strcmp("-print", *argv) == 0) {
		    print_par = TRUE;
                } else if (strcmp("-version", *argv) == 0) {
                    printf("%s\n", RCSID);
                    exit(0);
                } else {
		    gpiv_error("%s: unknown option: %s", RCSID, *argv);
		}
		argc_next = 1;
		break;
                
            default:
                gpiv_error(USAGE);
                break;
            }
    }

    if(argc != 2) { 
        gpiv_error("%s: %s", RCSID, USAGE);
    }

    strcpy(fname_mean, argv[0]);
    strcpy(fname_in, argv[1]);
}




void 
make_fname(char *fname, 
           char *fname_in,
           char *fname_out
           )
/* ----------------------------------------------------------------------------
 * define input and output filenames
 */
{     
    gpiv_io_make_fname(fname, GPIV_EXT_PIV, fname_in);
    if (print_par) printf("# Input data file: %s\n",fname_in);

/* 
 * Extension of filename with time-averaged reduced piv data
 */
#define EXT_SUB_TAVG_PIV  ".suta.piv" 
    gpiv_io_make_fname(fname, EXT_SUB_TAVG_PIV, fname_out);
    if (print_par) printf("# Output data file: %s\n",fname_out);

}



int 
main(int argc, 
     char *argv[]
     )
/* ----------------------------------------------------------------------------
 */
{
/*     FILE *fp_out = NULL; */
    gchar *err_msg = NULL;
    gchar fname[GPIV_MAX_CHARS],
        fname_in[GPIV_MAX_CHARS],
        fname_out[GPIV_MAX_CHARS],
        fname_mean[GPIV_MAX_CHARS];
    gint i, j, scale = 0;
    
    gchar  c_line[GPIV_MAX_LINES_C][GPIV_MAX_CHARS];
    gint nc_lines = 0;

    GpivPivData in_data, mean_data, out_data;

    in_data.point_x = NULL;
    in_data.point_y = NULL;
    in_data.dx = NULL;
    in_data.dy = NULL;
    in_data.snr = NULL;
    in_data.peak_no = NULL;

    mean_data.point_x = NULL;
    mean_data.point_y = NULL;
    mean_data.dx = NULL;
    mean_data.dy = NULL;
    mean_data.snr = NULL;
    mean_data.peak_no = NULL;

    out_data.point_x = NULL;
    out_data.point_y = NULL;
    out_data.dx = NULL;
    out_data.dy = NULL;
    out_data.snr = NULL;
    out_data.peak_no = NULL;



    command_args(argc, argv, fname, fname_mean);
    make_fname(fname, fname_in, fname_out);
    if (print_par) {
        g_message("fname_mean = %s", fname_mean);
    }


    if ((scale = gpiv_fcount_pivdata(fname_in, &in_data)) == -1) {
        gpiv_error("%s: Failure calling gpiv_count_pivdata", RCSID);
        return 1;
    }

    mean_data.nx = in_data.nx;
    mean_data.ny = in_data.ny;
    out_data.nx = in_data.nx;
    out_data.ny = in_data.ny;
    
    gpiv_alloc_pivdata(&in_data);
    gpiv_alloc_pivdata(&mean_data);
    gpiv_alloc_pivdata(&out_data);

    if ((err_msg = 
         gpiv_fread_pivdata(fname_in, &in_data, c_line, &nc_lines)) 
        != NULL)  gpiv_error ("%s: %s", RCSID, err_msg);

    if ((err_msg = 
         gpiv_fread_pivdata(fname_mean, &mean_data, c_line, &nc_lines)) 
        != NULL)  gpiv_error ("%s: %s", RCSID, err_msg);

    for (i = 0; i < in_data.ny; i++) {
        for (j = 0; j < in_data.nx; j++) {
            out_data.point_x[i][j] = in_data.point_x[i][j];
            out_data.point_y[i][j] = in_data.point_y[i][j];
            out_data.peak_no[i][j] = in_data.peak_no[i][j];
            if (mean_data.peak_no[i][j] > 0) {
                out_data.snr[i][j] = in_data.snr[i][j];
                out_data.dx[i][j] = in_data.dx[i][j] - mean_data.dx[i][j];
                out_data.dy[i][j] = in_data.dy[i][j] - mean_data.dy[i][j];
            } else {
                out_data.snr[i][j] = GPIV_SNR_DISABLE;
                out_data.dx[i][j] = 0.0;
                out_data.dy[i][j] = 0.0;
            }
        }
    }
    
    if ((err_msg =
         gpiv_fwrite_pivdata(fname_out, &out_data, c_line,
                             nc_lines, scale, RCSID))
        != NULL)  gpiv_error ("%s: %s", RCSID, err_msg);
    
    gpiv_free_pivdata(&in_data);
    gpiv_free_pivdata(&mean_data);
    gpiv_free_pivdata(&out_data);
    exit(0);
}