/**---------------------------------------------------------------------------
 ** 
 ** process_.c -- 
 ** 
 ** Author          : Pierre Jaccard
 ** Created On      : 1999/07/31 21:20:24
 ** Last Modified By: Pierre Jaccard
 ** Last Modified On: 1999/10/28 16:48:14
 ** Update Count    : 8
 ** Directory       : /pcdata1/jaccard/codas3c/gfi/src/apps/nbr/
 ** Version         : 0.0
 ** Status          : Unknown
 ** ---------------------------------------------------------------------- ** 
 ** DESCRIPTION: 
 ** 
 **    Undocumented.
 ** 
 ** ---------------------------------------------------------------------- ** 
 ** REVISIONS: 
 ** ---------------------------------------------------------------------- ** 
 ** CHANGES: 
 **------------------------------------------------------------------------**/

/* ##################################################################

	-Introduction:
	
	Utility functions fpr processing NB ADCP raw data. This file should be
	included in the main processing programme file.
	
	##################################################################
*/

#ifndef process__C_code_already_included
#define process__C_code_already_included

#include "process_.h"

/* ##################################################################
	 -Data Scaling:

	 The following describes how to convert data in NB ADCP raw format into
	 their physical format. The procedure to do this needs several lists, and is
	 described below:

	 +The function [[convert_nbr]] control the complete conversion process.

	 +The [[NBR_SCALING_TYPE]] list contains information about the scaling
	  parameters for each variable that is supposed to be converted.The
		different variables are identified by a conversion code (an
		integer). These codes are defined in [[conv_nbr.h]]. The list is
		initialized in the same file. An additionnal field [[name]] is present
		only for a better clarity of the list initialization, as well as may be
		output of the list content. 

		Scaling parameters with a [[BADDOUBLE]] value are assumed to be unknown at
		initialization. They must therefore be defined while the programme is
		running.

	 +The data to be converted is defined in an [[UNISTAT_LIST_TYPE]] list. Each
	  unistat element of this list will contain the converted data as
		[[double]]. This list is initialized in the main programme.

	 +In order to have a shorter code, a [[NBR_CONVERSION_TYPE]] list is
	  defined. Each element has a conversion code identifying one of the
	  variables defined in the [[NBR_SCALING_TYPE]] list. The data for each
	  variable can be reached through a pointer to that variable, and the type
	  associated to it. The code for the type must one of those defined in
	  [[dbext.h]]. The pointer must point to the NB ADCP raw data. 

		Hence, given a conversion code, the NB ADCP data can be accessed through
		this pointer. Using the same conversion code, the scaling parameters can
		be retrieved from the [[NBR_SCALING_TYPE]] list. Then, the NB ADCP raw
		data can be converted.

	 +In order to know where to look for the different NB ADCP raw data (this
	  depends on the options in the control file), a
		[[UNISTAT_VS_CONVERSION_TYPE]] list is defined. Each variable from the
		[[UNISTAT_LIST_TYPE]] has an entry in this list with the same code that
		identifies it. In addition, a conversion code which corresponds to an
		entry in the [[NBR_CONVERSION_TYPE]] and [[NBR_SCALING_TYPE]] lists is
		associated to each variable. A conversion code of [[BADINT]] in this list
		assumes that conversion is undefined at initialization. It must therefore
		be defined at startup, when the options from the control file are known.

	 +IN SHORT: 

	     +Each converted data is defined in the [[UNISTAT_LIST_TYPE]]. 
			 +The code identifying the variable in this list as an assciated
			  conversion code in the [[UNISTAT_VS_CONVERSION_TYPE]].
			 +This conversion code has associated elements in the
			  [[NBR_SCALING_TYPE]] and [[NBR_CONVERSION_TYPE]] lists. 
			 +The corresponding element in the [[NBR_CONVERSION_TYPE]] list defines
			  how to get access to the NB ADCP raw data.
			 +The corresponding elsemnt in the [[NBR_SCALING_TYPE]] list defines how
			  the scaling from NB ADCP raw data should be carried out.

	 WARNING: the structure containing all necessary control settings for
	          converting data must be initialized in the main programme (see the
						[[NBR_CNV_CONTROL_TYPE]] structure]] in [[conv_nbr.h]]).
   ################################################################## */

/* ##################################################################
	 -Beam Vectors and Rotations:
	 ################################################################## */

/* ------------------------------------------------------------------
	 --Beam Vector:

	 This function return a pointer to the beam vectors, in ADCP coordinate
	 system.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int fill_beam_vectors(double *v, PROCESS_CONTROL_TYPE *cnt) 
#else
int fill_beam_vectors(v, cnt)
#endif
{

	char func[]="fill_beam_vectors";

	BEAM_VECTOR_SIGN_TYPE *list;
	UBYTE and, xor;
	int done, i, j, k;
	double vec[3];
	double angle;

	/* Convert angle into radians */
	angle = (cnt->ba) * M_PI/180.0;

	/* Build the base beam vector */
	vec[0] = vec[1] = sin(angle);
	vec[2] = cos(angle);
	
	/* Get the sign factor for each beam, according to the adcp configuration */
	list = cnt->bvsl;
	done = 0;
	i = 0;
	while((i<4) && !done){
		and = (list+i)->and_mask;
		xor = (list+i)->xor_mask;
		done = !( (*(cnt->adcp_config) & and) ^ xor);
		i++;
	}
	i--;
	/* Check that something has been found */
	if(!done){
		error_msg(func, "No parameters defined for beam vectors");
		return(CONVERT_ERROR);
	}
		
	/* Apply the sign to the beam elements */
	for(j=0; j<4; j++)
		for(k=0; k<3; k++)
			*(v + j*3 + k) = vec[k] * ((double) *( ((list+i)->sign) + 3*j + k));
	
	return(0);
}

/* ------------------------------------------------------------------
	 --Get Rotation Matrix:

	 Given angles for heading, pitch and roll expressed in degrees, this
	 function fills the matrix that transforms vectors in earth coordinates.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int fill_rotation_matrix(double *mat, double h, double p, double r)
#else
int fill_rotation_matrix(mat, h, p, r)
#endif
{
		
	double R[3][3];

	/* Convert angles to radians */
	h *= M_PI/180.0;
	p *= M_PI/180.0;
	r *= M_PI/180.0;
	
	/* Build the rotation matrix */
	R[0][0] = cos(h)*cos(r)+sin(h)*sin(p)*sin(r);
	R[0][1] = cos(p)*sin(h);
	R[0][2] = cos(h)*sin(r) - cos(r)*sin(h)*sin(p);
	R[1][0] = cos(h)*sin(p)*sin(r) - cos(r)*sin(h);
	R[1][1] = cos(h)*cos(p);
	R[1][2] = -cos(h)*cos(r)*sin(p) - sin(h)*sin(r);
	R[2][0] = -cos(p)*sin(r);
	R[2][1] = sin(p);
	R[2][2] = cos(p)*cos(r);

	/* Copy data into the matrix */
	memcpy((void *) mat, (void *) R, 9*sizeof(double));
	
	return(0);

}
/* ------------------------------------------------------------------
	 --Matrix Multiplication:

	 mat3 = mat2 * mat1

	 mat1: n by p
	 mat2: m by n
	 mat3: m by p
	 
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int matrix_mult_d(double *mat3, double *mat2, double *mat1, 
									int m, int n, int p)
#else
int matrix_mult_d()
#endif
{

	int i, j, k;
	double mat1_jk, mat2_ij, *mat3_ik;
	double sum;

	/* mat3(i,k) = mat2(i,j) * mat1(j,k) */
	for(i=0; i<m; i++){
		for(k=0; k<p; k++){
			mat3_ik = mat3 + i*p + k;
			sum = 0.0;
			for(j=0; j<n; j++){
				mat1_jk = *(mat1 + j*p + k);
				mat2_ij = *(mat2 + i*n + j);
				if( (mat1_jk < ADJ_BADFLOAT) && (mat2_ij < ADJ_BADFLOAT)
						&& (sum < ADJ_BADFLOAT))
					sum += mat2_ij * mat1_jk;
				else
					sum = BADDOUBLE;
			}
			*(mat3_ik) = sum;
		}
	}

	return(0);

}

/* ------------------------------------------------------------------
	 --Fill Bin Mapping Index Table:

	 Given [[v]], a matrix whose columns are the beam_vectors in the ADCP
	 coordinate system, [[mat]] being the transformation matrix into earth
	 coordinate system, and finally [[angle]] being the beam angle, this
	 function returns a bin mapping index table for each bin cell.

	 The content of bin_map[i][j] is the index of the bin in beam [j] that is at
	 a depth corresponding to [i], the latter being the depth relative to the
	 earth vertical axis.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int fill_bin_map_table(double *v, double *mat, double angle, UBYTE *bin_map)
#else
int fill_bin_map_table(v, mat, angle, bin_map)
#endif
{

	char func[]="fill_bin_mapping_index";

	double v_adcp[3], v_earth[3];
	double r;
	int i, j, k;

	/* Convert beam angle into degrees */
	angle *= M_PI/180.0;

	/* Loop through the number of beams */
	for(i=0; i<4; i++){
		/* Get the beam vector */
		memcpy((void *) v_adcp, (void *) (v + i*3), 3*sizeof(double));
		/* Get the cell size ratio */
		matrix_mult_d(v_earth, mat, v_adcp, 3, 3, 1);
		if(abs_val(v_earth[2]) > MINFLOAT) 
			r = abs_val(cos(angle)/v_earth[2]);
		else{
			error_msg(func, "Inifinite cell size ratio");
			return(CONVERT_ERROR);
		}
		/* Loop through the bins */
		for(j=0; j<MAX_BINS; j++){
			/* Calculate the new index */
			k = round(j*r);
			/* Security checks */
			if(k < 0) k=0;
      /* 2000/12/05 PJ
         The max value was set to MAX_BINS instead of (MAX_BINS-1). Since the
         content of bin_map is used as indices in function apply_bin_mapping,
         there was a bound violation.
      */
			if(k >= MAX_BINS) k=MAX_BINS-1;

			/* Save the value */
			*(bin_map + j*4 + i)  = (UBYTE) k;
		}
	}

	return(0);

}

#if PROTOTYPE_ALLOWED
int apply_bin_mapping(UNISTAT_LIST_TYPE *list, int code, UBYTE *bin_map)
#else
int apply_bin_mapping(list, code, bin_map)
#endif
{

	char func[]="apply_bin_map";

	UNISTAT_LIST_TYPE *uni;

	double data[MAX_BINS][4];
	UBYTE  bin;
	int i, j, n;

	/* Get correct unistat variable */
	uni = (UNISTAT_LIST_TYPE *) get_unistat_element(list, code);
	if(!uni){
		error_msg(func, "No statistical data available for this code");
		return(1);
	}

	/* Loop through the bins */
	for(i=0; i<(uni->size); i+=4){
		
		/* Calculate bin number */
		n = i/4;

		/* Get bin mapping values and data */
		for(j=0; j<4; j++){
			bin = *(bin_map + n*4 + j);
			data[n][j] = *(uni->data + bin*4 + j);
		}
	}

	/* Put back data into the unistat list */
	for(i=0; i<(uni->size); i++)
		*(uni->data + i) = *(((double *) data) + i);

	return(0);

}

/* ##################################################################
	 -Beam Checking:

	 The following functions are used to determine the beam solution function,
	 depending on the beam data quality.
	 ################################################################## */

/* ------------------------------------------------------------------
	 --Select Good Beams:

	 This function looks at the beam values in the specified [[bins]] and
	 returns an [[ULONG]] specifiying which beams may be considered as good.

	 Several methods may be used to determine good beam values. These are
	 selected by setting the corresponding bits in the [[meth]] parameter. If
	 several methods are selected, then the composition of the results is
	 carried out according to the [[op]] operator (see the description below).

	 The different methods are applied in increasing bits number order.

	 The following methods are available:
	 
	    BIT0 : UNORDERED BAD DATA

				    +set only beams that are not bad
						+do not order beams

	    BIT1 : ORDERED SMALL DATA

				    +calculate mean value for the four beams,
						+beams that have a value larger than (<mean> - <sdev>) are
						 considered good
						+order beams according to their value

   The returned code should be interpreted as follows:

	    BIT0  to BIT3  : Beam is good if its corresponding bit is set.

			BIT4  to BIT5
			BIT6  to BIT7
			BIT8  to BIT9
			BIT10 to BIT11 : Order of goodness for each beam. The two bits
			                 assigned for each beam form a value inside [0,3],
											 which gives its order of goodness relative to the other
											 beams.

	 If several methods should be carried out, the operator is applied in the
	 following way:

	    BIT0  to BIT3  : 

			    0 : all tests must be successful 
				  1 : at least one test must be successful

			BIT4  to BIT11 : 

			    0 : the smallest order is kept
				  1 : the largest order is kept

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

#if PROTOTYPE_ALLOWED
ULONG check_beams(double *data, ULONG meth, BYTE op)
#else
ULONG check_beams(data, meth, op)
#endif
{

	char func[]="check_beams";

	UNISTAT_TYPE stat;
	LONG   ul_data[4];
	int index[4];

	int i, j;
	int order_flag;
	ULONG check_bits = BIT31;
	ULONG bits;
	ULONG good_mask, order_mask, min_mask, max_mask;

	/* Initializations */
	allocate_unistat(&stat, 1, "beam_data", U_ALL);
	zero_unistat(&stat);
	
	/* Apply each method successively */
	for(i=0; i<32; i++){
		/* Check if the corresponding method is set */
		if(meth & ((ULONG) to_bit(i))){
			/* Initialize the bits */
			bits = 0;
			/* Select the corresponding method */
			switch(i){
				
			case 0:				
				/* Test data values */
				for(j=0; j<4; j++){
					if(data[j] < ADJ_BADFLOAT)
						bits |= (ULONG) to_bit(j);
				}

				/* set the order flag */
				order_flag = 0;
				break;

			case 1:  
				/* Calculate the statistics */
				for(j=0; j<4; j++)
					update_unistat_d(&stat, data+j, 1);
				calculate_unistat(&stat);

				/* Select good beams */
				if(*(stat.sum) < ADJ_BADFLOAT){
					for(j=0; j<4; j++)
						if((data[j] < ADJ_BADFLOAT) 
							 && (data[j] >= (*(stat.sum) - sqrt(*(stat.sumsq)))))
							bits |= (ULONG) to_bit(j);
				}
				else
					bits = 0;

				/* Set the order flag */
				order_flag = 1;
				break;

			default:
				/* Nothing to do */
				continue;
				break;
			} /* switch(i) ; i=method */
			
			if(order_flag){
				/* Order the beams: convert first the array of data into an array of
					 ULONG, so that the already existing sort function can be used */
				double_to_long(ul_data, data, 4);
				sort_dir_index(4, index, (CHAR *) ul_data, sizeof(LONG));
				/* Fill the code for ordered beams */
				for(j=0; j<4; j++)
					bits |= ( ((ULONG) j) << (4 + 2*index[j]) );
			}
			else{
				switch(op){
				case 0:
					bits |= 0xFF0;
					break;
				}
			}

			/* Apply the composition */
			if(check_bits & BIT31){
				/* First pass */
				check_bits = bits;
				check_bits &= ~BIT31;
			}
			else{
				for(j=0; j<4; j++){
					/* Build the mask to get the good beam bit */
					good_mask = (ULONG) to_bit(j);
					/* Build the mask to get the order bits */
					order_mask = ((ULONG) 0x3) << (4 + j*2);
					switch(op){
					case 0:
						/* Good beam flag */
						if((bits & good_mask) && (check_bits & good_mask))
							check_bits |= good_mask;
						else
							check_bits &= ~good_mask;
						/* Order flag */
						min_mask = min_val((check_bits & order_mask), (bits & order_mask));
						check_bits |= ((check_bits & ~order_mask) | min_mask);
						break;
					case 1:
						if((bits & good_mask) || (check_bits & good_mask))
							check_bits |= good_mask;
						else
							check_bits &= ~good_mask;
						/* Order flag */
						max_mask = max_val((check_bits & order_mask), (bits & order_mask));
						check_bits |= ((check_bits & ~order_mask) | max_mask);
						break;
					default:
						warning_msg(func, "No composition defined for this operator");
						break;
					}
				}
			} /* check next beam */
		} /* if meth specified */
		/*Continue with the next possible method */
	}

	/* Release memory */
	free_unistat(&stat);

	return(check_bits);
		
}

/* ##################################################################
	 -Beam Solutions:

	 The following functions perform the different solution for transforming
	 velocity into the earth coordinate system. Their return value is either the
	 beam solution number, if successful, zero if failed, or negative if the
	 velocities component have been set to the bad value.
	 ################################################################## */

/* ------------------------------------------------------------------
	 --Four Beam Solution:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int four_beam_solution(double *vb, double *bv, ULONG checks, double *v)
#else
int four_beam_solution(vb, bv, checks, v)
#endif
{

	double bx, by, bz;
	
	/*Calculate the velocity */
	bx = *(bv + 0*3 + 0);
	by = *(bv + 2*3 + 1);
	bz = *(bv + 0*3 + 2);
	v[0] = (vb[0] - vb[1])/(2.0 * bx);
	v[1] = (vb[2] - vb[3])/(2.0 * by);
	v[2] = (vb[0] + vb[1] + vb[2] + vb[3])/(4.0 * bz);
	v[3] = (vb[0] + vb[1] - vb[2] - vb[3])/(4.0 * bz);

	return(4);

}
 
/* ------------------------------------------------------------------
	 --Three Beam Solution:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int three_beam_solution(double *vb, double *bv, ULONG checks, double *v)
#else
int three_beam_solution(vb, bv, checks, v)
#endif
{

	int index[4];
	int ix, iy;
	double bx, by, bz;

	/* Build the index table of selected good beams */
 	ix = 0;
	iy = 2;
	if(checks & BIT0) index[ix++] = 0;
	if(checks & BIT1) index[ix++] = 1;
	if(checks & BIT2) index[iy++] = 2;
	if(checks & BIT3) index[iy++] = 3;
		
	/* Calculate the velocity */
	if(ix == 2){
		/* Beam #1, #2 and (#3 or #4) */
		bx = *(bv + index[0]*3 + 0);
		by = *(bv + index[2]*3 + 1);
		bz = *(bv + index[0]*3 + 2);
		v[0] = (vb[0] - vb[1])/(2.0 * bx);
		v[2] = (vb[0] + vb[1])/(2.0 * bz);
		v[1] = (vb[index[2]] - (vb[0] + vb[1])/2.0)/by;
		v[3] = BADDOUBLE;
	}
	else if(iy == 4){
		/* Beam (#1 or #2) and #3, #4 */
		bx = *(bv + index[0]*3 + 0);
		by = *(bv + index[2]*3 + 1);
		bz = *(bv + index[2]*3 + 2);
		v[1] = (vb[2] - vb[3])/(2.0 * by);
		v[2] = (vb[2] + vb[3])/(2.0 * bz);
		v[0] = (vb[index[0]] - (vb[2] + vb[3])/2.0)/bx;
		v[3] = BADDOUBLE;
	}
	else
		return(0);

	return(3);

}

/* ------------------------------------------------------------------
	 --Two Beam Solution:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int two_beam_solution(double *vb, double *bv, ULONG checks, double *v)
#else
int two_beam_solution(vb, bv, checks, v)
#endif
{

	double bx, by, bz;
	int ix, iy;
	int index[4];

	/* Build the index table of selected good beams */
 	ix = 0;
	iy = 2;
	if(checks & BIT0) index[ix++] = 0;
	if(checks & BIT1) index[ix++] = 1;
	if(checks & BIT2) index[iy++] = 2;
	if(checks & BIT3) index[iy++] = 3;
	
	/* Calculate velocity */
	if(ix == 2){
		/* Beam #1 and #2 */
		bx = *(bv + index[0]*3 + 0);
		bz = *(bv + index[0]*3 + 2);
		v[0] = (vb[0] - vb[1])/(2.0 * bx);
		v[2] = (vb[0] + vb[1])/(2.0 * bz);
		v[1] = v[3] = BADDOUBLE;
	}
	else if(iy == 4){
		by = *(bv + index[2]*3 + 1);
		bz = *(bv + index[2]*3 + 2);
		v[1] = (vb[2] - vb[3])/(2.0 * by);
		v[2] = (vb[2] + vb[3])/(2.0 * bz);
		v[0] = v[3] = BADDOUBLE;
	}
	else
		return(0);

	return(2);

}

/* ------------------------------------------------------------------
	 --Bad Beam Solution:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int bad_beam_solution(double *v)
#else
int bad_beam_solution(v)
#endif
{

	int i;

	for(i=0; i<4; i++)
		v[i] = BADDOUBLE;

	return(-1);

}

/* ------------------------------------------------------------------
	 -Get PRH Data
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_prh(UNISTAT_LIST_TYPE *list, double *p, double *r, double *h)
#else
int get_prh()
#endif
{

	UNISTAT_LIST_TYPE *uni;
	
	uni = (UNISTAT_LIST_TYPE *) 
		get_unistat_element(list, HEADING_CODE);
	if(uni) *h = *(uni->data);
	else    return(CONVERT_ERROR);

	uni = (UNISTAT_LIST_TYPE *) 
		get_unistat_element(list, PITCH_CODE);
	if(uni) *p = *(uni->data);
	else    return(CONVERT_ERROR);

	uni = (UNISTAT_LIST_TYPE *) 
		get_unistat_element(list, ROLL_CODE);
	if(uni) *r = *(uni->data);
	else    return(CONVERT_ERROR);

	return(0);

}

/* ------------------------------------------------------------------
	 -Transform Velocity:

	 This function controls how velocity should be transformed.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int transform_one_velocity(double *vb, double *bv, ULONG checks, 
													 int sol, int force, double *v)
#else
int transform_one_velocity()
#endif
{

	char func[]="transform_one_velocity";

	int ngood, i;
	
	/* Check if the number of good beams corresponds to the solution number
		 specified */
	ngood = 0;
	for(i=0; i<4; i++)
		if(checks & ((ULONG) to_bit(i)))
			ngood++;
	/* In case of 2 beam solution, only paired beams can be used */
	if(sol == 2)
		if( ((checks & 0x11) ^ 0x3) && ((checks & 0xC) ^ 0xC) )
			ngood = 1;
	
	if(ngood < sol){
		if(force){
			i=bad_beam_solution(v);
			return(i);
		}
		else
			return(0);
	}
	
	/* Select among the possible solutions */
	switch(sol){
	case 4:
		i=four_beam_solution(vb, bv, checks, v);
		break;
	case 3:
		i=three_beam_solution(vb, bv, checks, v);
		break;
	case 2:
		i=two_beam_solution(vb, bv, checks, v);
		break;
	case 1:
		i=bad_beam_solution(v);
		break;
	default:
		warning_msg(func, "No solution defined for this code");
		break;
	}
	
	return(i);

}
	
/* ------------------------------------------------------------------
	 -Transform All Velocity:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int transform_all_velocity(LOAD_NBR_TYPE *rec, double *adcp_to_earth,
													 double *bv, int code) 
#else
int transform_all_velocity()
#endif
{

	UNISTAT_LIST_TYPE *vel, *bsol;
	FLAG_LIST_TYPE    *fsol;

	int i, j, n, done, flag;
	int last_good = 0;

	ULONG checks;
	double vb[4], v[MAX_BINS][4];
	double v0[3], v1[3];

	/* Get pointer to unistat velocity and pcg */
	vel = (UNISTAT_LIST_TYPE *) 
		get_unistat_element(rec->stat_list, code);
	
	/* Loop through the data */
	for(i=0; i<(vel->size); i+=4){
		
		/* Calculate bin number */
		n = i/4;

		/* Retrieve beam velocities */ 
		for(j=0; j<4; j++)
			vb[j] = *( vel->data + 4*n + j);
		
		/* Check the beams */
		checks = check_beams(vb, BIT0, 0);
		
		/* Try successively velocity transformations, starting with the highest
			 number of beams */
		done = 0;
		j = (rec->cnt->max_bs > 4) ? 4 : rec->cnt->max_bs;
		while((j > rec->cnt->min_bs) && (j>1) && !done){
			done = transform_one_velocity(vb, (double *) bv, checks, j, 0, v[n]);
			j--;
		}
		/* Force a solution if necessary */
		if(!done)
			done = transform_one_velocity(vb, (double *) bv, checks, j, 1, v[n]);
		
		/* Set the beam solution flag */
		if(done < 0) done = 0;
		set_flag(rec->flag_list, BEAM_SOLUTION_FLAG, n, done);

		/* Keep the last good bin */
		if(done > 0) last_good = n;
						
		/* Rotate the velocity data into earth coordinate system */
		for(j=0; j<3; j++)
			v0[j] = v[n][j];
		matrix_mult_d(v1, adcp_to_earth, v0, 3, 3, 1); 
		for(j=0; j<3; j++)
			v[n][j] = v1[j];

	}

	/* Paste the value of transformed velocity back into the velocity data
		 array */
	for(i=0; i<(vel->size); i+=4){
		n = i/4;
		for(j=0; j<4; j++)
			*(vel->data + n*4 + j) = v[n][j];
	}
	last_good++;
	last_good = (last_good > MAX_BINS) ? MAX_BINS : last_good;

	/* 98/08/28 PJ
		 Calculate the percentage of beam solutions over all bins. On output,
		 indexes have the following meaning:

		 0: % bad beams solutions
		 1: % 2 beams solutions
		 2: % 3 beams solutions
		 3: % 4 beams solutions 
		 
		 This is then averaged over the averaging time, and should give and idea
		 of solutions encountered during each procesed ensemble. However, over
		 shallow depths, the percentage of bad beam solutions will dominate */

  bsol = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, -code);
  fsol = (FLAG_LIST_TYPE *) 
		get_flag_element(rec->flag_list, BEAM_SOLUTION_FLAG);
	for(i=0; i<(bsol->size); i++)
			*(bsol->data + i) = 0.0;
	for(i=0; i<last_good; i++){
		flag = *(fsol->flag + i);
		flag = (flag > 0) ? flag-1 : 0;
		*(bsol->data + flag) += 1.0;
	}
	for(i=0; i<(bsol->size); i++)
		*(bsol->data + i) = *(bsol->data + i)/last_good;
	if(code == VELOCITY_CODE){
		bsol = (UNISTAT_LIST_TYPE *)
			get_unistat_element(rec->stat_list, LAST_GOOD_BIN_CODE);
		*(bsol->data) = (DOUBLE) last_good;
	}

	return(0);

}

/* ------------------------------------------------------------------
	 -Process Echo Intensity:
	 
	 This function tests for bad echo intensities and put these to the bad
	 value. It should be called before processing other variables, since its
	 purpose is also to check malfunctions of single beams.

	 In addition to check for bad beam intensities, it calculates the mean beam
	 echo intensity.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int process_echo(LOAD_NBR_TYPE *rec)
#else
int process_echo(rec)
#endif
{

	char func[]="process_echo";

	UNISTAT_LIST_TYPE *echo, *mn_echo;
	double data[4];
	double sum;
	ULONG checks;
	int i, n, j, m;

	/* Get unistat variable for raw echo */
  echo = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, RAW_ECHO_CODE);
	if(!echo){
		error_msg(func, "No statistics available for [[RAW_ECHO_CODE]]");
		return(1);
	}
  mn_echo = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, RAW_MEAN_AMP_CODE);
	if(!mn_echo){
		error_msg(func, "No statistics available for [[RAW_MEAN_AMPLITUDE]]");
		return(1);
	}

	/* Loop through the bins */
	for(i=0; i<(echo->size); i+=4){
		
		/* Get bin number */
		n = i/4;
		
		/* Get echo data */
		for(j=0; j<4; j++)
			data[j] = *(echo->data + n*4 + j);
			
		/* Check the value of each beam */
		checks = check_beams(data, BIT1, 0);
		
		/* Set selected bad data */
		for(j=0; j<4; j++){
			if(!(checks & (ULONG) to_bit(j)))
				*(echo->data + n*4 + j) = BADDOUBLE;
		}

		/* Calculate mean beam value */
		m = 0;
		sum = 0.0;
		for(j=0; j<4; j++){
			if(*(echo->data + n*4 + j) < ADJ_BADFLOAT){
				m++;
				sum += *(echo->data + n*4 + j);
			}
		}
		if(m > 0)
			*(mn_echo->data + n) = sum/m;
		else
			*(mn_echo->data + n) = BADDOUBLE;

	}

	return(0);

}

#if PROTOTYPE_ALLOWED
int process_pcg(LOAD_NBR_TYPE *rec)
#else
int process_pcg(rec)
#endif
{

	char func[]="process_pcg";

	NBR_ENSEMBLE_TYPE *e;
	NBR_LEADER_TYPE *l;
	UNISTAT_LIST_TYPE *pcg, *mn_pcg, *echo;
	double sum;
	int i, j, k, n;

	/* Assign pointers */
	e = &(rec->ens);
	l = &(rec->ens.leader);

	/* Get unistat variables for pcg and mean pcg */
	pcg = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, PCG_CODE);
	if(!pcg){
		error_msg(func, "No statistics available for [[PCG]]");
		return(1);
	}
	mn_pcg = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, MEAN_PCG_CODE);
	if(!mn_pcg){
		error_msg(func, "No statistics available for [[MEAN_PCG]]");
		return(1);
	}
	echo = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, RAW_ECHO_CODE);
	if(!echo){
		error_msg(func, "No statistics available for [[RAW_ECHO]]");
		return(1);
	}

	/* 98/08/25 PJ
		 If percent good has been recorded, just update data with these values. If
		 not, one could count the good velocities, however, it is propbably better
		 to look at echo intensisties */
	if(e->presence_bits & PERCENT_GOOD_PRESENT){
		/* Loop through the bins */
		for(i=0; i<(pcg->size); i+=4){
		
			/* Calculate bin number */
			n = i/4;

			/* Calculate mean beam value of pcg */
      /* 1999/09/21 Pierre Jaccard
         Must test the value of each pcg before calculating the its mean
         value. */
			sum = 0.0;
      k = 0;
			for(j=0; j<4; j++){
        if( *(pcg->data + 4*n + j) < ADJ_BADFLOAT ){
          sum += *(pcg->data + 4*n + j);
          k++;
        }
      }
      if(k > 0)
        *(mn_pcg->data + n) = sum/k;
		}
	}
	else{
		/* Echo intensities have been processed in `process_echo()' */
		for(i=0; i<(echo->size); i+=4){
		
			/* Get bin number */
			n = i/4;
		
			/* Create pcg data */
			sum = 0.0;
			for(j=0; j<4; j++){
				if(*(echo->data + n*4 + j) < ADJ_BADFLOAT)
					sum += 1.0;
			}
			*(mn_pcg->data + n) = (sum/4)*100.0;
		}
	}

	/* Update statistics of these variables to take into account the number of
		 pings averaged inside the ADCP itself */
	for(j=0; j<(int)(l->pings_per_ensemble - 1); j++){
		update_unistat_list(pcg, PCG_CODE);
		update_unistat_list(mn_pcg, MEAN_PCG_CODE);
	}

	return(0);

}

#if PROTOTYPE_ALLOWED
int process_sw(LOAD_NBR_TYPE *rec)
#else
int process_sw(rec)
#endif
{

	char func[]="process_sw";

	NBR_ENSEMBLE_TYPE *e;
	NBR_LEADER_TYPE *l;
	UNISTAT_LIST_TYPE *sw, *mn_sw, *echo;
	ULONG checks;
	double sum, data[4];
	int i, j, n, npts;

	/* Assign pointers */
	e = &(rec->ens);
	l = &(rec->ens.leader);

	/* Get unistat variables for sw and mean sw */
	sw = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, SPEC_WIDTH_CODE);
	if(!sw){
		error_msg(func, "No statistics available for [[SPECTRAL_WIDTH]]");
		return(1);
	}
	mn_sw = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, MEAN_SPEC_WIDTH_CODE);
	if(!mn_sw){
		error_msg(func, "No statistics available for [[MEAN_SPEC_WIDTH]]");
		return(1);
	}
	echo = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, RAW_ECHO_CODE);
	if(!echo){
		error_msg(func, "No statistics available for [[RAW_ECHO]]");
		return(1);
	}
	
	if(e->presence_bits & SPECTRAL_WIDTH_PRESENT){
		/* Loop through the bins */
		for(i=0; i<(sw->size); i+=4){
		
			/* Calculate bin number */
			n = i/4;

			/* Get echo data */
			for(j=0; j<4; j++)
				data[j] = *(echo->data + n*4 + j);

			/* Check the value of each beam */
			checks = check_beams(data, BIT1, 0);
		
			/* Set selected bad data */
			for(j=0; j<4; j++){
				if(!(checks & (ULONG) to_bit(j)))
					*(sw->data + n*4 + j) = BADDOUBLE;
				data[j] = *(sw->data + n*4 + j);
			}
			
			/* Calculate mean beam value of pcg */
			sum = 0.0;
			npts = 0;
			for(j=0; j<4; j++){
				if(data[j]  < ADJ_BADFLOAT){
					sum += data[j];
					npts++;
				}
			}
			if(npts > 0)
				*(mn_sw->data + n) = sum/npts;
			else
				*(mn_sw->data + n) = BADDOUBLE;
				
		}
	}
	else{
		/* Set bad data everywhere */
		for(i=0; i<(sw->size); i+=4){
			
			/* Get bin number */
			n = i/4;
		
			/* Set data to bad */
			*(mn_sw->data + n) = BADDOUBLE;
			for(j=0; j<4; j++)
				*(sw->data + n*4 + j) = BADDOUBLE;
		}
	}

	return(0);

}

#if PROTOTYPE_ALLOWED
int process_bt(LOAD_NBR_TYPE *rec)
#else
int process_bt(rec)
#endif
{

	UNISTAT_LIST_TYPE *vel, *rng;
	int i;
	ULONG checks;

	/* Get pointers to bt range and velocities */
	vel = (UNISTAT_LIST_TYPE *) 
		get_unistat_element(rec->stat_list, BT_VELOCITY_CODE);
	rng = (UNISTAT_LIST_TYPE *)
		get_unistat_element(rec->stat_list, BT_RANGE_CODE);

	/* Check validity of velocity */
	checks = check_beams(vel->data, BIT0, 0);

	/* Process bt range values */
	for(i=0; i<4; i++)
		if(!(checks & (ULONG) to_bit(i)))
			*(rng->data + i) = BADDOUBLE;
			
	return(0);

}

/* ------------------------------------------------------------------
	 -Process Data:

	 This is the main processing function.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int process_data(LOAD_NBR_TYPE *rec, double *adcp_to_ship)
#else
int process_data()
#endif
{

	char func[]="process_data";
	
	int *bin_map_var;
	int error;

	double p, r, h;
	double ship_to_earth[3][3], adcp_to_earth[3][3];

	double bv[4][3];
	UBYTE  bin_map[MAX_BINS][4];
	
	/* Conversion into physical format */
	error = nbr_cnv_all_data(rec->stat_list, (char *) rec, rec->cnv);
	if(error){
		error_msg(func, "While converting raw data");
		return(error);
	}

	/* Calculate the rotation matrix for adcp-to-earth transformation */
	if( get_prh(rec->stat_list, &p, &r, &h) )
		return(CONVERT_ERROR);

	/* Calculate full rotation matrix */
	fill_rotation_matrix((double *) ship_to_earth, h, p, r);
	matrix_mult_d((double *) adcp_to_earth, (double *) ship_to_earth, 
								(double *) adcp_to_ship, 3, 3, 3);

	/* Get beam vectors */
	fill_beam_vectors((double *) bv, rec->cnt); 

	/* Apply bin mapping */
	if(option_bits & BIN_MAPPING_BIT){
		fill_bin_map_table((double *) bv, (double *) adcp_to_earth, rec->cnt->ba, 
											 (UBYTE *) bin_map);
		bin_map_var = rec->cnt->bml;
		while(*bin_map_var != BADINT){
			apply_bin_mapping(rec->stat_list, *bin_map_var, (UBYTE *) bin_map);
			bin_map_var++;
		}
	}

	/* Check for bad beams */
	process_echo(rec);

	/* Check if velocity should be transformed */
	if( !(*(rec->cnt->adcp_config) & 0x2)){ 
	
		if( transform_all_velocity(rec, (double *) adcp_to_earth, (double *) bv,
															 VELOCITY_CODE) ){
			error_msg(func, "While transforming [[VELOCITY]]");
			return(CONVERT_ERROR);
		}

		if( transform_all_velocity(rec, (double *) adcp_to_earth, (double *) bv, 
															 BT_VELOCITY_CODE) ){
			error_msg(func, "While transforming [[BT_VELOCITY]]");
			return(CONVERT_ERROR);
		}
	}

	/* Process other variables */
	process_bt(rec);
	process_pcg(rec);
	process_sw(rec);

	return(0);
}

#endif /* process__C_code_already_included */
