/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Copyright by The HDF Group.                                               *
 * All rights reserved.                                                      *
 *                                                                           *
 * This file is part of H4H5TOOLS. The full H4H5TOOLS copyright notice,      *
 * including terms governing use, modification, and redistribution, is       *
 * contained in the COPYING file, which can be found at the root of the      *
 * source code distribution tree. The full H4H5TOOLS copyright notice can    *
 * also be found at https://www.hdfgroup.org/licenses.  If you do not have   *
 * access to either file, you may request a copy from help@hdfgroup.org.     *
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */


#include "netcdf.h"
#include "hdfeos2.h"
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <algorithm>
#ifndef _WIN32
#include <unistd.h>
#endif


#define suicide()  do { _suicide(__FILE__, __LINE__); _exit(1); } while (0)
void _suicide(const char *fname, int line) { std::cerr << "suicide at " << fname << ":" << line << std::endl; _exit(1); }
#define nomem()  do { _nomem(__FILE__, __LINE__); _exit(1); } while (0)
void _nomem(const char *fname, int line) { std::cerr << "memory exhausted at " << fname << ":" << line << std::endl; _exit(1); }

struct reject {
};

#define lack0_throw()  do { _error0(__FILE__, __LINE__, ctx.eosfullname(), "lack"); throw reject(); } while (0)
#define lack1_throw(msg)  do { _error1(__FILE__, __LINE__, ctx.eosfullname(), "lack", msg); throw reject(); } while (0)
#define lack2_throw(msg, a1)  do { _error2(__FILE__, __LINE__, ctx.eosfullname(), "lack", msg, a1); throw reject(); } while (0)
#define junk1_throw(msg)  do { _error1(__FILE__, __LINE__, ctx.eosfullname(), "junk", msg); throw reject(); } while (0)
#define junk2_throw(msg, a1)  do { _error2(__FILE__, __LINE__, ctx.eosfullname(), "junk", msg, a1); throw reject(); } while (0)
#define incons2_throw(msg, a1)  do { _error2(__FILE__, __LINE__, ctx.eosfullname(), "inconsistency", msg, a1); throw reject(); } while (0)
#define susp2_throw(msg, a1)  do { _error2(__FILE__, __LINE__, ctx.eosfullname(), "suspicious", msg, a1); throw reject(); } while (0)
#define unknown1_throw(msg)  do { _error1(__FILE__, __LINE__, ctx.eosfullname(), "unknown", msg); throw reject(); } while (0)
#define unknown2_throw(msg, a1)  do { _error2(__FILE__, __LINE__, ctx.eosfullname(), "unknown", msg, a1); throw reject(); } while (0)
#define diff1_throw(msg, lv, rv)  do { _error3(__FILE__, __LINE__, ctx.eosfullname(), "difference", msg, lv, rv); throw reject(); } while (0)
#define diff2_throw(msg, a1, lv, rv)  do { _error4(__FILE__, __LINE__, ctx.eosfullname(), "difference", msg, a1, lv, rv); throw reject(); } while (0)
#define diffname1_throw(msg, lv, rv)  do { _error3(__FILE__, __LINE__, ctx.eosfullname(), "different name", msg, lv, rv); throw reject(); } while (0)
#define unimplemented2_throw(msg, a1)  do { _error2(__FILE__, __LINE__, ctx.eosfullname(), "unimplemented", msg, a1); throw reject(); } while (0)

void _error0(const char *fname, int line, const std::string &path, const char *kind) { std::cerr << fname << ':' << line << ' ' << path << ": " << kind << std::endl; }
void _error1(const char *fname, int line, const std::string &path, const char *kind, const char *msg) { std::cerr << fname << ':' << line << ' ' << path << ": " << kind << ": " << msg << std::endl; }
template<typename T> void _error2(const char *fname, int line, const std::string &path, const char *kind, const char *msg, const T &a1) { std::cerr << fname << ':' << line << ' ' << path << ": " << kind << ": " << msg << "; " << a1 << std::endl; }
template<typename T, typename U> void _error3(const char *fname, int line, const std::string &path, const char *kind, const char *msg, const T &a1, const U &a2) { std::cerr << fname << ':' << line << ' ' << path << ": " << kind << ": " << msg << "; " << a1 << "; " << a2 << std::endl; }
template<typename T, typename U, typename V> void _error4(const char *fname, int line, const std::string &path, const char *kind, const char *msg, const T &a1, const U &a2, const V &a3) { std::cerr << fname << ':' << line << ' ' << path << ": " << kind << ": " << msg << "; " << a1 << "; " << a2 << "; " << a3 << std::endl; }

#define PROGRESS(expr)  do { if (ctx.log.verified) { ctx.indent(); std::cout << expr << std::endl; } } while (0)
#define DETECTED(expr)  do { if (ctx.log.verified) { ctx.indent(); std::cout << "detected: " << expr << std::endl; } } while (0)
#define VERIFIED(expr)  do { if (ctx.log.verified) { ctx.indent(); std::cout << "verified: " << expr << std::endl; } } while (0)
#define TOLERATED(expr)  do { if (ctx.log.verified) { ctx.indent(); std::cout << "tolerated: " << expr << std::endl; } } while (0)

static bool lengthcmp(int32 hetype, size_t nctype)
{
	long l = hetype;
	long r = nctype;
	return l != r;
}

static std::string nc4safename(const std::string eosname)
{
	std::string mangled;

	for (std::string::const_iterator i = eosname.begin(); i != eosname.end(); ++i) {
#define IS_DIGIT(c)  ('0' <= (c) && (c) <= '9')
#define IS_ALPHA(c)  (('a' <= (c) && (c) <= 'z') || ('A' <= (c) && (c) <= 'Z'))
#define IS_OTHER(c)  ((c) == '_')
		if (mangled.empty()) {
			if (IS_DIGIT(*i))
				mangled += '_';
		}

		if (IS_DIGIT(*i) || IS_ALPHA(*i) || IS_OTHER(*i))
			mangled += *i;
		else
			mangled += '_';
	}
	return mangled;
}

static int namecmp(const char *eosname, const char *ncname)
{
	std::string mangled = nc4safename(eosname);
	return strcmp(mangled.c_str(), ncname);
}

struct context
{
	std::vector<std::string> eospath;

#define PROGRESS_INTERNAL(expr)  do { if (this->log.progress) { this->indent(); std::cout << expr << std::endl; } } while (0)
	void push_eospath(const std::string &obj)
	{
		PROGRESS_INTERNAL("Verifying " << obj << " {{{");
		this->eospath.push_back(obj);
	}
	void pop_eospath()
	{
		this->eospath.pop_back();
		PROGRESS_INTERNAL("}}}");
	}
	const std::string eosname() const
	{
		return this->eospath.back();
	}
	const std::string eosfullname() const
	{
		std::ostringstream oss;
		for (unsigned int i = 0; i < this->eospath.size(); ++i) {
			oss << this->eospath[i];
			if (i != this->eospath.size() - 1)
				oss << ":";
		}
		return oss.str();
	}

	std::vector<int> netcdfid;
	intalloc netcdfeosdatasets;
	int netcdfgeogroup;
	int netcdfdatagroup;
	int netcdfattrgroup;

	void push_ncid(int id)
	{
		this->netcdfid.push_back(id);
	}
	void pop_ncid()
	{
		this->netcdfid.pop_back();
	}
	int ncid() const
	{
		return this->netcdfid.back();
	}

	void set_ncgeogroupid(int grpid)
	{
		this->netcdfgeogroup = grpid;
	}
	int ncgeogroupid() const
	{
		return this->netcdfgeogroup;
	}
	void set_ncdatagroupid(int grpid)
	{
		this->netcdfdatagroup = grpid;
	}
	int ncdatagroupid() const
	{
		return this->netcdfdatagroup;
	}
	void set_ncattrgroupid(int grpid)
	{
		this->netcdfattrgroup = grpid;
	}
	int ncattrgroupid() const
	{
		return this->netcdfattrgroup;
	}

	context()
	{
		intalloc_init(&this->netcdfeosdatasets);
		this->netcdfgeogroup = -1;
		this->netcdfdatagroup = -1;
		this->netcdfattrgroup = -1;
	}
	~context()
	{
		intalloc_free(&this->netcdfeosdatasets);
	}

	struct log_t
	{
		bool progress;
		bool verified;

		log_t()
		{
			this->progress = true;
			this->verified = true;
		}
	};
	log_t log;

	void indent() const
	{
		for (unsigned int i = 0; i < this->eospath.size(); ++i)
			std::cout << "   ";
	}
};

static int inq_ncgroups(int ncid, intalloc *groups)
{
	int numgrps;
	if (nc_inq_grps(ncid, &numgrps, NULL) != NC_NOERR) suicide();
	if (!intalloc_ready(groups, numgrps)) nomem();
	if (nc_inq_grps(ncid, NULL, groups->s) != NC_NOERR) suicide();
	groups->len = numgrps;
	return numgrps;
}

static void inq_unlimdims(int ncid, std::vector<int> &unlimdims)
{
	int numunlimdim;
	if (nc_inq_unlimdims(ncid, &numunlimdim, NULL) != NC_NOERR) suicide();
	unlimdims.resize(numunlimdim);
	if (nc_inq_unlimdims(ncid, NULL, &unlimdims[0]) != NC_NOERR) suicide();
}

static bool inq_coordinateattr(context &ctx, int ncid, int ncvarid, std::string *value)
{
	nc_type nctype;
	size_t nclen;
	char *buffer = NULL;

	if (nc_inq_att(ncid, ncvarid, "coordinates", &nctype, &nclen) != NC_NOERR) return false;
	if (nctype != NC_CHAR) diff1_throw("type mismatch", NC_CHAR, nctype);
	if ((buffer = (char *)malloc(nclen + 1)) == NULL) nomem();
	buffer[nclen] = 0;
	if (nc_get_att_text(ncid, ncvarid, "coordinates", buffer) != NC_NOERR) suicide();

	if (value) *value = buffer;

	free(buffer);
	return true;
}

static const char *GEOFIELD_VGROUP = "Geolocation Fields";
static const char *DATAFIELD_VGROUP = "Data Fields";
static const char *GRIDATTR_VGROUP = "Grid Attributes";
static const char *SWATHATTR_VGROUP = "Swath Attributes";

static void verify_attribute(context &ctx, const attribute_t *eosattr, int ncvarid)
{
	void *values = NULL;
	nc_type nctype;
	size_t nclen;
	std::ostringstream msgss;

	std::string mangledeosname = nc4safename(eosattr->name.s);
	if (nc_inq_att(ctx.ncid(), ncvarid, mangledeosname.c_str(), &nctype, &nclen) != NC_NOERR) suicide();

	if (nctype == NC_CHAR) {
		if (lengthcmp(eosattr->value.len, nclen)) diff1_throw("length", eosattr->value.len, nclen);
		if ((values = malloc(nclen + 1)) == NULL) nomem();
		char *typednc = (char *)values;
		typednc[nclen] = 0;
		if (nc_get_att_text(ctx.ncid(), ncvarid, mangledeosname.c_str(), typednc) != NC_NOERR) suicide();
		if (strncmp(eosattr->value.s, typednc, eosattr->value.len)) diff1_throw("value", eosattr->value.s, typednc);
		msgss << "length: " << eosattr->value.len << ", value: \"" << eosattr->value.s << "\"";
	}
	else {
		switch (nctype) {
#define GETATTR(tid, suffix, ctype)																											\
			case tid:																																\
			{																																			\
				int bufflen = nclen * sizeof(ctype);																						\
				if (lengthcmp(eosattr->value.len, bufflen)) diff1_throw("length", eosattr->value.len, bufflen);			\
				if ((values = malloc(bufflen)) == NULL) nomem();																		\
				ctype *typedeos = (ctype *)eosattr->value.s;																				\
				ctype *typednc = (ctype *)values;																							\
				if (nc_get_att_##suffix(ctx.ncid(), ncvarid, mangledeosname.c_str(), typednc) != NC_NOERR) suicide();			\
				for (int i = 0; i < (int)nclen; ++i) {																						\
					if (typedeos[i] != typednc[i]) diff2_throw("element", i, typedeos[i], typednc[i]);						\
				}																																		\
				msgss << "count: " << nclen << ", length: " << eosattr->value.len;												\
				break;																																\
			}
			GETATTR(NC_BYTE, schar, signed char)
			GETATTR(NC_SHORT, short, short)
			GETATTR(NC_INT, int, int)
			GETATTR(NC_FLOAT, float, float)
			GETATTR(NC_DOUBLE, double, double)
			GETATTR(NC_UBYTE, ubyte, unsigned char)
			GETATTR(NC_USHORT, ushort, unsigned short)
			GETATTR(NC_UINT, uint, unsigned int)
			GETATTR(NC_INT64, longlong, long long)
			GETATTR(NC_UINT64, ulonglong, unsigned long long)
			default:
				unimplemented2_throw("unknown type", nctype);
				break;
		}
	}

	free(values);

	VERIFIED("attribute " << eosattr->name.s << "(" << mangledeosname << "); " << msgss.str());
}

static void verify_attributes(context &ctx, const attributealloc &attrs)
{
	int ncnumattrs = 0;
	if (nc_inq(ctx.ncid(), NULL, NULL, &ncnumattrs, NULL) != NC_NOERR) suicide();

	for (int i = 0; i < attrs.len; ++i) {
		const attribute_t *eosattr = &attrs.s[i];

		ctx.push_eospath(eosattr->name.s);
		verify_attribute(ctx, eosattr, -1);
		ctx.pop_eospath();
	}
}

static void verify_dimension(context &ctx, const dimension_t *eosdim, int ncdimid, const std::string &ncdimname, bool unlimited)
{
	std::ostringstream msgss;

	size_t nclength;
	if (nc_inq_dimlen(ctx.ncid(), ncdimid, &nclength) != NC_NOERR) suicide();
	if (unlimited) {
		if (eosdim->dimsize != 0) diff1_throw("dimsize (unlimited)", eosdim->dimsize, 0);
	}
	else {
		if (lengthcmp(eosdim->dimsize, nclength)) diff1_throw("dimsize", eosdim->dimsize, nclength);
	}
	msgss << "length: " << eosdim->dimsize << (unlimited ? " (unlimited)" : "");

	VERIFIED("dimension " << eosdim->name.s << "(" << ncdimname << "); " << msgss.str());
}

static void verify_dimensions(context &ctx, const grid_t *grid)
{
	// Unlimited dimension
	std::vector<int> unlimdims;
	inq_unlimdims(ctx.ncid(), unlimdims);

	bool explicitxdim = false, explicitydim = false;
	for (int i = 0; i < grid->dim.len; ++i) {
		const dimension_t *eosdim = &grid->dim.s[i];
		std::string mangledeosname = nc4safename(eosdim->name.s);
		int ncdimid;

		// One unusual file uses nlon and nlat in lieu of standard XDim and YDim.
		// The h4toh5 converts more smartly, but the verifier does not need to be that smart.
		bool isxdim = strcmp(eosdim->name.s, "XDim") == 0 || strcmp(eosdim->name.s, "nlon") == 0;
		bool isydim = strcmp(eosdim->name.s, "YDim") == 0 || strcmp(eosdim->name.s, "nlat") == 0;

		ctx.push_eospath(eosdim->name.s);
		ctx.push_ncid(ctx.ncdatagroupid());
		if (nc_inq_dimid(ctx.ncid(), mangledeosname.c_str(), &ncdimid) != NC_NOERR) {
			// XDim and YDim are implicit, but some files explicitly define them.
			// It's fine that h4toh5 does not convert them.
			if (!isxdim && !isydim) lack0_throw();
		}
		else {
			bool unlimited = std::find(unlimdims.begin(), unlimdims.end(), ncdimid) != unlimdims.end();
			verify_dimension(ctx, eosdim, ncdimid, mangledeosname, unlimited);
		}
		ctx.pop_ncid();
		ctx.pop_eospath();

		if (isxdim) explicitxdim = true;
		else if (isydim) explicitydim = true;
	}

	int ncnumdims = 0;
	// The current h4toh5 is storing all dimensions under Data Fields.
	if (nc_inq(ctx.ncdatagroupid(), &ncnumdims, NULL, NULL, NULL) != NC_NOERR) suicide();

	// XDim and YDim may not be captured by HDF-EOS2 API, but implied.
	int impliciteosdim = 2;
	if (explicitxdim) --impliciteosdim;
	if (explicitydim) --impliciteosdim;

	if (grid->dim.len + impliciteosdim != ncnumdims) diff1_throw("# dimension", grid->dim.len + impliciteosdim, ncnumdims);
}

static void verify_dimensions(context &ctx, const swath_t *swath)
{
	int ncnumdims = 0;

	// Dimensions are stored in the group corresponding to an HDF-EOS2 Data Set.
	// Also, they should be accessible from Data Fields group and Geolocation Fields group.
	if (nc_inq(ctx.ncid(), &ncnumdims, NULL, NULL, NULL) != NC_NOERR) suicide();
	if (swath->dim.len != ncnumdims) diff1_throw("# dimension", swath->dim.len, ncnumdims);

	// Unlimited dimension
	std::vector<int> unlimdims;
	inq_unlimdims(ctx.ncid(), unlimdims);

	for (int i = 0; i < swath->dim.len; ++i) {
		const dimension_t *eosdim = &swath->dim.s[i];
		std::string mangledeosname = nc4safename(eosdim->name.s);
		int ncdimid;

		ctx.push_eospath(eosdim->name.s);
		if (nc_inq_dimid(ctx.ncid(), mangledeosname.c_str(), &ncdimid) != NC_NOERR) lack0_throw();
		bool unlimited = std::find(unlimdims.begin(), unlimdims.end(), ncdimid) != unlimdims.end();
		verify_dimension(ctx, eosdim, ncdimid, mangledeosname, unlimited);
		ctx.pop_eospath();
	}
}

static const char * converted_griddimname(const char *eosdimname, bool orthogonal)
{
	if (orthogonal) {
		if (strcmp(eosdimname, "XDim") == 0) return "lon";
		if (strcmp(eosdimname, "YDim") == 0) return "lat";

		// I could see one unusual file that uses nlon and nlat instead of XDim and YDim.
		if (strcmp(eosdimname, "nlon") == 0) return "lon";
		if (strcmp(eosdimname, "nlat") == 0) return "lat";
	}
	return eosdimname;
}

static void verify_field_coordinates(context &ctx, int ncvarid, const char *coordinates, bool isgrid)
{
	int varnumdims = 0;
	int vardimids[NC_MAX_VAR_DIMS];
	if (nc_inq_var(ctx.ncid(), ncvarid, NULL, NULL, &varnumdims, vardimids, NULL) != NC_NOERR) suicide();

	int geogroup = isgrid ? ctx.ncdatagroupid() : ctx.ncid();
	int length = strlen(coordinates);
	for (int i = 0, j = 0; j <= length; ++j) {
		if ((j == length && length) || coordinates[j] == ' ') {
			std::string dimscale(coordinates + i, j - i);
			std::ostringstream msgss;

			int dimscaleid;
			if (nc_inq_varid(geogroup, dimscale.c_str(), &dimscaleid) != NC_NOERR) incons2_throw("coordinate", dimscale);

			int dimscalenumdims = 0;
			int dimscaledimids[NC_MAX_VAR_DIMS];
			if (nc_inq_var(geogroup, dimscaleid, NULL, NULL, &dimscalenumdims, dimscaledimids, NULL) != NC_NOERR) suicide();
			// If a dimscale refers dim0, dim1, ..., dimn, the field should also refer those dimensions.
			for (int k = 0; k < dimscalenumdims; ++k) {
				int numfound = 0;
				char dimname[NC_MAX_NAME + 1];
				for (int l = 0; l < varnumdims; ++l) {
					if (dimscaledimids[k] == vardimids[l])
						numfound++;
				}
				if (numfound == 0) lack2_throw("dimension referred by coordinate", dimscale);
				if (numfound != 1) susp2_throw("dimension referred twice or more", dimscale);
				if (nc_inq_dimname(geogroup, dimscaledimids[k], dimname) != NC_NOERR) suicide();
				msgss << (k == 0 ? "" : ", ") << dimname;
			}

			VERIFIED("coordinate " << dimscale << "; " << msgss.str());
			i = j + 1;
		}
	}
}

static void verify_field(context &ctx, const fieldinfo_t *eosfield, int ncvarid, const std::string &ncvarname, bool forceorthogonal, bool isgrid, bool isdatafield)
{
	std::ostringstream msgss;
	nc_type nctype;
	int ncnumdims = 0;
	int ncdimids[NC_MAX_VAR_DIMS];
	int ncnumattrs = 0;
	if (nc_inq_var(ctx.ncid(), ncvarid, NULL, &nctype, &ncnumdims, ncdimids, &ncnumattrs) != NC_NOERR) suicide();

	if (eosfield->dims.len != ncnumdims) diff1_throw("# dimension", eosfield->dims.len, ncnumdims);
	int numpoints = 1;
	bool xdimexist = false, ydimexist = false;
	for (int i = 0; i < eosfield->dims.len; ++i) {
		const dimension_t *eosdim = &eosfield->dims.s[i];
		char ncdimname[NC_MAX_NAME + 1];
		size_t ncdimlength;

		if (nc_inq_dim(ctx.ncid(), ncdimids[i], ncdimname, &ncdimlength) != NC_NOERR) suicide();
		if (namecmp(converted_griddimname(eosdim->name.s, forceorthogonal), ncdimname)) diffname1_throw("dimension", eosdim->name.s, ncdimname);
		if (lengthcmp(eosdim->dimsize, ncdimlength)) diff1_throw("dimsize", eosdim->dimsize, ncdimlength);
		numpoints *= eosdim->dimsize;
		msgss << (i == 0 ? "dimension: " : " * ") << eosdim->name.s << '[' << eosdim->dimsize << ']';
		if (strcmp(eosdim->name.s, "XDim") == 0) xdimexist = true;
		if (strcmp(eosdim->name.s, "YDim") == 0) ydimexist = true;
	}

	void *values;
	switch (nctype) {
#define GETVAR(tid, suffix, ctype, usefp)																				\
		case tid:																												\
		{																															\
			int bufflen = numpoints * sizeof(ctype);																	\
			if ((values = malloc(bufflen)) == NULL) nomem();														\
			ctype *typedeos = (ctype *)eosfield->data.s;																\
			ctype *typednc = (ctype *)values;																			\
			if (nc_get_var_##suffix(ctx.ncid(), ncvarid, typednc)	!= NC_NOERR) suicide();					\
			for (int i = 0; i < numpoints; ++i) {																		\
				if (usefp && isnan(typedeos[i]) && isnan(typednc[i])) continue;								\
				if (typedeos[i] != typednc[i]) diff2_throw("element", i, typedeos[i], typednc[i]);		\
			}																														\
			break;																												\
		}
		GETVAR(NC_BYTE, schar, signed char, false)
		GETVAR(NC_CHAR, text, char, false)
		GETVAR(NC_SHORT, short, short, false)
		GETVAR(NC_INT, int, int, false)
		GETVAR(NC_FLOAT, float, float, true)
		GETVAR(NC_DOUBLE, double, double, true)
		GETVAR(NC_UBYTE, ubyte, unsigned char, false)
		GETVAR(NC_USHORT, ushort, unsigned short, false)
		GETVAR(NC_UINT, uint, unsigned int, false)
		GETVAR(NC_INT64, longlong, long long, false)
		GETVAR(NC_UINT64, ulonglong, unsigned long long, false)
		default:
			unimplemented2_throw("unknown type", nctype);
			break;
#undef GETVAR
	}
	free(values);

	bool shouldhavecoordinates = false;
	if (isdatafield) {
		if (isgrid) {
			// XDim and YDim are almost always implied, but one can fabricate bizarre Grid
			if (xdimexist && ydimexist) {
				// Unless forceorthogonal is true, "coordinates" should be given.
				shouldhavecoordinates = !forceorthogonal;
			}
		}
		else {
			// Even if the rank is larger than 1, Swath data field does not always have "coordinates" attribute
			// because those dimensions have nothing to do with geolocations. Since this case is too complicated
			// to handle here, another routine will check this.
		}
	}
	if (shouldhavecoordinates) {
		std::string value;
		if (!inq_coordinateattr(ctx, ctx.ncid(), ncvarid, &value)) lack1_throw("coordinates attribute");
		verify_field_coordinates(ctx, ncvarid, value.c_str(), isgrid);
	}

	VERIFIED("field " << eosfield->name.s << "(" << ncvarname << "); " << msgss.str());
}

static void verify_generated_1d_lonlat(context &ctx, const char *name, int32 eosdimsize, int ncvarid)
{
	std::ostringstream msgss;
	nc_type nctype;
	int ncnumdims = 0;
	int ncdimids[NC_MAX_VAR_DIMS];
	if (nc_inq_var(ctx.ncid(), ncvarid, NULL, &nctype, &ncnumdims, ncdimids, NULL) != NC_NOERR) suicide();

	if (1 != ncnumdims) diff1_throw("# dimension", 1, ncnumdims);
	{
		char ncdimname[NC_MAX_NAME + 1];
		size_t ncdimlength;

		if (nc_inq_dim(ctx.ncid(), ncdimids[0], ncdimname, &ncdimlength) != NC_NOERR) suicide();
		if (strcmp(name, ncdimname)) diffname1_throw("dimension", name, ncdimname);
		if (lengthcmp(eosdimsize, ncdimlength)) diff1_throw("dimsize", eosdimsize, ncdimlength);
		msgss << "length: " << eosdimsize;
	}

	VERIFIED("generated dimscale " << name << "; " << msgss.str());
}

static void verify_generated_2d_lonlat(context &ctx, const char *name, int32 eosxdimsize, int32 eosydimsize, int ncvarid, int &ydimmajor)
{
	std::ostringstream msgss;
	nc_type nctype;
	int ncnumdims = 0;
	int ncdimids[NC_MAX_VAR_DIMS];
	if (nc_inq_var(ctx.ncid(), ncvarid, NULL, &nctype, &ncnumdims, ncdimids, NULL) != NC_NOERR) suicide();

	if (2 != ncnumdims) diff1_throw("# dimension", 2, ncnumdims);
	{
		char ncdimnames[2][NC_MAX_NAME + 1];
		size_t ncdimlengths[2];
		int ym = 0;

		if (nc_inq_dim(ctx.ncid(), ncdimids[0], ncdimnames[0], &ncdimlengths[0]) != NC_NOERR) suicide();
		if (nc_inq_dim(ctx.ncid(), ncdimids[1], ncdimnames[1], &ncdimlengths[1]) != NC_NOERR) suicide();
		if (strcmp(ncdimnames[0], "YDim") == 0) ym = 1;
		if (ydimmajor == -1) ydimmajor = ym;
		else if (ydimmajor != ym) diff1_throw("different major row", ydimmajor, ym);

		int xi = ydimmajor ? 1 : 0;
		int yi = ydimmajor ? 0 : 1;
		if (strcmp(ncdimnames[xi], "XDim")) diffname1_throw("dimension", ncdimnames[xi], "XDim");
		if (strcmp(ncdimnames[yi], "YDim")) diffname1_throw("dimension", ncdimnames[yi], "YDim");
		if (lengthcmp(eosxdimsize, ncdimlengths[xi])) diff1_throw("XDim dimsize", eosxdimsize, ncdimlengths[xi]);
		if (lengthcmp(eosydimsize, ncdimlengths[yi])) diff1_throw("YDim dimsize", eosydimsize, ncdimlengths[yi]);

		msgss << "major: " << (ydimmajor ? 'y' : 'x') << ", length: " << eosxdimsize << ", " << eosydimsize;
	}

	VERIFIED("generated 2d dimscale " << name << "; " << msgss.str());
}

static void verify_fields(context &ctx, const fieldinfoalloc &fields, bool forceorthogonal, bool isgrid, bool isdatafield, int ncgroupid)
{
	for (int i = 0; i < fields.len; ++i) {
		const fieldinfo_t *eosfield = &fields.s[i];
		std::string mangledeosname = nc4safename(eosfield->name.s);
		int ncvarid;

		ctx.push_eospath(eosfield->name.s);
		ctx.push_ncid(ncgroupid);
		if (nc_inq_varid(ctx.ncid(), mangledeosname.c_str(), &ncvarid) != NC_NOERR) lack0_throw();
		verify_field(ctx, eosfield, ncvarid, mangledeosname, forceorthogonal, isgrid, isdatafield);
		ctx.pop_ncid();
		ctx.pop_eospath();
	}
}

static void verify_datafields(context &ctx, const fieldinfoalloc &fields, bool forceorthogonal, bool isgrid)
{
	return verify_fields(ctx, fields, forceorthogonal, isgrid, true, ctx.ncdatagroupid());
}

static void verify_geofields(context &ctx, const fieldinfoalloc &fields, bool isgrid)
{
	return verify_fields(ctx, fields, false, isgrid, false, ctx.ncgeogroupid());
}

static void verify_grid_coordvar(context &ctx, const grid_t *grid, bool orthogonal)
{
	ctx.push_ncid(ctx.ncdatagroupid());
	{
		int ncvarid[2];
		if (nc_inq_varid(ctx.ncid(), "lon", &ncvarid[0]) != NC_NOERR) lack1_throw("1D lon");
		if (nc_inq_varid(ctx.ncid(), "lat", &ncvarid[1]) != NC_NOERR) lack1_throw("1D lat");

		if (orthogonal) {
			verify_generated_1d_lonlat(ctx, "lon", grid->info.xdim, ncvarid[0]);
			verify_generated_1d_lonlat(ctx, "lat", grid->info.ydim, ncvarid[1]);
		}
		else {
			int ydimmajor = -1;
			verify_generated_2d_lonlat(ctx, "lon", grid->info.xdim, grid->info.ydim, ncvarid[0], ydimmajor);
			verify_generated_2d_lonlat(ctx, "lat", grid->info.xdim, grid->info.ydim, ncvarid[1], ydimmajor);
		}
	}
	ctx.pop_ncid();
}

static const dimmap_t * find_dimmap(const dimmapalloc *dimmaps, const char *datadimname, const char *geodimname)
{
	for (int i = 0; i < dimmaps->len; ++i) {
		const dimmap_t *dimmap = &dimmaps->s[i];
		if (strcmp(dimmap->data.s, datadimname)) continue;
		if (strcmp(dimmap->geo.s, geodimname)) continue;
		return dimmap;
	}
	return 0;
}

static void verify_adjusted_lonlat(context &ctx, const fieldinfo_t *eosdatafield, const fieldinfo_t *eosgeofield, const char *mangledname)
{
	int ncsrcvarid, ncsrcnumattrs;
	int ncadjustedvarid, ncnumdims, ncdimids[NC_MAX_VAR_DIMS];

	if (nc_inq_varid(ctx.ncgeogroupid(), mangledname, &ncsrcvarid) != NC_NOERR) suicide();
	if (nc_inq_varid(ctx.ncgeogroupid(), mangledname, &ncadjustedvarid) != NC_NOERR) lack2_throw("adjusted coordinate variable", mangledname);
	if (nc_inq_var(ctx.ncgeogroupid(), ncsrcvarid, NULL, NULL, NULL, NULL, &ncsrcnumattrs) != NC_NOERR) suicide();
	if (nc_inq_var(ctx.ncgeogroupid(), ncadjustedvarid, NULL, NULL, &ncnumdims, ncdimids, NULL) != NC_NOERR) suicide();

	// The data field and the adjusted coordinate variable should have the same shape.
	if (eosdatafield->dims.len < ncnumdims) diff1_throw("# dimension", eosdatafield->dims.len, ncnumdims);
	for (int i = 0; i < ncnumdims; ++i) {
		char ncdimname[NC_MAX_NAME + 1];
		size_t ncdimlength;

		if (nc_inq_dim(ctx.ncgeogroupid(), ncdimids[i], ncdimname, &ncdimlength) != NC_NOERR) suicide();

		bool found = false;
		for (int j = 0; j < eosdatafield->dims.len; ++j) {
			const dimension_t *eosdim = &eosdatafield->dims.s[j];

			if (namecmp(eosdim->name.s, ncdimname) == 0) {
				if (lengthcmp(eosdim->dimsize, ncdimlength)) diff1_throw("dimsize", eosdim->dimsize, ncdimlength);
				found = true;
				break;
			}
		}
		if (!found) unknown2_throw("dimension", ncdimname);
	}

	// The adjusted coordinate variable should have all attributes the source coordinate variable has.
	for (int i = 0; i < ncsrcnumattrs; ++i) {
		char attrname[NC_MAX_NAME + 1];
		nc_type ncsrctype, ncadjustedtype;
		size_t ncsrclen, ncadjustedlen;

		if (nc_inq_attname(ctx.ncgeogroupid(), ncsrcvarid, i, attrname) != NC_NOERR) suicide();
		if (nc_inq_att(ctx.ncgeogroupid(), ncsrcvarid, attrname, &ncsrctype, &ncsrclen) != NC_NOERR) suicide();
		if (nc_inq_att(ctx.ncgeogroupid(), ncadjustedvarid, attrname, &ncadjustedtype, &ncadjustedlen) != NC_NOERR) lack2_throw("adjusted coordinate variable", attrname);
		if (ncsrctype != ncadjustedtype) diff1_throw("adjusted coordinate variable type", ncsrctype, ncadjustedtype);
		if (ncsrclen != ncadjustedlen) diff1_throw("adjusted coordinate variable size", ncsrclen, ncadjustedlen);
	}
}

static void verify_swath_adjustedcoordvar(context &ctx, const swath_t *swath)
{
	// For each datafield, check if this refers to geolocation.
	for (int i = 0; i < swath->datafield.len; ++i) {
		const fieldinfo_t *eosdatafield = &swath->datafield.s[i];
		std::string mangledeosname = nc4safename(eosdatafield->name.s);
		int ncdatavarid;

		ctx.push_eospath(eosdatafield->name.s);
		if (nc_inq_varid(ctx.ncdatagroupid(), mangledeosname.c_str(), &ncdatavarid) != NC_NOERR) lack0_throw();

		std::vector<std::string> eoscoordinates;

		for (int j = 0; j < swath->geofield.len; ++j) {
			const fieldinfo_t *eosgeofield = &swath->geofield.s[j];
			std::ostringstream mangledname;
			mangledname << eosgeofield->name.s;

			bool dimmapused = false;
			bool associated = true;
			for (int k = 0; k < eosgeofield->dims.len; ++k) {
				const dimension_t *eosgeodim = &eosgeofield->dims.s[k];

				bool matched = false;
				for (int l = 0; l < eosdatafield->dims.len; ++l) {
					const dimension_t *eosdatadim = &eosdatafield->dims.s[l];

					if (strcmp(eosgeodim->name.s, eosdatadim->name.s) == 0) {
						mangledname << "_0:1";
//						DETECTED(eosgeofield->name.s << "[" << eosgeodim->name.s << "] is associated");
						matched = true;
						break;
					}

					const dimmap_t *eosdimmap = find_dimmap(&swath->dimmap, eosdatadim->name.s, eosgeodim->name.s);
					if (eosdimmap) {
						dimmapused = true;
						mangledname << "_" << eosdimmap->offset << ":" << eosdimmap->increment;
//						DETECTED(eosgeofield->name.s << "[" << eosgeodim->name.s << "] is associated through dimmap from " << eosdatadim->name.s);
						matched = true;
						break;
					}
				}
				if (!matched) {
					associated = false;
					break;
				}
			}
			if (!associated) continue;

//			DETECTED(eosgeofield->name.s << " is associated as " << (dimmapused ? mangledname.str().c_str() : "it is"));
			std::string mangledncname = nc4safename(mangledname.str());
			if (strcmp(eosgeofield->name.s, "Latitude") && strcmp(eosgeofield->name.s, "Longitude")) continue;
			if (dimmapused)
				verify_adjusted_lonlat(ctx, eosdatafield, eosgeofield, mangledncname.c_str());
			eoscoordinates.push_back(dimmapused ? mangledncname.c_str() : eosgeofield->name.s);
		}

		if (eoscoordinates.empty()) {
			if (inq_coordinateattr(ctx, ctx.ncdatagroupid(), ncdatavarid, 0)) junk1_throw("coordinates attribute");
		}
		else {
			std::string nccoordinates;
			if (!inq_coordinateattr(ctx, ctx.ncdatagroupid(), ncdatavarid, &nccoordinates)) lack1_throw("coordinates attribute");
			int length = (int)nccoordinates.size();
			for (int j = 0, k = 0; k <= length; ++k) {
				if ((k == length && length) || nccoordinates[k] == ' ') {
					std::string ncdimscale = nccoordinates.substr(j, k - j);

					std::vector<std::string>::iterator l = std::find(eoscoordinates.begin(), eoscoordinates.end(), ncdimscale);
					if (l == eoscoordinates.end()) junk2_throw("coordinates", ncdimscale);

					VERIFIED("coordinate " << ncdimscale);
					eoscoordinates.erase(l);

					j = k + 1;
				}
			}
			for (std::vector<std::string>::const_iterator j = eoscoordinates.begin(); j != eoscoordinates.end(); ++j) {
				lack2_throw("coordinates", *j);
			}
		}

		ctx.pop_eospath();
	}
}

static bool get_ncdatasetgroup(context &ctx, int *groupid)
{
	int matchedindex = -1;
	for (int i = 0; i < ctx.netcdfeosdatasets.len; ++i) {
		char groupname[NC_MAX_NAME + 1];

		if (nc_inq_grpname(ctx.netcdfeosdatasets.s[i], groupname) != NC_NOERR) suicide();
		if (namecmp(ctx.eosname().c_str(), groupname) == 0) {
			matchedindex = i;
			break;
		}
	}
	if (matchedindex == -1) return false;
	*groupid = ctx.netcdfeosdatasets.s[matchedindex];
	return true;
}

static bool set_ncfieldgroups(context &ctx, bool isgrid)
{
	// HDF-EOS fields are stored in either Geolocation Fields or Data Fields
	intalloc fieldgroups;
	intalloc_init(&fieldgroups);

	inq_ncgroups(ctx.ncid(), &fieldgroups);
	for (int i = 0; i < fieldgroups.len; ++i) {
		char groupname[NC_MAX_NAME + 1];

		if (nc_inq_grpname(fieldgroups.s[i], groupname) != NC_NOERR) suicide();

		if (namecmp(GEOFIELD_VGROUP, groupname) == 0) ctx.set_ncgeogroupid(fieldgroups.s[i]);
		else if (namecmp(DATAFIELD_VGROUP, groupname) == 0) ctx.set_ncdatagroupid(fieldgroups.s[i]);
		else if (isgrid && namecmp(GRIDATTR_VGROUP, groupname) == 0) ctx.set_ncattrgroupid(fieldgroups.s[i]);
		else if (!isgrid && namecmp(SWATHATTR_VGROUP, groupname) == 0) ctx.set_ncattrgroupid(fieldgroups.s[i]);
		else unknown1_throw(groupname);
	}
	intalloc_free(&fieldgroups);
	return true;
}

static void verify_grid(context &ctx, const grid_t *grid)
{
	int ncgroupid = -1;
	if (!get_ncdatasetgroup(ctx, &ncgroupid)) lack0_throw();

	ctx.push_ncid(ncgroupid);
	set_ncfieldgroups(ctx, true);

	bool orthogonal = true;
	{
		// detect if the projection is two-dimensional or one-dimensional
		int ignoredid;
		do {
			if (nc_inq_dimid(ctx.ncdatagroupid(), "XDim", &ignoredid) != NC_NOERR) break;
			if (nc_inq_dimid(ctx.ncdatagroupid(), "YDim", &ignoredid) != NC_NOERR) break;
			if (nc_inq_varid(ctx.ncdatagroupid(), "lon", &ignoredid) != NC_NOERR) break;
			if (nc_inq_varid(ctx.ncdatagroupid(), "lat", &ignoredid) != NC_NOERR) break;

			orthogonal = false;
		} while (0);

		DETECTED((orthogonal ? "one-dimensional projection" : "two-dimensional projection"));
	}

	verify_attributes(ctx, grid->attr);
	verify_dimensions(ctx, grid);
	verify_datafields(ctx, grid->field, orthogonal, true);
	verify_grid_coordvar(ctx, grid, orthogonal);

	ctx.pop_ncid();

	VERIFIED("grid " << grid->name.s);
}

static void verify_swath(context &ctx, const swath_t *swath)
{
	int ncgroupid = -1;
	if (!get_ncdatasetgroup(ctx, &ncgroupid)) lack0_throw();

	ctx.push_ncid(ncgroupid);
	set_ncfieldgroups(ctx, false);

	verify_attributes(ctx, swath->attr);
	verify_dimensions(ctx, swath);
	verify_datafields(ctx, swath->datafield, false, false);
	verify_geofields(ctx, swath->geofield, false);
	verify_swath_adjustedcoordvar(ctx, swath);

	ctx.pop_ncid();

	VERIFIED("swath " << swath->name.s);
}

static void verify(context &ctx, const char *eosfname, const char *ncfname)
{
	{
		int ncid;
		if (nc_open(ncfname, 0, &ncid) != NC_NOERR) suicide();
		ctx.push_ncid(ncid);
	}

	inq_ncgroups(ctx.ncid(), &ctx.netcdfeosdatasets);
	{
		namelistalloc gridnames;
		namelistalloc_init(&gridnames);
		if (parse_gridnames(eosfname, &gridnames) == -1) suicide();

		int32 readfd;
		if ((readfd = open_grids(eosfname)) == -1) suicide();
		for (int i = 0; i < gridnames.len; ++i) {
			grid_t grid;
			init_grid_t(&grid);
			if (read_grid(readfd, gridnames.s[i].name.s, &grid) == -1) suicide();

			ctx.push_eospath(gridnames.s[i].name.s);
			verify_grid(ctx, &grid);
			ctx.pop_eospath();
			free_grid_t(&grid);
		}
		close_grids(readfd);
		namelistalloc_free(&gridnames);
	}

	{
		namelistalloc swathnames;
		namelistalloc_init(&swathnames);
		if (parse_swathnames(eosfname, &swathnames) == -1) suicide();

		int32 readfd;
		if ((readfd = open_swaths(eosfname)) == -1) suicide();
		for (int i = 0; i < swathnames.len; ++i) {
			swath_t swath;
			init_swath_t(&swath);
			if (read_swath(readfd, swathnames.s[i].name.s, &swath) == -1) suicide();

			ctx.push_eospath(swathnames.s[i].name.s);
			verify_swath(ctx, &swath);
			ctx.pop_eospath();
			free_swath_t(&swath);
		}
		close_swaths(readfd);
		namelistalloc_free(&swathnames);
	}

	ctx.pop_ncid();

	VERIFIED(eosfname);
}

static void usage()
{
	std::cout
		<< "USAGE: [HDF-EOS2 file] [netCDF4 file]\n"
		<< " -p  see verbose progress messages\n"
		<< " -v  see verbose verified objects\n"
		<< std::endl;
}

int main(int argc, char **argv)
{
	bool logprogress = false;
	bool logverified = false;
#ifdef _WIN32
	int optind = 1;
	{
		for ( ; optind < argc; ++optind) {
#define OPTION(s, e)  if (strcmp(argv[optind], s) == 0) { e; continue; }
			OPTION("-p", logprogress = true)
			OPTION("-v", logverified = true)
			break;
#undef OPTION
		}
		if (argc != optind + 2) {
			usage();
			return 1;
		}
	}
#else
	{
		int c;
		while ((c = getopt(argc, argv, "pv")) != -1) {
			switch (c) {
				case 'p': logprogress = true; break;
				case 'v': logverified = true; break;
				default:
					usage();
					return 1;
			}
		}
		if (argc != optind + 2) {
			usage();
			return 1;
		}
	}
#endif

	const char *eos2file = argv[optind];
	const char *nc4file = argv[optind + 1];

	context ctx;
	ctx.log.progress = logprogress;
	ctx.log.verified = logverified;

	PROGRESS(eos2file << " : " << nc4file);
	ctx.push_eospath("/");
	bool succeeded = false;
	try {
		verify(ctx, eos2file, nc4file);
		succeeded = true;
	}
	catch (reject r) {
	}
	ctx.pop_eospath();

	return succeeded ? 0 : 1;
}

// vim:ts=3:sw=3:cindent
