#include "dccrad.h"

/************************************************************************
 * CR_DCOD								*
 *									*
 * This function decodes Canadian radar data into BUFR format.		*
 *									*
 * CR_DCOD ( CRADFL, LUNOUT )						*
 *									*
 * Input parameters:							*
 *	CRADF		char*		Canadian radar data file	*
 *	LUNOUT		int		FORTRAN logical unit number for	*
 *					BUFR output file		*
 **									*
 * Log:									*
 * J. Ator/NCEP		03/11						*
 * J. Ator/NCEP		12/11	Added HSALG to BUFR output		*
 * J. Ator/NCEP		01/12	Only decode HREF from CONVOL files	*
 * J. Ator/NCEP		02/14	Check for NULL return from		*
 *				RSL_anyformat_to_radar			*
 ***********************************************************************/
void cr_dcod( char *cradfl, int lunout ) {

    Radar *radar;
    Volume *volume;
    Sweep *sweep;
    Ray *ray, *swray;

    struct alg_tableEntry {  
	char sid[4];	/* Station ID */
	f77r8 alg;	/* Height of antenna above ground level
			   at station */
    }
    /*  Note that the following entries MUST be sorted in ascending
	order with respect to the sid field. */
    alg_table[] = {
	{ "WBI", 27.5 },
	{ "WGJ", 23.1 },
	{ "WHK", 14.1 },
	{ "WHN", 17.1 },
	{ "WKR", 30.5 },
	{ "WMB", 20 },
	{ "WMN", 30 },
	{ "WSO", 23.5 },
	{ "WTP", 13.9 },
	{ "WUJ", 20.4 },
	{ "WVY", 20.3 },
	{ "WWW", 29 },
	{ "XAM", 15.5 },
	{ "XBE", 15.1 },
	{ "XBU", 14 },
	{ "XDR", 29 },
	{ "XFT", 21.5 },
	{ "XFW", 14.1 },
	{ "XGO", 21.5 },
	{ "XLA", 26.5 },
	{ "XMB", 20.4 },
	{ "XME", 15.5 },
	{ "XNC", 23.1 },
	{ "XNI", 23.1 },
	{ "XPG", 14.1 },
	{ "XRA", 14.1 },
	{ "XSI", 31.5 },
	{ "XSM", 16 },
	{ "XSS", 26.1 },
	{ "XTI", 22.1 },
	{ "XWL", 15.5 }
    };

    char wkary[NUMELEM(alg_table)][4];
    char (*psid)[4];  /* psid is a pointer to an array of 4 chars */

    float x, y;

    int i, j, k, m, ier;
    int dattim;
    int narray, nlevel;
    int iforce = 0;
    int sbtypb;

    unsigned long nsubset = 0;

    char *pc;
    char logmsg[100];
    char subtyp[9];
    char mnemstr1[9];
    char mnemstr2[9];
    char mnemstr3[9];
    char idstr[7] = "canrad";
    char craddn[DCMXLN];
    char cradbn[DCMXLN];
    
    double rlat, rlon;

    f77r8 r8array1[MAX_DATABINS];
    f77r8 r8array2[MAX_DATABINS];
    f77r8 r8array3[MAX_DATABINS];
    f77r8 hsalg;

    /*  Declare a union to use in converting the radar identifer from char to f77r8 */
    union {
	char crid[8];
	f77r8 frid;
    } cnvid;

    /*
    **  Start a log entry for this file.
    */
    cfl_path( cradfl, craddn, cradbn, &ier );
    strcpy( logmsg, "Canadian radar filename: ");
    dc_wclg( 2, "DC", 2, strcat( logmsg, cradbn ), &ier );

    /*
    **  Open the file.
    */
    radar = RSL_anyformat_to_radar(cradfl, NULL);
    if ( radar == NULL ) {
	sprintf( logmsg, "%s", "ERROR: bad radar file - no BUFR subsets generated!" );
	dc_wclg( 2, "DC", 2, logmsg, &ier );
	return;
    }

    /*
    **  Calculate the site latitude and longitude.
    */
    rlat = ( ( ( (double) radar->h.lats / 60. ) + (double) radar->h.latm ) / 60. ) +
		(double) radar->h.latd;
    rlon = ( ( ( (double) radar->h.lons / 60. ) + (double) radar->h.lonm ) / 60. ) +
		(double) radar->h.lond;

    /*
    **  Get the radar identifier from the file name.
    */
    if ( ( pc = strstr( cradfl, "URP:" ) ) != NULL ) {
	cnvid.crid[0] = pc[4];
	cnvid.crid[1] = pc[5];
	cnvid.crid[2] = pc[6];
	cnvid.crid[3] = '\0';
    }

    /*
    **  Get the height of the antenna above ground level.
    */
    for ( i = 0; i < NUMELEM(alg_table); i++ ) {
	strcpy( wkary[i], alg_table[i].sid );
	/* the sid values need to be in a contiguous array for use by bsearch */
    }
    psid = ( char (*)[4] ) bsearch( cnvid.crid, wkary, NUMELEM(alg_table), 4,
			( int (*) ( const void *, const void * ) ) strcmp );
    hsalg = ( ( psid == NULL ) ? BMISS : alg_table[psid-wkary].alg );

    /*
    **  Loop through the radar structure.  It consists of one or more volumes, each
    **  of which consists of one or more sweeps, each of which consists of one or
    **  more rays, each of which consists of one or more data bins.
    */

    for ( i = 0; i < radar->h.nvolumes; i++ ) {
	if ( ( volume = radar->v[i] ) == NULL ) continue;  /* skip NULL volumes */
	if ( ( i == 4 ) && ( strstr( cradfl, "CONVOL" ) != NULL ) ) {  /* only decode reflectivity from CONVOL files */
	    sbtypb = 80;
	    strcpy ( mnemstr2, "HREF");
	}
	else if ( ( i == 1 ) && ( strstr( cradfl, "DOPVOL" ) != NULL ) ) {  /* only decode velocity from DOPVOL files */
	    sbtypb = 110;
	    strcpy ( mnemstr2, "DMVR");
	    strcpy ( mnemstr3, "DVSW");
	}
	else {
	    continue;  /* skip everything else */
	}

	for ( j = 0; j < volume->h.nsweeps; j++ ) {
	    if ( ( sweep = volume->sweep[j] ) == NULL ) continue;  /* skip NULL sweeps */

	    for ( k = 0; k < sweep->h.nrays; k++ ) {
		if ( ( ray = sweep->ray[k] ) == NULL ) continue;  /* skip NULL rays */

		/*
		**  Assume a non-NULL ray at this point, and begin a new BUFR
		**  subset to hold the data values.
		*/
		dattim = ( ray->h.year * 1000000 ) + ( ray->h.month * 10000 ) +
			( ray->h.day * 100 ) + ray->h.hour;
	        sprintf ( subtyp, "%s%.3d", "NC006", ( sbtypb + ray->h.hour ) );
		openmb ( &lunout, subtyp, &dattim, strlen( subtyp ) );

		narray = 3;
		nlevel = 1;

		/*
		**  Date.
		*/
		strcpy ( mnemstr1, "YYMMDD" );
		r8array1[0] = ( f77r8 ) ray->h.year;
		r8array1[1] = ( f77r8 ) ray->h.month;
		r8array1[2] = ( f77r8 ) ray->h.day;
		ufbseq ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Time.
		*/
		strcpy ( mnemstr1, "HHMMSS" );
		r8array1[0] = ( f77r8 ) ray->h.hour;
		r8array1[1] = ( f77r8 ) ray->h.minute;
		r8array1[2] = ( f77r8 ) ray->h.sec;
		ufbseq ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		narray = 2;

		/*
		**  Latitude and longitude.
		*/
		strcpy ( mnemstr1, "LTLONC" );
		r8array1[0] = ( f77r8 ) rlat;
		r8array1[1] = ( f77r8 ) rlon;
		ufbseq ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		narray = 1;

		/*
		**  Radar identifier.
		*/
		if ( pc != NULL ) {
		    strcpy ( mnemstr1, "SSTN" );
		    r8array1[0] = cnvid.frid;
		    ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 
		}

		/*
		**  Station elevation.
		*/
		strcpy ( mnemstr1, "HSMSL" );
		r8array1[0] = ( f77r8 ) radar->h.height;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Height of antenna above station ground.
		*/
		strcpy ( mnemstr1, "HSALG" );
		r8array1[0] = ( f77r8 ) hsalg;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Restrictions on redistribution.
		*/
		strcpy ( mnemstr1, "RSRD" );
		r8array1[0] = ( f77r8 ) 128;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Antenna elevation.
		*/
		strcpy ( mnemstr1, "ANEL" );
		r8array1[0] = ( f77r8 ) ( ray->h.elev < 180. ? ray->h.elev : ray->h.elev - 360. );
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Antenna azimuth.
		*/
		strcpy ( mnemstr1, "ANAZ" );
		r8array1[0] = ( f77r8 ) ray->h.azimuth;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Pulse repetition frequency.
		*/
		strcpy ( mnemstr1, "PRFR" );
		r8array1[0] = ( f77r8 ) ray->h.prf;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Pulse width.  These values are stored in micro-seconds in the RSL output, but
		**  need to be stored in seconds in the BUFR output.
		*/
		strcpy ( mnemstr1, "PWID" );
		r8array1[0] = ( f77r8 ) ray->h.pulse_width / 1000000.;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Beam width.
		*/
		strcpy ( mnemstr1, "BEAMW" );
		r8array1[0] = ( f77r8 ) ray->h.beam_width;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Frequency bandwidth.  These values are stored in MHz in the RSL output, but
		**  need to be stored in Hz in the BUFR output.
		*/
		strcpy ( mnemstr1, "IFRB" );
		r8array1[0] = ( f77r8 ) ray->h.frequency * 1000000.;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Nyquist velocity.
		*/
		strcpy ( mnemstr1, "HNQV" );
		r8array1[0] = ( f77r8 ) ray->h.nyq_vel;
		ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 

		/*
		**  Data values along the ray.
		*/
		nlevel = 0;
		strcpy ( mnemstr1, "DIST" );

		if ( i == 1 ) {
		    /*
		    **  Look for a spectral width ray in volume 2 which has the same metadata
		    **  as this radial velocity ray.
		    */
		    swray = cr_gswr( radar, ray );
		}

		for ( m = 0; m < ray->h.nbins; m++ ) {
		    x = ray->h.f(ray->range[m]);
		    if ( ( x < MISSING_VALUE ) && ( nlevel < MAX_DATABINS ) ) {
			r8array1[nlevel] = ( f77r8 ) ray->h.range_bin1 + ( (m+1) * ray->h.gate_size );
			r8array2[nlevel] = ( f77r8 ) x;
			if ( ( i == 1 ) && ( swray != NULL ) ) {
			  /*
			  ** Check whether this bin has a spectral width.  Note that spectral widths
			  ** are stored as m**2/s**2 in the RSL library output, but they need to be
			  ** stored as m/s in the BUFR output.
			  */
			  if ( ( m < swray->h.nbins ) &&
			      ( y = swray->h.f(swray->range[m]) ) < MISSING_VALUE ) {
			    r8array3[nlevel] = sqrt( ( f77r8 ) y );
			  }
			  else {
			    r8array3[nlevel] = BMISS;
			  }
			}
			nlevel++;
		    }
		}

		if ( nlevel > 0 ) {
		    nsubset++;
		    ufbint ( &lunout, r8array1, &narray, &nlevel, &ier, mnemstr1, strlen ( mnemstr1 ) ); 
		    ufbint ( &lunout, r8array2, &narray, &nlevel, &ier, mnemstr2, strlen ( mnemstr2 ) ); 
		    if ( i == 1 ) {
		      ufbint ( &lunout, r8array3, &narray, &nlevel, &ier, mnemstr3, strlen ( mnemstr3 ) ); 
		    }
		    /*
		    **  Write out the BUFR subset.
		    */
		    ut_wbfr ( &lunout, idstr, &iforce, &ier, strlen ( idstr ) );
		}
	    }
	}
    }

    /*
    **  Make sure all BUFR output has been written out.
    */
    iforce = 1;
    ut_wbfr ( &lunout, idstr, &iforce, &ier, strlen ( idstr ) );

    /*
    **  Write a subset count to the log.
    */
    sprintf( logmsg, "%s %lu %s", "generated", nsubset, "BUFR subsets" );
    dc_wclg( 2, "DC", 2, logmsg, &ier );
}