/**---------------------------------------------------------------------------
 ** 
 ** dnmi_.c -- 
 ** 
 ** Author          : Pierre Jaccard
 ** Created On      : 1999/08/02 09:26:14
 ** Last Modified By: Pierre Jaccard
 ** Last Modified On: 1999/08/02 09:26:29
 ** Update Count    : 1
 ** Directory       : /home/pego/pcd1/codas3c/gfi/src/apps/misc/
 ** Version         : 0.0
 ** Status          : Unknown
 ** ---------------------------------------------------------------------- ** 
 ** DESCRIPTION: 
 ** 
 **    Undocumented.
 ** 
 ** ---------------------------------------------------------------------- ** 
 ** REVISIONS: 
 ** ---------------------------------------------------------------------- ** 
 ** CHANGES: 
 **------------------------------------------------------------------------**/

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

	-Introduction:

	Functions for reading and converting data from DNMI files.

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

#include "dnmi_.h"

/* ------------------------------------------------------------------
	 -Initialize DNMI Structure:
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int init_dnmi(DNMI_PROFILE_TYPE *dnmi)
#else
int ()
#endif
{

	/* Initialize data pointers */
	dnmi->s.data = dnmi->u.data = dnmi->v.data = NULL;
	dnmi->lon = dnmi->lat = NULL;
	dnmi->db.prf.u = dnmi->db.prf.v = dnmi->db.prf.s = NULL;
	dnmi->db.prf.f = NULL;
	dnmi->db.blk.x = dnmi->db.blk.y = NULL;

	return(0);

}

/* ------------------------------------------------------------------
	 -DNMI Free Pointer
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_free_pointer(char **p)
#else
int dnmi_free_pointer()
#endif
{

	if(*p){
		free((void *) (*p));
		*p = NULL;
	}
	
	return(0);

}

/* ------------------------------------------------------------------
	 -DNMI Free Data:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_free_data(DNMI_PROFILE_TYPE *dnmi)
#else
int ()
#endif
{

	dnmi_free_pointer((char **) &(dnmi->s.data));
	dnmi_free_pointer((char **) &(dnmi->u.data));
	dnmi_free_pointer((char **) &(dnmi->v.data));
	dnmi_free_pointer((char **) &(dnmi->lon));
	dnmi_free_pointer((char **) &(dnmi->lat));
	dnmi_free_pointer((char **) &(dnmi->db.prf.u));
	dnmi_free_pointer((char **) &(dnmi->db.prf.v));
	dnmi_free_pointer((char **) &(dnmi->db.prf.s));
	dnmi_free_pointer((char **) &(dnmi->db.prf.f));
	dnmi_free_pointer((char **) &(dnmi->db.blk.x));
	dnmi_free_pointer((char **) &(dnmi->db.blk.y));
		
	return(0);

}

/* ------------------------------------------------------------------
	 -Check DNMI Header:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_chk_hdr(DNMI_HEADER_TYPE *hdr, int code)
#else
int ()
#endif
{

	/* Check if this is the expected variable */
	if(error_found( code != (int) hdr->data_type, "Invalid data code\n"))
		return(-1);
	
	return(0);

}

/* ------------------------------------------------------------------
	 -Compare DNMI Header:

	 This is an additionnal test to check header consistency.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_cmp_hdr(DNMI_HEADER_TYPE *h0, DNMI_HEADER_TYPE *h1)
#else
int ()
#endif
{

	int error;

	error = 
		(h0->producer                  != h1->producer                  ) || 
		(h0->grid_area                 != h1->grid_area                 ) || 
		(h0->prognosis                 != h1->prognosis                 ) || 
		(h0->vertical_coord_type       != h1->vertical_coord_type       ) || 
		(h0->vertical_reference        != h1->vertical_reference        ) || 
		(h0->grid_type                 != h1->grid_type                 ) || 
		(h0->x_grid_number             != h1->x_grid_number             ) || 
		(h0->y_grid_number             != h1->y_grid_number             ) || 
		(h0->origin_year               != h1->origin_year               ) || 
		(h0->origin_month_day          != h1->origin_month_day          ) || 
		(h0->origin_hour               != h1->origin_hour               ) || 
		(h0->x_coord_N_pole            != h1->x_coord_N_pole            ) || 
		(h0->y_coord_N_pole            != h1->y_coord_N_pole            ) || 
		(h0->grid_spacing              != h1->grid_spacing              ) || 
		(h0->y_axis_parallel_longitude != h1->y_axis_parallel_longitude );

	if(error)
		return(-1);

	return(0);

}

/* ------------------------------------------------------------------
	 -Read DNMI Data:

	 This function reads <n> bytes of data from the binary DNMI file and
	 converts it into the current host format. The value of <code> specifies the
	 type in which the data should be converted.

	 It returns the number of converted variables, or -1 on error.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int read_dnmi_data(FILE *fp, char *data, int n, int from_host, int code)
#else
int ()
#endif
{

	unsigned int nv;
	int ierr = 1;
	
	char *buff;

	/* Initialize memory for buff */
	buff = calloc((size_t) n, sizeof(char));
	if(error_found(!buff, "Insufficient memory (buff)\n"))
		return(-1);
	
	/* Read data from the file */
	if((int) fread(buff, sizeof(char), (size_t) n, fp) != n){
		if(feof(fp)) return(-1);
		error_found(1, "Reading binary data from file\n");
		goto on_error;
	}
	
	/* Convert data into correct machine format */
	switch(code)
		{
		case LONG_VALUE_CODE:
		case ULONG_VALUE_CODE:
			nv = (unsigned int) round(((float) n)/((float) sizeof(ULONG)));
			if(error_found((int) convert_array_long(data, buff, from_host,
																							HOST_ENVIRONMENT, nv) != n,
										 "Converting binary data (ULONG)\n"))
				goto on_error;
			break;
		case SHORT_VALUE_CODE:
		case USHORT_VALUE_CODE:
			nv = (unsigned int) round(((float) n)/((float) sizeof(USHORT)));
			if(error_found((int) convert_array_short(data, buff, from_host,
																							 HOST_ENVIRONMENT, nv) != n,
										 "Converting binary data (ULONG)\n"))
				goto on_error;
			break;
		default:
			printf("\nERROR: Unknown type code\n");
			goto on_error;
			break;
		}

	ierr = 0;

on_error:
	if(buff) free(buff);

	if(ierr) return(-1);
	
	return((int) nv);
}

/* ------------------------------------------------------------------
	 -Get DNMI Record:

	 Read a full record from a DNMI 2D barotropic tidal currents file.

	 It seems that header and data array are delimited on both sides with a
	 specific ULONG value. This was not the case in the first file I got from
	 DNMI.

	 Returns the number of grid points.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_dnmi_record(FILE *fp, DNMI_DATA_TYPE *dnmi, int host)
#else
int ()
#endif
{

	ULONG delim;
	int i, n_data, n;

	/* Check if data pointer is null */
	dnmi_free_pointer((char **) &(dnmi->data));
	
	/* Initialize header structure */
	set_byte((UBYTE *) &(dnmi->hdr), (UBYTE) '\0', sizeof(DNMI_HEADER_TYPE));

	/* Get header start delimiter */
	if(read_dnmi_data(fp, (char *) &delim, (int) sizeof(ULONG), host, 
										ULONG_VALUE_CODE) != 1){
		if(feof(fp)) return(-1);
		error_found(1, "Reading delimiter (header start)\n");
		return(-1);
	}
	/* Check if start of header */
	if(error_found(delim != DNMI_HEADER_DELIM,
								 "Not at beginning of header\n"))
		return(-1);

	/* Get header */
	n = read_dnmi_data(fp, (char *) &(dnmi->hdr), (int) sizeof(DNMI_HEADER_TYPE), 
										 host, USHORT_VALUE_CODE);
	if(error_found(n != DNMI_HEADER_SIZE, 
								 "Did not read correct number of data from header\n"))
		return(-1);
	
	/* Check if end of header */
	if(error_found(read_dnmi_data(fp, (char *) &delim, (int) sizeof(ULONG), 
																host, ULONG_VALUE_CODE) != 1,
								 "Reading delimiter (header stop)\n"))
		return(-1);
	if(error_found(delim != DNMI_HEADER_DELIM,
								 "Not at end of header\n"))
		return(-1);

	/* Check if start of data */
	if(error_found(read_dnmi_data(fp, (char *) &delim, (int) sizeof(ULONG), 
																host, ULONG_VALUE_CODE) != 1,
								 "Reading delimiter (data start)\n"))
		return(-1);
	if(error_found(delim != DNMI_ARRAY_DELIM,
								 "Not at beginning of data array\n"))
		return(-1);
	
	/* Get the number of grid data */
	n_data = ((int) dnmi->hdr.x_grid_number) * ((int) dnmi->hdr.y_grid_number);
	if(error_found(!(n_data > 0), "No data present in array\n"))
		return(-1);
	dnmi->data = (SHORT *) calloc((size_t) n_data, sizeof(SHORT));
	if(error_found(!(dnmi->data), "Insufficient memory for DNMI data\n"))
		return(-1);
	for(i=0; i<n_data; i++)
		(dnmi->data)[i] = BADSHORT;
	
	/* Get the data from the file */
	n = read_dnmi_data(fp, (char *) dnmi->data, n_data*sizeof(SHORT),
										 host, SHORT_VALUE_CODE);
	if(error_found(n != n_data, "Did not read correct number of data\n"))
		return(-1);

	/* Check if end of data */
	if(error_found(read_dnmi_data(fp, (char *) &delim, (int) sizeof(ULONG),
																host, ULONG_VALUE_CODE) != 1,
								 "Reading delimiter (data stop)\n"))
		return(-1);
	if(error_found(delim != DNMI_ARRAY_DELIM,
								 "Not at end of data array\n"))
		return(-1);
	
	return(n_data);

}

/* ------------------------------------------------------------------
	 -Get DNMI 2-Dimensional Tide Profile:

	 Returns the number of grid points.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_dnmi_profile(FILE *fp, DNMI_PROFILE_TYPE *dnmi) 
#else
int ()
#endif
{
	
	int n0, n1, host;

	/* Initializations */
	host = dnmi->misc->from_host;

	/* Get surface elevation */
	n0 = get_dnmi_record(fp, &(dnmi->s), host);
	if(n0 < 0){
		if(feof(fp)) return(-1);
		error_found(READ_ERROR, "Reading surface topography\n");
		return(-1);
	}
	/* Check header */
	if(error_found(dnmi_chk_hdr(&(dnmi->s.hdr), DNMI_S_CODE),
								 "Checking header for surface topography\n"))
		return(-1);
	
	/* Get U component */
	n1 = get_dnmi_record(fp, &(dnmi->u), host);
	if(error_found((n1<0) || (n1 != n0), "Reading U component\n"))
		return(-1);
	/* Check header */
	if(error_found(dnmi_chk_hdr(&(dnmi->u.hdr), DNMI_U_CODE),
									"Checking header for U component\n"))
		return(-1);
	/* Compare header */
	if(error_found( dnmi_cmp_hdr(&(dnmi->s.hdr), &(dnmi->u.hdr)),
									"Header for U has changed since S\n"))
		return(-1);

	/* Get V component */
	n1 = get_dnmi_record(fp, &(dnmi->v), host);
	if(error_found((n1<0) || (n1 != n0), "Reading V component\n"))
		return(-1);
	if(error_found( dnmi_chk_hdr(&(dnmi->v.hdr), DNMI_V_CODE),
									"Checking header for V component\n"))
		return(-1);
	/* Compare header */
	if(error_found( dnmi_cmp_hdr(&(dnmi->u.hdr), &(dnmi->v.hdr)),
									"Header for V has changed since S\n"))
		return(-1);
	
	return(n0);

}
	 
/* ------------------------------------------------------------------
	 -DNMI to YMDHMH Time
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_to_ymdhmh_time(YMDHMS_TIME_TYPE *time, DNMI_HEADER_TYPE *hdr)
#else
int ()
#endif
{

	double yd;

	/* Set origin time */
	time->year = hdr->origin_year;
	time->month = (USHORT) floor(((double) hdr->origin_month_day)/100.0);
	time->day = (USHORT) (hdr->origin_month_day - (100*time->month)); 
	time->hour = (USHORT) (((double) hdr->origin_hour)/100.0);
	time->minute = 0;
	time->second = set_hs(0.0);
	
	/* Get current profile time */
	yd = year_day(time, (int) hdr->origin_year);
	yd += hdr->relative_hour/24.0;
	yd_to_ymdhmh_time(yd, (int) hdr->origin_year, time);

	if(error_found(invalid_time(time), "Invalid DNMI time\n"))
		return(-1);

	return(0);

}

/* ------------------------------------------------------------------
	 -DNMI to Lat Lon:
	 
	 Convert index values into DOUBLE latitude and longitude positions. Returns
	 the number of converted grid points, or negative if an error occured.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_to_gpos(DNMI_PROFILE_TYPE *dnmi)
#else
int ()
#endif
{

	int *x, *y;
	int i, n, nx, ny;
	int error = -1;

	/* Release memory */
	dnmi_free_pointer((char **) &(dnmi->lon));
	dnmi_free_pointer((char **) &(dnmi->lat));

	/* Initializations */
	nx = (int) dnmi->s.hdr.x_grid_number;
	ny = (int) dnmi->s.hdr.y_grid_number;
	n  = dnmi->misc->n_data;
	x = y = NULL;

	/* Allocate memory for grid indexes */
	x = (int *) calloc((size_t) nx, sizeof(int));
	if(error_found(!x, "Insufficient memory (x)\n")) goto on_error;
	y = (int *) calloc((size_t) ny, sizeof(int));
	if(error_found(!y, "Insufficient memory (y)\n")) goto on_error;

	/* Allocate memory for geographical positions */
	dnmi->lon = (double *) calloc((size_t) n , sizeof(double));
	if(error_found(!(dnmi->lon), "Insufficient memory (lon)\n")) goto on_error;
	dnmi->lat = (double *) calloc((size_t) n, sizeof(double));
	if(error_found(!(dnmi->lat), "Insufficient memory (lat)\n")) goto on_error;

	/* Select the conversion to use, according to the grid type */
	switch(dnmi->s.hdr.grid_type)
		{
		case POLAR_STEREO_PROJ_CODE:
			/* Rectangular grid */
			for(i=0; i<nx; i++)
				x[i] = i+1;
			for(i=0; i<ny; i++)
				y[i] = i+1;
			/* Convert grid indexes to latitude and longitudes */
			error = 
				inv_polar_stereo_proj((double) (dnmi->s.hdr.x_coord_N_pole), 
															(double) (dnmi->s.hdr.y_coord_N_pole), 
															(double) (dnmi->misc->grid_points), 
															(double) (dnmi->s.hdr.y_axis_parallel_longitude),
															x, y, nx, ny, dnmi->lon, dnmi->lat);
			break;
		default:
			error_found(1, "Unknown projection\n");
			break;
		}

on_error:

	/* Release memory */
	dnmi_free_pointer((char **) &x);
	dnmi_free_pointer((char **) &y);

	if(error) return(error);

	return(n);

}

/* ------------------------------------------------------------------
	 -DNMI Set xy Index:

	 Returns the number of good indexes.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int dnmi_set_xy(DNMI_PROFILE_TYPE *dnmi)
#else
int dnmi_set_xy()
#endif
{

	DNMI_POS_RNG_TYPE *p;
	DNMI_RNG_TYPE     *r;
	DNMI_MISC_TYPE    *m;
	
	DMSH_LL_TYPE pos;
	LONG hun;
	int i, lon_180, n;

	/* Initializations */
	m = dnmi->misc;
	r = dnmi->rng;
	p = &(r->pos);
	lon_180 = (m->bits & LONGITUDE_180_BIT) ? -1 : 0;

	/* Release memory */
	dnmi_free_pointer((char **) &(r->index));

	/* Allocate index memory */
	r->index = calloc((size_t) m->n_data, sizeof(int));
	if(error_found(!(r->index), "Insufficient memory (index)\n"))
		return(-1);
	for(i=0; i<(m->n_data); i++)
		r->index[i] = -1;
	
	/* Loop through the positions */
	n = 0;
	for(i=0; i<(m->n_data); i++){
		hun = (LONG)(360000.0 * (dnmi->lon)[i]);
		HUNPOS(&(pos.lon), &hun);
		hun = (LONG)(360000.0 * (dnmi->lat)[i]);
		HUNPOS(&(pos.lat), &hun);
		if(position_inside_range(&(p->min.lon), &(p->max.lon), 
														 &(pos.lon), lon_180) &&
			 position_inside_range(&(p->min.lat), &(p->max.lat), 
														 &(pos.lat), 0)){
			(r->index)[n] = i;
			n++;
		}
	}
	
	return(n);
	
}
