#include "nq_basic.h"
#include "nq_structure.h"

#include <netcdf.h>
#include <math.h>

#include <visu_tools.h>
#include <visu_elements.h>
#include <visu_data.h>
#include <extraFunctions/scalarFields.h>
#include <coreTools/toolFileFormat.h>
#include <coreTools/toolMatrix.h>

/* Local methods */
static gboolean loadNQETSF(VisuData *data, const gchar* filename,
			   FileFormat *format, GError **error);

RenderingFormatLoad* nqStructuralInit()
{
  char *type[] = {"*.nc", "*-etsf.nc", (char*)0};
  char *descr = _("ETSF file format");
  RenderingFormatLoad *meth;
  
  meth = g_malloc(sizeof(RenderingFormatLoad));
  meth->name = "ETSF (Nanoquanta) file format";
  meth->fmt = fileFormatNew(descr, type);
  if (!meth->fmt)
    g_error("Can't initialize the Nanoquanta loading method, aborting...\n");
  meth->priority = 5;
  meth->load = loadNQETSF;

  return meth;
}

static gboolean loadNQETSF(VisuData *data, const gchar* filename,
			   FileFormat *format, GError **error)
{
  gboolean res;
  int netcdfId, status, i, j, k, ret;
  int varId;
  size_t dimSize;
  char* varNbElement = "number_of_atom_species";
  int nbEle;
  char* varNbNode = "number_of_atoms";
  int nbNode;
  char* varAtomsNames = "atom_species_names";
  int varIdSymbols;
  size_t dimsSymbols[2];
  char* symbols;
  char* varNodeToEle = "atom_species";
  int varIdNodeToEle;
  size_t dimsNodeToEle;
  int *nodeToEle;
  char* varCoord = "reduced_atom_positions";
  int varIdCoord;
  size_t dimsCoord[2];
  double* coord;
  char *varRprimd = "primitive_vectors";
  int varIdRprimd;
  size_t dimsRprimd[2];
  double rprimd[3][3];
  size_t start[] = {0, 0};
  char *varTitle = "title";
  size_t sizeTitle;
  char title[256];
  int *nbNodesPerEle;
  VisuElement **visuEle;
  float boxGeometry[6];
  float vect[3], mat[3][3], vectB[3];
  nc_type ncType;
  gchar *infoUTF8, *name;
  gboolean *flag;

  g_return_val_if_fail(error && *error == (GError*)0, FALSE);
  g_return_val_if_fail(data && filename, FALSE);

  res = nqOpen_netcdfFile(filename, &netcdfId, error);
  if (!res)
    return FALSE;

  /* From now on, the file is a NETCDF valid file, so
     we return TRUE, but errrors can still occurs. */

  /********************************/
  /* Retrieve the numeral values. */
  /********************************/
  /* Grep the number of elements. */
  if (!nqGetDim(netcdfId, error, varNbElement, &varId, &dimSize))
    {
      nqClose_netcdfFile(netcdfId);
      return TRUE;
    }
	nbEle = (int)dimSize;
  /* Grep the number of nodes. */
  if (!nqGetDim(netcdfId, error, varNbNode, &varId, &dimSize))
    {
      nqClose_netcdfFile(netcdfId);
      return TRUE;
    }
  nbNode = (int)dimSize;

  /**************************************************************/
  /* Check the conformance of the arrays to the specifications. */
  /**************************************************************/
  /* Check the names of elements. */
  dimsSymbols[0] = nbEle;
  dimsSymbols[1] = 80;
  if (!nqCheckVar(netcdfId, error, varAtomsNames, &varIdSymbols,
		 NC_CHAR, 2, dimsSymbols))
		{
      nqClose_netcdfFile(netcdfId);
      return TRUE;
    }
  /* Check the 'node to element'array. */
  dimsNodeToEle = nbNode;
  if (!nqCheckVar(netcdfId, error, varNodeToEle, &varIdNodeToEle,
		 NC_INT, 1, &dimsNodeToEle))
		{
      nqClose_netcdfFile(netcdfId);
      return TRUE;
    }
  /* Check the reduce coordinates array. */
  dimsCoord[0] = nbNode;
  dimsCoord[1] = 3;
  if (!nqCheckVar(netcdfId, error, varCoord, &varIdCoord,
		 NC_DOUBLE, 2, dimsCoord))
		{
      nqClose_netcdfFile(netcdfId);
      return TRUE;
    }
  /* Check the rprimd matrix. */
  dimsRprimd[0] = 3;
  dimsRprimd[1] = 3;
  if (!nqCheckVar(netcdfId, error, varRprimd, &varIdRprimd,
		 NC_DOUBLE, 2, dimsRprimd))
		{
      nqClose_netcdfFile(netcdfId);
      return TRUE;
    }

  /****************************/
  /* Retrieve now the arrays. */
  /****************************/
  /* Grep the names of elements. */
  symbols = g_malloc(sizeof(char) * 80 * nbEle);
  status = nc_get_vara_text(netcdfId, varIdSymbols, start, dimsSymbols, symbols);
  if (status != NC_NOERR)
    {
      *error = g_error_new(NQ_ERROR, NQ_ERROR_FILE_FORMAT,
			   _("Retrieve value for variable '%s': %s."),
			   varAtomsNames, nc_strerror(status));
      nqClose_netcdfFile(netcdfId);
      g_free(symbols);
      return TRUE;
    }
  for (i = 0; i < nbEle; i++)
    symbols[(i + 1) * 80 - 1] = '\0';
  /* Grep table from 'node to element'. */
  nodeToEle = g_malloc(sizeof(int) * nbNode);
  status = nc_get_vara_int(netcdfId, varIdNodeToEle, start, &dimsNodeToEle, nodeToEle);
  if (status != NC_NOERR)
    {
      *error = g_error_new(NQ_ERROR, NQ_ERROR_FILE_FORMAT,
			   _("Retrieve value for variable '%s': %s."),
			   varNodeToEle, nc_strerror(status));
      nqClose_netcdfFile(netcdfId);
      g_free(symbols);
      g_free(nodeToEle);
      return TRUE;
    }
  for (i = 0; i < nbNode; i++)
    {
      /* Specifications 1.3 says that the values range from 1.
         So we remove 1 to be coherent with C numbering. */
      nodeToEle[i] -= 1;
      if (nodeToEle[i] < 0 || nodeToEle[i] >= nbEle)
	{
	  *error = g_error_new(NQ_ERROR, NQ_ERROR_FILE_FORMAT,
			       _("Error in indexing array '%s', index out of bounds."),
			       varNodeToEle);
	  nqClose_netcdfFile(netcdfId);
	  g_free(symbols);
	  g_free(nodeToEle);
	  return TRUE;
	}
    }
  /* Grep the reduced coordinates. */
  coord = g_malloc(sizeof(double) * 3 * nbNode);
  status = nc_get_vara_double(netcdfId, varIdCoord, start, dimsCoord, coord);
  if (status != NC_NOERR)
    {
      *error = g_error_new(NQ_ERROR, NQ_ERROR_FILE_FORMAT,
			   _("Retrieve value for variable '%s': %s."),
			   varCoord, nc_strerror(status));
      nqClose_netcdfFile(netcdfId);
      g_free(symbols);
      g_free(nodeToEle);
      g_free(coord);
      return TRUE;
    }
  /* Grep the box definition. */
  status = nc_get_vara_double(netcdfId, varIdRprimd, start, dimsRprimd, &rprimd[0][0]);
  if (status != NC_NOERR)
    {
      *error = g_error_new(NQ_ERROR, NQ_ERROR_FILE_FORMAT,
			   _("Retrieve value for variable '%s': %s."),
			   varRprimd, nc_strerror(status));
      nqClose_netcdfFile(netcdfId);
      g_free(symbols);
      g_free(nodeToEle);
      g_free(coord);
      return TRUE;
    }

  /* Prepare the arrays for memory allocation in the VisuData object. */
  res = matrix_reducePrimitiveVectors(boxGeometry, rprimd);
  if (!res)
    {
      *error = g_error_new(NQ_ERROR, NQ_ERROR_FILE_FORMAT,
			   _("The variable '%s' is not well formed (the basis is not 3D)."),
			   varRprimd);
      nqClose_netcdfFile(netcdfId);
      g_free(symbols);
      g_free(nodeToEle);
      g_free(coord);
      return TRUE;
    }
  nbNodesPerEle = g_malloc(sizeof(int) * nbEle);
  for (i = 0; i < nbEle; i++)
    nbNodesPerEle[i] = 0;
  for (i = 0; i < nbNode; i++)
    nbNodesPerEle[nodeToEle[i]] += 1;
  visuEle = g_malloc(sizeof(VisuElement*) * nbEle);
  for (i = 0; i < nbEle; i++)
    {
      name = g_strstrip(symbols + i * 80);
      visuEle[i] = visuElementGet_fromName(name);
      if (!visuEle[i])
        {
          visuEle[i] = visuElementNewWithName(name);
          visuElementAdd(visuEle[i]);
        }
    }

  DBG_fprintf(stderr, "NQ structure : File '%s' parsed.\n", filename);
  DBG_fprintf(stderr, " | box definition : ( %f %f %f )\n",
	      rprimd[0][0], rprimd[0][1], rprimd[0][2]);
  DBG_fprintf(stderr, " |                  ( %f %f %f )\n",
	      rprimd[1][0], rprimd[1][1], rprimd[1][2]);
  DBG_fprintf(stderr, " |                  ( %f %f %f )\n",
	      rprimd[2][0], rprimd[2][1], rprimd[2][2]);
  DBG_fprintf(stderr, " | number of nodes per element (%d):\n", nbEle);
  for (i = 0; i < nbEle ; i++)
    DBG_fprintf(stderr, " |  %d nodes for '%s' (%p)\n", nbNodesPerEle[i],
		visuEle[i]->name, (gpointer)visuEle[i]);

  g_free(symbols);

  visuDataSet_boxGeometry(data, boxGeometry, TRUE);
  k = 0;
  for (j = 0; j < 3; j++)
    for (i = 0; i < 3; i++)
      {
	if (i > j)
	  mat[i][j] = 0.;
	else
	  mat[i][j] = boxGeometry[k++];
      }
  ret = visuDataSet_population(data, nbEle, nbNodesPerEle, visuEle);
  if (!ret)
    g_error("Can't store the nodes in the VisuData object.\n");
  for (i = 0; i < nbNode; i++)
    {
      for (j = 0; j < 3; j++)
	vectB[j] = (float)coord[3 * i + j];
      matrix_productVector(vect, mat, vectB);
      visuDataAdd_VisuElement(data, visuEle[nodeToEle[i]], vect[0], vect[1], vect[2], i);
    }

  g_free(nodeToEle);
  g_free(coord);
  g_free(nbNodesPerEle);
  g_free(visuEle);

  /**********************/
  /* Optional elements. */
  /**********************/
  status = nc_inq_att(netcdfId, NC_GLOBAL, varTitle, &ncType, &sizeTitle);
  if (status == NC_NOERR && ncType == NC_CHAR && sizeTitle < 255)
    {
  	  status = nc_get_att_text(netcdfId, NC_GLOBAL, varTitle, title);
  	  if (status == NC_NOERR)
  	    {
  	      title[sizeTitle] = '\0';
  	      infoUTF8 = getStringInUTF8(title);
  	      visuDataSet_fileCommentary(data, infoUTF8);
  	      g_free(infoUTF8);
  	    }
    }

  nqClose_netcdfFile(netcdfId);
  
  /******************************************/
  /* Add some flags to the VisuData object. */
  /******************************************/
  flag = g_malloc(sizeof(gboolean));
  *flag = TRUE;
  visuDataSet_property(data, SCALAR_FIELD_DEFINED_IN_STRUCT_FILE, (gpointer)flag);

  return TRUE;
}