/**---------------------------------------------------------------------------
 ** 
 ** tr_ctd.c -- 
 ** 
 ** Author          : Pierre Jaccard
 ** Created On      : 1999/08/04 13:52:42
 ** Last Modified By: Pierre Jaccard
 ** Last Modified On: 1999/08/04 14:06:53
 ** Update Count    : 5
 ** Directory       : /home/pego/pcd1/codas3c/gfi/src/apps/util/
 ** Version         : 0.0
 ** Status          : Unknown
 ** ---------------------------------------------------------------------- ** 
 ** DESCRIPTION: 
 ** 
 **    Undocumented.
 ** 
 ** ---------------------------------------------------------------------- ** 
 ** REVISIONS: 
 ** ---------------------------------------------------------------------- ** 
 ** CHANGES: 
 **------------------------------------------------------------------------**/

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

	-Introduction:

	[[tr_ctd.c]] utility for interpolating CTD data to transducer depth at ADCP
	times or positions.

	1998/11/27 PJ
	           Originally, this programme used to read simultaneously data from
	           ADCP and CTD databases. However, a segmentation fault occured
	           when reading ADCP data from a CODAS database generated with
	           loadtnbp (unexplained bug). Therefore I rewrote it for reading
	           ADCP positions from the output of getnav. The only loss by doing
	           so is that the transducer depth must be constant during all
	           specified time ranges (can this be a problem for someone? - use
	           xtract to list the transducer depth along the selected axis and
	           modify the programme for reading transducer depth from an
	           external file).

						 Previous name was [[lst_ctd]]
	
	##################################################################

*/

#include "tr_ctd.h"

/* ------------------------------------------------------------------
	 -Initializations of Database Extraction
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int init_xvar(EXTRACT_VARIABLE_TYPE *xvar)
#else
int init()
#endif
{

	int i;

	/* Initialze xtraction structures */
	for(i=0; i<MAX_VARS; i++)
		init_xtract_prf_var(xvar+i);

	/* Initialize the name of xtraction vars */
	strcpy(xvar[T_INDEX].name  , "TIME"       );
	strcpy(xvar[LON_INDEX].name, "LONGITUDE"  );
	strcpy(xvar[LAT_INDEX].name, "LATITUDE"   );
	strcpy(xvar[P_INDEX].name  , "P"          );
	strcpy(xvar[TEM_INDEX].name, "TEMPERATURE");
	strcpy(xvar[SAL_INDEX].name, "SALINITY"   );

	return(0);
}

/* ------------------------------------------------------------------
	 -Free Data:
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int free_data(EXTRACT_VARIABLE_TYPE *list)
#else
int ()
#endif
{

	int i;

	for(i=0; i<MAX_VARS; i++)
		free_xtract_prf_var(list+i);
	
	return(0);

}

/* ------------------------------------------------------------------
	 -CTD Transducer Data:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int ctd_tr_data(TR_CTD_TYPE *ctd, TR_PARAM_TYPE *tr, FLOAT *data, 
								YMDHMS_TIME_TYPE *start, YMDHMS_TIME_TYPE *end)
#else
int ctd_tr_data()
#endif
{
	
	YMDHMS_TIME_TYPE this_time;
	EXTRACT_VARIABLE_TYPE xvar[MAX_VARS];

	int dbid, dbam, dbmm, ierr, search;
	int n, i, error, n_data, do_fill;
	char str[40];
	
	/* Initializations */
	ierr = 0;
	error = 0;
	init_xvar(xvar);

	/* Initialize CTD database */
	dbid = CTD_DB;
	dbam = READ_ONLY;
	dbmm = DIR_ON_DISK;
	DBOPEN(&dbid, ctd->dbname, &dbam, &dbmm, &ierr);
	DBERROR(&ierr, "DBOPEN CTD DB");
	if(ierr) return(-1);
	
	/* Position CTD database */
	search = TIME_SEARCH;
	DBSRCH(&search, (char *) start, &ierr);
	if(ierr && (ierr != SEARCH_BEFORE_BEGINNING)){
		DBERROR(&ierr, "DBSRCH CTD DB");
		return(-1);
	}
	/* Loop through the given time range */
	n = 0;
	ierr = 0;
	error = 0;
	printf("\n\n");
	while(!ierr && check_time(&this_time, start, end)){

		/* Print current time */
		ymdhms_to_str(str, &this_time, 1);
		printf("CTD: %s\n", str);

		/* Extract the data */
		for(i=0; i<MAX_VARS; i++){
			n_data = extract_prf_var_f(xvar+i, tr->year_base);
			if( error_found(n_data < 0, "Extracting transducer data")){
				error = n_data;
				goto on_error;
			}
		}
	
		/* Special conversions */
		if(tr->bits & LON_80_BIT){
			if(*xvar[LON_INDEX].data > 180.0)
				*xvar[LON_INDEX].data -= 360.0;
		}
		
		/* Fill gaps inside data profile */
		do_fill = do_fill_gaps(xvar[P_INDEX].data, xvar[TEM_INDEX].data, 
													 tr->tr_depth, tr->z_gap, (int) xvar[TEM_INDEX].n);
		if(error_found(do_fill < 0, "While in do_fill_gaps(T)\n")){
			error = -1;
			goto on_error;
		}
		if(do_fill)
			fill_gaps((float *) xvar[TEM_INDEX].data, (int) xvar[TEM_INDEX].n);
		do_fill = do_fill_gaps(xvar[P_INDEX].data, xvar[SAL_INDEX].data, 
													 tr->tr_depth, tr->z_gap, (int) xvar[SAL_INDEX].n);
		if(error_found(do_fill < 0, "While in do_fill_gaps(S)\n")){
			error = -1;
			goto on_error;
		}
		if(do_fill)
			fill_gaps((float *) xvar[SAL_INDEX].data, (int) xvar[SAL_INDEX].n);
		
		/* Save current values and regrid to transducer depth */
		*(data + (T_INDEX*NPF_MAX) + n)   = xvar[T_INDEX].data[0];
		*(data + (LON_INDEX*NPF_MAX) + n) = xvar[LON_INDEX].data[0];
		*(data + (LAT_INDEX*NPF_MAX) + n) = xvar[LAT_INDEX].data[0];
		*(data + (TEM_INDEX*NPF_MAX) + n) = BADFLOAT;
		*(data + (SAL_INDEX*NPF_MAX) + n) = BADFLOAT;

		s_regridli((float *) xvar[P_INDEX].data, (float *) xvar[TEM_INDEX].data, 
							 (float *) &(tr->tr_depth), 
							 (float *) data + (TEM_INDEX*NPF_MAX) + n, 
							 (int) xvar[TEM_INDEX].n, 1);
		s_regridli((float *) xvar[P_INDEX].data, (float *) xvar[SAL_INDEX].data, 
							 (float *) &(tr->tr_depth), 
							 (float *) data + (SAL_INDEX*NPF_MAX) + n,
							 (int) xvar[SAL_INDEX].n, 1);		

		/* Release memory from extraction */
		free_data(xvar);

		/* Move to next profile */
		DBMOVE(&(ctd->steps), &ierr);
		DBERROR(&ierr, "DBMOVE CTD DB");
		n++;

		/* Control profile number */
		if(error_found(n >= NPF_MAX, "Too many profiles\n")){
			error = INSUFFICIENT_MEMORY;
			goto on_error;
		}
	}

on_error:

	/* Release memory */
	free_data(xvar);

	/* Close CTD database */
	DBCLOSE(&ierr);
	DBERROR(&ierr, "DBCLOSE CTD DB");

	/* Return value */
	if(error || ierr)
		return(-1);
	
	return(n);

}

/* ------------------------------------------------------------------
	 -ADCP Transducer Data:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int adcp_tr_data(TR_ADCP_TYPE *adcp, FLOAT *data)
#else
int adcp_tr_data()
#endif
{
	
	FILE *fp;
	int error, n, n_data;
	char buff[150];

	/* Initializations */
	n = 0;
	error = 0;
	fp = check_fopen(adcp->fname, "r");

	/* Loop through the ADCP position file */
	while(getline_nc(fp, buff, 120) >= 0){
		
		/* Get ADCP positions */
		n_data = sscanf(buff, "%f %f %f", 
										data + (T_INDEX*NPF_MAX) + n,
										data + (LON_INDEX*NPF_MAX) + n,
										data + (LAT_INDEX*NPF_MAX) + n);
		if(error_found(n_data != 3, "Reading ADCP time and positions\n")){
			error = READ_ERROR;
			goto on_error;
		}

		n++;
		if(error_found(n >= NPF_MAX, "Too many profiles\n")){
			error = INSUFFICIENT_MEMORY;
			goto on_error;
		}
	}

on_error:

	fclose(fp);
	if(error)
		return(-1);

	return(n);

}

/* ------------------------------------------------------------------
	 MAIN:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
void do_it(FILE *fp_cnt)
#else
void do_it(fp_cnt)
FILE *fp_cnt;
#endif
{

	YMDHMS_TIME_TYPE start, end;
	FILE_NAME_TYPE t_fname, s_fname;
	FILE *fp_t, *fp_s;

	FLOAT tmp_f[MAX_VARS][NPF_MAX];
	int error, ind[NPF_MAX];
	int i, j, k;

	/* Initialize variables and options */
	tr.z_gap = DEFAULT_Z_GAP;
	tr.bits = 0x0;
	ctd.steps = 1;
	error = -1;

	/* Get the parameters and options */
	if(get_parameters(fp_cnt, parameters, name_list))  exit(-1);
	print_parameters(stdout, parameters, name_list, "\n");
	if(execute_options(fp_cnt, options, 0)<=0){
		printf("\nERROR: while reading options.\n\n");
		exit(-1);
	}
	 
	/* Initialize output files */
	new_extension(t_fname, fout, ".tem");
	fp_t = check_fopen(t_fname,"w");
	fprintf(fp_t, "%% TR_CTD:");
	print_parameters(fp_t, parameters, name_list, "\n%");
	fprintf(fp_t,
					"\n%%\n%% Time        ctd_tr_temperature\n");
	
	new_extension(s_fname, fout, ".sal");
	fp_s = check_fopen(s_fname,"w");	
	fprintf(fp_s, "%% LST_CTD:");
	print_parameters(fp_s, parameters, name_list, "\n%");
	fprintf(fp_s,
					"\n%%\n%% Time        ctd_tr_salinity\n");
	
	/* Loop through CTD time ranges */
	while(get_time_range(fp_cnt, &start, &end) == 1){
		
		print_time_range(stdout, &start, &end);
		
		/* Initialize data arrays */
		for(i=0; i<NPF_MAX; i++){
			for(j=0; j<MAX_VARS; j++){
				ctd.data[j][i] = BADFLOAT;
				adcp.data[j][i] = BADFLOAT;
			}
		}
		
		/* Get transducer CTD data for this time range */
		ctd.n = ctd_tr_data(&ctd, &tr, (FLOAT *) ctd.data[0], &start, &end);
		if(error_found(ctd.n <= 0, "Getting CTD data\n"))
			goto on_error;
		
		/* Get ADCP positions */
		adcp.n = adcp_tr_data(&adcp, (FLOAT *) adcp.data[0]);
		if(error_found(adcp.n <= 0, "Getting ADCP data\n"))
			goto on_error;
		
		/* Remove bad CTD axis data */
		 j = 0;
		 for(i=0; i<ctd.n; i++){
			 if(ctd.data[tr.axis][i] < ADJ_BADFLOAT){
				 for(k=0; k<MAX_VARS; k++){
					 tmp_f[k][j] = ctd.data[k][i];
				 }
				 j++;
			 }
		 }
		 if(error_found(j < 2, "No CTD data available\n"))
			 goto on_error;
		 ctd.n = j;
		 for(i=0; i<ctd.n; i++)
			 for(j=0; j<MAX_VARS; j++)
				 ctd.data[j][i] = tmp_f[j][i];

		 /* Remove bad ADCP axis data */
		 j = 0;
		 for(i=0; i<adcp.n; i++){
			 if(adcp.data[tr.axis][i] < ADJ_BADFLOAT){
				 for(k=0; k<MAX_VARS; k++){
					 tmp_f[k][j] = adcp.data[k][i];
				 }
				 j++;
			 }
		 }
		 if(error_found(j == 0, "No valid ADCP data available\n"))
			 goto on_error;
		 adcp.n = j;
		 for(i=0; i<adcp.n; i++)
			 for(j=0; j<MAX_VARS; j++)
				 adcp.data[j][i] = tmp_f[j][i];

		 /* Sort CTD data by axis element */		 
		 sort_float(ctd.n, ind, (FLOAT *) ctd.data[tr.axis]);
		 for(i=0; i<ctd.n; i++){
			 for(j=0; j<MAX_VARS; j++)
				 tmp_f[j][i] = ctd.data[j][ind[i]];
		 }
		 for(i=0; i<ctd.n; i++){
			 for(j=0; j<MAX_VARS; j++)
				 ctd.data[j][i] = tmp_f[j][i];
		 }
		 

		 /* Sort ADCP data by first key element */
		 sort_float(adcp.n, ind, adcp.data[tr.axis]);
		 for(i=0; i<adcp.n; i++){
			 for(j=0; j<MAX_VARS; j++)
				 tmp_f[j][i] = adcp.data[j][ind[i]];
		 }
		 for(i=0; i<adcp.n; i++){
			 for(j=0; j<MAX_VARS; j++)
				 adcp.data[j][i] = tmp_f[j][i];
		 }
		 
		 /* Fill gaps of CTD data */
		 fill_gaps_on_grid((float *) ctd.data[TEM_INDEX], 
											 (float *) ctd.data[tr.axis], ctd.n);
		 fill_gaps_on_grid((float *) ctd.data[SAL_INDEX], 
											 (float *) ctd.data[tr.axis], ctd.n);

		 /* Interpolate CTD data to ADCP positions */
		 s_regridli((float *) ctd.data[tr.axis], (float *) ctd.data[TEM_INDEX],
								(float *) adcp.data[tr.axis], (float *) adcp.data[TEM_INDEX],
								ctd.n, adcp.n);
		 s_regridli((float *) ctd.data[tr.axis], (float *) ctd.data[SAL_INDEX],
								(float *) adcp.data[tr.axis], (float *) adcp.data[SAL_INDEX],
								ctd.n, adcp.n);
		 
		 /* Fill gaps for missing values at beginning and end */
		 fill_gaps((float *) adcp.data[TEM_INDEX], adcp.n);
		 fill_gaps((float *) adcp.data[SAL_INDEX], adcp.n);

		 /* Sort Final ADCP data in ascending time */
		 sort_float(adcp.n, ind, adcp.data[T_INDEX]);
		 for(i=0; i<adcp.n; i++){
			 for(j=0; j<MAX_VARS; j++)
				 tmp_f[j][i] = adcp.data[j][ind[i]];
		 }
		 for(i=0; i<adcp.n; i++){
			 for(j=0; j<MAX_VARS; j++)
				 adcp.data[j][i] = tmp_f[j][i];
		 }		 

		 /* Write data: time and CTD data in first two columns for fix_temp
				            programme. Added in addition columns with positions in
				            case checking of interpolation is required */
		 for(i=0; i<adcp.n; i++){
			 fprintf(fp_t, "%12.6f  %7.4f  %12.6f  %12.6f\n", 
							 adcp.data[T_INDEX][i], adcp.data[TEM_INDEX][i],
							 adcp.data[LON_INDEX][i], adcp.data[LAT_INDEX][i]);
			 fprintf(fp_s, "%12.6f  %7.4f  %12.6f  %12.6f\n",
							 adcp.data[T_INDEX][i], adcp.data[SAL_INDEX][i],
							 adcp.data[LON_INDEX][i], adcp.data[LAT_INDEX][i]);
		 }

	 }
			
	 error = 0;

on_error:
   fclose(fp_t);
	 fclose(fp_s);
	 if(error)
		 printf("\n\nLST_CTD terminated on ERROR\n\n");
	 else
		 printf("\n\nLST_CTD terminated normally\n\n");
		 
   return;
}

