/**---------------------------------------------------------------------------
 ** 
 ** vector.c -- Manipulate data in vectors
 ** 
 ** Author          : Pierre Jaccard
 ** Created On      : 1999/07/15 22:05:37
 ** Last Modified By: Pierre Jaccard
 ** Last Modified On: 1999/07/16 20:09:53
 ** Update Count    : 2
 ** Directory       : /home/pego/pcd1/codas3c/gfi/src/libs/math/
 ** Version         : 0.0
 ** Status          : Unknown
 ** ---------------------------------------------------------------------- ** 
 ** DESCRIPTION: 
 ** 
 **    Various data manipulation in vectors.
 ** 
 ** ---------------------------------------------------------------------- ** 
 ** REVISIONS: 
 ** ---------------------------------------------------------------------- ** 
 ** CHANGES: 
 **------------------------------------------------------------------------**/

#include "vector.h"

/* ------------------------------------------------------------------
	 This function does the same as fill_gaps(), except that a grid for
	 interpolating points over gaps of bad values must be specified.

   RETURNS: the number of interpolated bad data
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int fill_gaps_on_grid(/* Array of data */
                      float *f,
                      /* Array of grid< points */
                      float *g, 
                      /*< Number of elements in arrays */
                      int n)
#else
int ()
#endif
{

	float *good_f, *good_g;
	int n_good, i, j;
	
	/* Count the number of good values in f */
	n_good = 0;
	for(i=0; i<n; i++)
		if(f[i] < ADJ_BADFLOAT)
			n_good++;

	/* Allocate memory for good values */
	good_f = (float *) calloc((size_t) n_good, sizeof(float));
	if(error_found(!good_f, "Cannot allocate memory for good data\n"))
		return(n);
	good_g = (float *) calloc((size_t) n_good, sizeof(float));
	if(error_found(!good_g, "Cannot allocate memory for good grid\n"))
		return(n);
	
	/* Copy good data */
	j = 0;
	for(i=0; i<n; i++){
		if(f[i] < ADJ_BADFLOAT){
			good_f[j] = f[i];
			good_g[j] = g[i];
			j++;
		}
	}
	if(error_found(j > n_good, "While making copy of good data\n"))
		return(n);

	/* Interpolate on grid */
	s_regridli(good_g, good_f, g, f, n_good, n); 
	
	/* Fill gaps at beginning and end of vector */
	fill_gaps(f, n);

	/* Release memory */
	free(good_f);
	free(good_g);

	return(n - n_good);
	
}

/* ------------------------------------------------------------------
	 This function checks if gaps should be filled or not. It looks to the
	 value of z below and above z0. If at least one of these is defined, it
	 checks if the maximal distance between these and z0 is smaller then the
	 specified gap value.

	 RETURNS: -1 on error, 0 if no filling should be performed and 1 filling can
            be carried out 
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int do_fill_gaps(/* Array of grid points */
                 float *z, 
                 /* Array of data at grid points */
                 float *f, 
                 /* Reference grid point */
                 float z0, 
                 /* Maximal gap distance allowed */
                 float max_gap, 
                 /* Number of elements in arrays */
                 int n)
#else
int do_fill_gaps()
#endif
{

	float *gz, *gf, gap;
	int i, j;
	int i0, i1, error, n_good;

	/* Initializations */
	i0 = -1;
	i1 = -1;
	error = 1;
	gz = (float *) NULL;
	gf = (float *) NULL;
	gap = MAXFLOAT;

	/* Get the number of good grid points */
	n_good = 0;
	for(i=0; i<n; i++)
		if((z[i] < ADJ_BADFLOAT) && (f[i] < ADJ_BADFLOAT))
			n_good++;
	
	/* Check if there is at least one good grid point */
	if(n_good<2)
		return(0);

	/* Make a copy of good grid points only */
	gz = (float *) calloc((size_t) n_good, sizeof(float));
	if(error_found(!gz, "Cannot allocate memory for good z\n"))
		goto on_error;
	gf = (float *) calloc((size_t) n_good, sizeof(float));
	if(error_found(!gf, "Cannot allocate memory for good f\n"))
		goto on_error;
	j = 0;
	for(i=0; i<n; i++){
		if((z[i] < ADJ_BADFLOAT) && (f[i] < ADJ_BADFLOAT)){
			gz[j] = z[i];
			gf[j] = f[i];
			j++;
		}
	}
	if(error_found(j>n_good, "Error while making copy of good grid points\n"))
		goto on_error;
		
	/* Find the index of z on both sides of z0 */
	i = 1;
	while((i<n_good) && (gz[i] < z0)) i++;
	if((i == 1) && (gz[0] >= z0))
		i1 = 0;
	else if((i == n_good) && (gz[n_good-1] < z0))
		i0 = n_good-1;
	else{
		i0 = i-1;
		i1 = i;
	}
	
	/* Get the maximal gap */
	if(i0 >= 0)
		gap = max_val(gap, abs_val(gz[i0]-z0));
	if(i1 >= 0)
		gap = max_val(gap, abs_val(gz[i1]-z0));
	
	error = 0;

on_error:

	if(gz) free(gz);
	if(gf) free(gf);

	/* Check error */
	if(error)
		return(-1);

	/* Check the gap distance */
	if(gap > max_gap)
		return(0);
	
	return(1);
	
}

/* ------------------------------------------------------------------
   Function s_regrid() expects vectors to be of type float. This
   function is an interface for other input vectors of type double.

	 RETURNS: return value from s_regrid(), or EOF if an error occured.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int s_regrid_d(/* Array of old grid points */
               double *zo, 
               /* Array of old data values */
               double *fo, 
               /* Array of new desired grid points */
               double *zn, 
               /* Array of new data values */
               double *fn, 
               /* Number of elements in old arrays */
               int no, 
               /* Number of elements in new arrays */
               int nn)
#else
int ()
#endif
{

	float *fzo, *ffo, *fzn, *ffn;
	int n_good, i, ierr;

	/* Allocate memory for floats */
	fzo = (float *) calloc((size_t) no, sizeof(float));
	ffo = (float *) calloc((size_t) no, sizeof(float));
	fzn = (float *) calloc((size_t) nn, sizeof(float));
	ffn = (float *) calloc((size_t) nn, sizeof(float));
	ierr = (!fzo || !ffo || !fzn || !ffn) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): float vectors"))
		return(EOF);
	
	/* Fill float vectors */
	for(i=0; i<no; i++){
		fzo[i] = ffo[i] = BADFLOAT;
		if(zo[i] < ADJ_BADFLOAT) fzo[i] = (float) zo[i];
		if(fo[i] < ADJ_BADFLOAT) ffo[i] = (float) fo[i];
	}
	for(i=0; i<nn; i++){
		fzn[i] = ffn[i] = BADFLOAT;
		if(zn[i] < ADJ_BADFLOAT) fzn[i] = (float) zn[i];
		if(fn[i] < ADJ_BADFLOAT) ffn[i] = (float) fn[i];
	}
	
	/* Regrid data */
	n_good = s_regrid(fzo, ffo, fzn, ffn, no, nn);

	/* Paste interpolated values into double vector */
	for(i=0; i<nn; i++){
		fn[i] = BADDOUBLE;
		if(ffn[i] < ADJ_BADFLOAT) fn[i] = (double) ffn[i];
	}
	
	/* Release memory */
	free((void *) fzo);
	free((void *) ffo);
	free((void *) fzn);
	free((void *) ffn);

	return(n_good);

}

