
//
// This source code is based in the 'wgrib' package from the
// US NOAA-NWS / National Weather Service
//
// http://www.cpc.ncep.noaa.gov/products/wesley/wgrib.html
//

#define EXPORT __declspec(dllexport) __stdcall

// ibm2flt       wesley ebisuzaki
//
// v1.1 .. faster
// v1.1 .. if mant == 0 -> quick return
//
double ibm2flt(unsigned char *ibm) {

	int positive, power, mant;
	unsigned int abspower;
	double value, exp;

	mant = (ibm[1] << 16) + (ibm[2] << 8) + ibm[3];
	if(mant == 0)
		return 0.0;

	positive = (ibm[0] & 0x80) == 0;
	power = (int) (ibm[0] & 0x7f) - 64;
	abspower = power > 0 ? power : -power;

	exp = 16.0;
	value = 1.0;
	while(abspower)
	{
		if (abspower & 1)
			value *= exp;
		exp = exp * exp;
		abspower >>= 1;
	}

	if(power < 0)
		value = 1.0 / value;
	value = value * mant / 16777216.0;
	if(positive == 0)
		value = -value;

	return value;

}

//
// w. ebisuzaki
//
//  return x**y
//
//
//  input: double x
//	   int y
//
double int_power(double x, int y) {
	
	double value;
	
	if (y < 0)
	{
		y = -y;
		x = 1.0 / x;
	}
	
	value = 1.0;
	while (y)
	{
		if (y & 1)
			value *= x;
		x *= x;
		y >>= 1;
	}
	
	return value;

}

// Integer conversion 

#define int1(a)				((int) ((a)[0]))
#define int2(a)				((int) (((a)[0] << 8) + ((a)[1])))
#define int3(a)				((int) (((a)[0] << 16) + ((a)[1] << 8) + ((a)[2])))

// GRIB message parameters

#define IS_GRIB(p)			(((p)[0] == 'G') && ((p)[1] == 'R') && ((p)[2] == 'I') && ((p)[3] == 'B'))
#define IS_END(p)			(((p)[0] == '7') && ((p)[1] == '7') && ((p)[2] == '7') && ((p)[3] == '7'))
#define GRIB_VER(p)			(int1(p + 7))
#define SEC_LEN(p)			(int3(p))

// PDS section parameters

#define PDS_LEN(pds)		(int3(pds))
#define PDS_Vsn(pds)		(int1(pds + 3))
#define PDS_Center(pds)		(int1(pds + 4))
	//007	US National Weather Service - NCEP (WMC)
#define PDS_Model(pds)		(int1(pds + 5))
	//081	Analysis from GFS (Global Forecast System)
	//082	Analysis from GDAS (Global Data Assimilation System)
#define PDS_Grid(pds)		(int1(pds + 6))
	//002	 10512-point (144x 73) 2.5 deg global longitude-latitude grid. (1,1) at (0E, 90N), matrix layout. N.B.: prime meridian not duplicated
	//003	 65160-point (360x181) 1.0 deg global longitude-latitude grid. (1,1) at (0E, 90N), matrix layout. N.B.: prime meridian not duplicated
	//004	259920-point (720x361) 0.5 deg global longitude-latitude grid. (1,1) at (0E, 90N), matrix layout, N.B.: prime meridian not duplicated
#define PDS_HAS_GDS(pds)	(((pds)[7] & 128) != 0)
#define PDS_HAS_BMS(pds)	(((pds)[7] & 64) != 0)
#define PDS_PARAM(pds)		(int1(pds + 8))
	//007	HGT		Geopotential height				gpm
	//011	TMP		Temperature						K
	//033	UGRD	u-component of wind				m/s
	//034	VGRD	v-component of wind				m/s
	//039	VVEL	Vertical velocity (pressure)	Pa/s
	//041	ABSV	Absolute vorticity				1/s
	//052	RH		Relative humidity				%
	//153	CLWMR	Cloud water						kg/kg
	//154	O3MR	Ozone mixing ratio				kg/kg
#define PDS_LevelType(pds)	(int1(pds + 9))
	//220	PBLRI	Planetary Boundary Layer (derived from Richardson number)
	//100	ISBL	Isobaric level pressure in hectoPascals (hPa) (2 octets)
#define PDS_Level(pds)		(int2(pds + 10))
#define PDS_Year(pds)		(int1(pds + 12))
#define PDS_Month(pds)		(int1(pds + 13))
#define PDS_Day(pds)		(int1(pds + 14))
#define PDS_Hour(pds)		(int1(pds + 15))
#define PDS_Minute(pds)		(int1(pds + 16))
#define PDS_ForecastTimeUnit(pds) (int1(pds + 17))
	//000	Minute
	//001	Hour
	//002	Day
	//003	Month
	//004	Year
#define PDS_Period(pds)		(int2(pds + 18))
#define PDS_Period1(pds)	(int1(pds + 18))
#define PDS_Period2(pds)	(int1(pds + 19))
#define PDS_TimeRange(pds)	(int1(pds + 20))
	//010	Period1 occupies bytes +18,2; product valid at reference time + Period1
#define PDS_NumAve(pds)		(int2(pds + 21))
#define PDS_NumMissing(pds)	(int1(pds + 23))
#define PDS_Century(pds)	(int1(pds + 24))
#define PDS_Subcenter(pds)	(int1(pds + 25))
#define PDS_DScale(pds)		(int_power(10.0, -int2(pds + 26)))
#define PDS_Year4(pds)		((PDS_Century(pds) - 1) * 100 + PDS_Year(pds))

// GDS section parameters

#define GDS_LEN(gds)		(int3(gds))
#define GDS_NV(gds)			(int1(gds + 3))
#define GDS_PV(gds)			((gds[3] == 0) ? -1 : int2(gds + 4) - 1)
#define GDS_PL(gds)			((gds[4] == 255) ? -1 : (int2(gds + 3) << 2) + int2(gds + 4) - 1)
#define GDS_LatLon(gds)		(gds[5] == 0)
	//000	LatLon
	//001	Mercator
	//002	Gnomonic
	//003	Lambert
	//004	Gaussian
	//005	Polar
	//008	Albers
	//010	RotLL
	//050	Harmonic
	//192	Triangular
#define GDS_NLon(gds)		(int2(gds + 6))
#define GDS_NLat(gds)		(int2(gds + 8))
#define GDS_Lat1(gds)		(int3(gds + 10))
#define GDS_Lon1(gds)		(int3(gds + 13))
#define GDS_Mode(gds)		(int1(gds + 16))
#define GDS_Lat2(gds)		(int3(gds + 17))
#define GDS_Lon2(gds)		(int3(gds + 20))
#define GDS_DLon(gds)		(GDS_Mode(gds) & 128 ? int2(gds + 23) : 0)
#define GDS_DLat(gds)		(GDS_Mode(gds) & 128 ? int2(gds + 25) : 0)
#define GDS_Scan(gds)		(int1(gds + 27))

// BMS section parameters

#define BMS_UnusedBits(bms)	(int1(bms + 3))
#define BMS_StdMap(bms)		(int2(bms + 4))
#define BMS_NxNy(bms)		(BMS_StdMap(bms) ? 0 : ((BMS_Len(bms) << 3) - 48 - BMS_UnusedBits(bms))

// BDS section parameters

#define BDS_BScale(bds)		(int_power(2.0, int2(bds + 4)))
#define BDS_RefValue(bds)	(ibm2flt(bds + 6))
#define BDS_NumBits(bds)	(int1(bds + 10))
#define BDS_MoreFlags(bds)	((bds[3] & 16) != 0)
#define BDS_DataStart(bds)	((int) (11 + (BDS_MoreFlags(bds) ? 3 : 0)))

// Unpack functions

void EXPORT GRIBFill(int npoints, float *data, float value)
{
	for(;npoints;npoints--)
		*data = value;
}

int EXPORT GRIBUnPack(unsigned char *msg, int *npoints, float *data)
{
	
	unsigned char *ptr, *pds, *bms, *gds, *bds, *eds, *nxt, *bits;
	int nlon, nlat, len, n_bits, t_bits, tbits, bbits, jmask, i, n;
	double ref, dscale, bscale, scale;
	
	// Parse message sections
	if(!IS_GRIB(msg) || GRIB_VER(msg) != 1)
		return 1; // No GRIB signature found, GRIB version is not 1
	len = SEC_LEN(msg+4);
	ptr = msg + 8;
	pds = ptr;
	ptr += SEC_LEN(pds);
	if(PDS_HAS_GDS(pds))
	{
		gds = ptr;
		ptr += SEC_LEN(gds);
	}
	else
		gds = 0;
	if(PDS_HAS_BMS(pds))
	{
		bms = ptr;
		ptr += SEC_LEN(bms);
	}
	else
		bms = 0;
	bds = ptr;
	ptr += SEC_LEN(bds);
	eds = ptr;
	ptr += 4;
	nxt = ptr;
	if((gds == 0) || (bms != 0) || (len != (nxt - msg)) || !IS_END(eds))
		return 1; // No grid section found, bitmap is present, inconsistent section size, no end section found
	
	// Calculate grid size and reference values
	if(!GDS_LatLon(gds))
		return 1; // No LatLon grid type
	nlon = GDS_NLon(gds);
	nlat = GDS_NLat(gds);
	n = nlon * nlat;
	dscale = PDS_DScale(pds);
	bscale = BDS_BScale(bds);
	scale = bscale * dscale;
	ref = BDS_RefValue(bds) * dscale;
	n_bits = BDS_NumBits(bds);
	
	// Unpack bitstream and scale values
	if(npoints == 0)
		return 1; // No data count pointer is specified
	else if((*npoints != n) || (data == 0))
		*npoints = n; // If data count does not match or a null data pointer is specified, return data count
	else
	{
		if(*npoints < n)
			n = *npoints;
		bits = bds + BDS_DataStart(bds);
		tbits = bbits = 0;
		if (n_bits > 25)
			return 1;
		jmask = (1 << n_bits) - 1;
		t_bits = 0;
		for (i = 0; i < n; i++)
		{
			if (n_bits - t_bits > 8)
			{
				tbits = (tbits << 16) | (bits[0] << 8) | (bits[1]);
				bits += 2;
				t_bits += 16;
			}
			while (t_bits < n_bits)
			{
				tbits = (tbits << 8) + bits[0];
				bits++;
				t_bits += 8;
			}
			t_bits -= n_bits;
			data[i] = (float) ((double) ((tbits >> t_bits) & jmask) * scale + ref);
		}
	}
	
	return 0;

}
