/*   EXTRAITS DE LA LICENCE
	Copyright CEA, contributeurs : Damien CALISTE, laboratoire L_Sim, (2001-2008)
  
	Adresse ml :
	CALISTE, damien P caliste AT cea P fr.

	Ce logiciel est un programme informatique servant  visualiser des
	structures atomiques dans un rendu pseudo-3D. 

	Ce logiciel est rgi par la licence CeCILL soumise au droit franais et
	respectant les principes de diffusion des logiciels libres. Vous pouvez
	utiliser, modifier et/ou redistribuer ce programme sous les conditions
	de la licence CeCILL telle que diffuse par le CEA, le CNRS et l'INRIA 
	sur le site "http://www.cecill.info".

	Le fait que vous puissiez accder  cet en-tte signifie que vous avez 
	pris connaissance de la licence CeCILL, et que vous en avez accept les
	termes (cf. le fichier Documentation/licence.fr.txt fourni avec ce logiciel).
*/

/*   LICENCE SUM UP
	Copyright CEA, contributors: Damien CALISTE, L_Sim laboratory, (2001-2008)

	E-mail address:
	CALISTE, damien P caliste AT cea P fr.

	This software is a computer program whose purpose is to visualize atomic
	configurations in 3D.

	This software is governed by the CeCILL  license under French law and
	abiding by the rules of distribution of free software.  You can  use, 
	modify and/ or redistribute the software under the terms of the CeCILL
	license as circulated by CEA, CNRS and INRIA at the following URL
	"http://www.cecill.info". 

	The fact that you are presently reading this means that you have had
	knowledge of the CeCILL license and that you accept its terms. You can
	find a copy of this licence shipped with this software at
        Documentation/licence.en.txt.
*/

#include <visu_basic.h>
#include <visu_rendering.h>
#include <visu_data.h>
#include <renderingMethods/renderingAtomic_ascii.h>
#include <coreTools/toolMatrix.h>
#include <coreTools/toolPhysic.h>
#include <extraFunctions/scalarFields.h>

#define CUBE_DESCRIPTION _("<span size=\"smaller\">"			\
			   "This plug-in reads <b>Cube</b>"		\
			   " files (see \n<u><span color=\"blue\">"	\
			   "http://local.wasp.uwa.edu.au/\n"		\
			   "~pbourke/dataformats/cube/"			\
			   "</span></u>),\nboth structural and"		\
			   " and volumetric data.</span>")
#define CUBE_AUTHORS     "Caliste Damien"

/* http://local.wasp.uwa.edu.au/~pbourke/dataformats/cube/ */

struct cube_reader
{
  const gchar  *filename;
  GString      *line;
  gchar        *comment;
  double        box[3][3];
  float         shift[3];
  int           ntype, natom;
  GHashTable   *elements;
  VisuElement **nodeTypes;
  float        *coords;
  int           mesh[3];
  double       *density;
  GIOStatus     status;
  GIOChannel   *flux;
};

static gboolean cubeAtomicLoad(VisuData *data, const gchar* filename,
			       ToolFileFormat *format, int nSet, GError **error);
static gboolean cubeDensityLoad(const gchar *filename, GList **fieldList,
				GError **error, GHashTable *table);
/* static gboolean xvAtomicLoad(VisuData *data, const gchar* filename, */
/* 			     ToolFileFormat *format, int nSet, GError **error); */

/* Local variables */
static gchar *iconPath;

gboolean cubeInit()
{
  char *typeCube[] = {"*.cube", "*.cub", (char*)0};
  char *nameCube   = "Cube format (from Gaussian)";
/*   char *typeXV[] = {"*.XV", "*.STRUCT_OUT", (char*)0}; */
/*   char *nameXV   = "Siesta output format"; */
/*   char *descrXV  = _("Siesta geometric output format"); */
  ToolFileFormat *fmt;

  fmt = tool_file_format_new(_("Gaussian structural/volumetric format"),
                             typeCube);
  visu_rendering_addFileFormat(visu_rendering_getByName(VISU_RENDERING_ATOMIC_NAME), 0,
                               fmt, 60, cubeAtomicLoad, (GDestroyNotify)0);

  scalarFieldAdd_loadMethod(nameCube, cubeDensityLoad,
			    fmt, G_PRIORITY_HIGH - 5);

/*   fmt = tool_file_format_new(descrXV, typeXV); */

/*   meth = g_malloc(sizeof(VisuRenderingFormatLoad)); */
/*   meth->name = nameXV; */
/*   meth->fmt = fmt; */
/*   meth->priority = 61; */
/*   meth->load = xvAtomicLoad; */
/*   visu_renderingAtomic_addLoadMethod(meth); */

  iconPath = g_build_filename(V_SIM_PIXMAPS_DIR, "cube.png", NULL);

  return TRUE;
}
const char* cubeGet_description()
{
  return CUBE_DESCRIPTION;
}

const char* cubeGet_authors()
{
  return CUBE_AUTHORS;
}

const char* cubeGet_icon()
{
  return iconPath;
}


static void reader_free(struct cube_reader *rd)
{
  DBG_fprintf(stderr, "Cube: free the reader structure.\n");

  g_string_free(rd->line, TRUE);
  if (rd->comment)
    g_free(rd->comment);
  g_hash_table_destroy(rd->elements);
  if (rd->nodeTypes)
    g_free(rd->nodeTypes);
  if (rd->coords)
    g_free(rd->coords);
  if (rd->density)
    g_free(rd->density);
  
  if (g_io_channel_shutdown(rd->flux, FALSE, (GError**)0) != G_IO_STATUS_NORMAL)
    g_warning("Cube: can't close file.");
  g_io_channel_unref(rd->flux);
}

static gboolean read_line(struct cube_reader *rd, GError **error)
{
  rd->status = g_io_channel_read_line_string(rd->flux, rd->line, NULL, error);
  return (rd->status == G_IO_STATUS_NORMAL);
}

static gboolean read_comment(struct cube_reader *rd, GError **error)
{
  /* Read the two comment line. */
  if (!read_line(rd, error))
    {
      *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			   _("Wrong Cube format, missing comment lines.\n"));
      return FALSE;
    }
  g_strstrip(rd->line->str);
  rd->comment = g_strdup(rd->line->str);
  if (!read_line(rd, error))
    {
      *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			   _("Wrong Cube format, missing comment lines.\n"));
      return FALSE;
    }
  g_strstrip(rd->line->str);
  g_string_prepend(rd->line, " ");
  g_string_prepend(rd->line, rd->comment);
  g_free(rd->comment);
  rd->comment = g_strdup(rd->line->str);
  DBG_fprintf(stderr, "Cube: comment '%s'.\n", rd->comment);
  return TRUE;
}

static gboolean read_natom(struct cube_reader *rd, GError **error)
{
  /* Read the number of atoms. */
  if (!read_line(rd, error))
    {
      *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			   _("Wrong Cube format, missing the atom line.\n"));
      return FALSE;
    }
  if (sscanf(rd->line->str, "%d %f %f %f", &rd->natom,
	     rd->shift, rd->shift + 1, rd->shift + 2) != 4)
    {
      *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			   _("Wrong Cube format, missing the number of"
			     " atoms on the third line.\n"));
      return FALSE;
    }
  rd->natom = ABS(rd->natom);
  DBG_fprintf(stderr, "Cube: number of atoms %d.\n", rd->natom);
  return TRUE;
}

static gboolean read_box(struct cube_reader *rd, GError **error)
{
  int i;
  gchar dir[3] = "XYZ";

  /* Read the box definition. */
  for (i = 0; i < 3; i++)
    {
      if (!read_line(rd, error))
	{
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			       _("Wrong Cube format, missing the %d box line.\n"),
			       i + 1);
	  return FALSE;
	}
      if (sscanf(rd->line->str, "%d %lf %lf %lf\n",
		 rd->mesh + i, rd->box[i], rd->box[i] + 1, rd->box[i] + 2) != 4)
	{
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			       _("Wrong Cube format, missing float values"
				 " for box definition.\n"));
	  return FALSE;
	}
      if (rd->mesh[i] <= 0)
	{
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			       _("Wrong Cube format, wrong mesh"
				 " size in %c direction.\n"), dir[i]);
	  return FALSE;
	}
      rd->box[i][0] *= rd->mesh[i];
      rd->box[i][1] *= rd->mesh[i];
      rd->box[i][2] *= rd->mesh[i];
    }
  DBG_fprintf(stderr, "Cube: read box %8g %8g %8g\n"
	      "               %8g %8g %8g\n"
	      "               %8g %8g %8g\n",
	      rd->box[0][0], rd->box[0][1], rd->box[0][2],
	      rd->box[1][0], rd->box[1][1], rd->box[1][2],
	      rd->box[2][0], rd->box[2][1], rd->box[2][2]);
  return TRUE;
}

static gboolean read_coords(struct cube_reader *rd, GError **error)
{
  int i, zele, nb;
  float rcov, trash;
  gchar *name, *ptChar;
  VisuElement *type;
  struct VisuMethAscii *infos;
  gboolean new;

  rd->ntype     = 0;
  rd->nodeTypes = g_malloc(sizeof(VisuElement*) * rd->natom);
  rd->coords    = g_malloc(sizeof(float) * 3 * rd->natom);

  DBG_fprintf(stderr, "Cube: read coordinates.\n");
  for  (i = 0; i < rd->natom; i++)
    {
      if (!read_line(rd, error))
	{
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			       _("Wrong Cube format, missing"
				 " line %d with coordinates.\n"),
			       i + 1);
	  return FALSE;
	}
      nb = sscanf(rd->line->str, "%d %f %f %f %f\n", &zele, &trash,
		  rd->coords + 3 * i, rd->coords + 3 * i + 1, rd->coords + 3 * i + 2);
      if (nb != 5)
	{
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			       _("Wrong Cube format, can't read coordinates.\n"));
	  return FALSE;
	}
      rd->coords[3 * i + 0] -= rd->shift[0];
      rd->coords[3 * i + 1] -= rd->shift[1];
      rd->coords[3 * i + 2] -= rd->shift[2];
      /* Try to find a z number instead of a name. */
      rcov = -1.f;
      if (!tool_physic_getElementFromZ(&ptChar, &rcov, zele))
        {
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			       _("Wrong Cube format, one atomic number is 0 or lower,"
                                 " or greater than 103, check your input file.\n"));
	  return FALSE;
        }
      name = g_strdup(ptChar);
      /* adding nomloc to the hashtable */
      type = visu_element_retrieveFromName(name, &new);
      if (new && rcov > 0.f)
	visu_renderingAtomic_setRadius(type, rcov * 0.52917720859);
      g_return_val_if_fail(type, FALSE);
      rd->nodeTypes[i] = type;
      infos = (struct VisuMethAscii*)g_hash_table_lookup(rd->elements,
						     (gconstpointer)type);
      if (!infos)
	{
	  infos = g_malloc(sizeof(struct VisuMethAscii));
	  infos->ele = type;
	  infos->pos = rd->ntype;
	  infos->nbNodes = 1;
	  g_hash_table_insert(rd->elements, (gpointer)type, (gpointer)infos);
	  rd->ntype += 1;
	}
      else
	infos->nbNodes += 1;
    }
  DBG_fprintf(stderr, " | read OK.\n");
  return TRUE;
}

static gboolean read_density(struct cube_reader *rd, GError **error)
{
  int i, j, n;
  gchar **tokens;
#if DEBUG == 1
  GTimer *timer;
  gulong fractionTimer;
#endif

#if DEBUG == 1
  timer = g_timer_new();
  g_timer_start(timer);
#endif

  n = rd->mesh[0] * rd->mesh[1] * rd->mesh[2];
  rd->density = g_malloc(sizeof(double) * n);

  DBG_fprintf(stderr, "Cube: read density (%gMo).\n",
              sizeof(double) * n / 1024. / 1024.);
  i = 0;
  do
    {
      if (!read_line(rd, error))
	{
	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
	        	       _("Wrong Cube format, missing density lines.\n"));
	  return FALSE;
	}
      tokens = g_strsplit(rd->line->str, " ", 0);
      for (j = 0; tokens[j] && i < n; j++)
        if (tokens[j][0])
          {
            if (sscanf(tokens[j], "%lf", rd->density + i) != 1)
              {
                *error = g_error_new(VISU_ERROR_RENDERING,
                                     RENDERING_ERROR_FORMAT,
                                     _("Wrong Cube format, unreadable float value"
                                       " %d for density.\n"), i);
                g_strfreev(tokens);
                return FALSE;
              }
            i += 1;
          }
      g_strfreev(tokens);
    }
  while(i < n && rd->status == G_IO_STATUS_NORMAL);
#if DEBUG == 1
  g_timer_stop(timer);
  fprintf(stderr, "Cube: density parsed in %g milli-s.\n", g_timer_elapsed(timer, &fractionTimer)/1e-3);
  g_timer_destroy(timer);
#endif
  if (rd->status != G_IO_STATUS_NORMAL)
    {
      *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT,
			   _("Wrong Cube format, missing float values"
			     " for density (%d read, %d awaited).\n"),
			   i, rd->mesh[0] * rd->mesh[1] * rd->mesh[2]);
      return FALSE;
    }
  return TRUE;
}

static int read_cube_file(struct cube_reader *rd, GError **error,
			  gboolean withDensity)
{
  /* We read every line that corresponds to this schema : "%s %f %f %f" */
  DBG_fprintf(stderr, "Cube: reading file as a cube file.\n");

  /* The storage for read line. */
  rd->line      = g_string_new("");
  rd->comment   = (gchar*)0;

  /* Storage of number of elements per types. */
  rd->ntype     = 0;
  rd->elements  = g_hash_table_new_full(g_direct_hash, g_direct_equal, NULL, g_free);
  rd->nodeTypes = (VisuElement**)0;

  rd->coords    = (float*)0;
  rd->density   = (double*)0;

  /* We read the all the file. */
  if (!read_comment(rd, error))
    {
      DBG_fprintf(stderr, "Cube: no comment, not a cube file.\n");
      return -1;
    }
  if (!read_natom(rd, error))
    return -1;
  if (!read_box(rd, error))
    return -1;

  if (!read_coords(rd, error))
    return 1;
  if (withDensity)
    {
      if (!read_density(rd, error))
	return 1;
    }
  return 0;
}

static gboolean cubeAtomicLoad(VisuData *data, const gchar* filename,
			       ToolFileFormat *format _U_, int nSet _U_, GError **error)
{
  int res, iNodes, i;
  struct cube_reader rd;
  VisuElement **types;
  unsigned int *nattyp;
  gboolean *flag;

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

  rd.filename = filename;
  rd.flux = g_io_channel_new_file(filename, "r", error);
  if (!rd.flux)
    return FALSE;

  res = read_cube_file(&rd, error, FALSE);

  if (res == 0)
    {
      /* Allocate the space for the nodes. */
      types  = g_malloc(sizeof(VisuElement*) * rd.ntype);
      nattyp = g_malloc(sizeof(int) * rd.ntype);
      g_hash_table_foreach(rd.elements, (GHFunc)visu_meth_asciiValToType, (gpointer)&types);
      g_hash_table_foreach(rd.elements, (GHFunc)visu_meth_asciiValToNb, (gpointer)&nattyp);

      DBG_fprintf(stderr, " | begin to transfer data to VisuData.\n");
      /* Begin the storage into VisuData. */
      visu_data_setNSubset(data, 1);

      res = visu_data_setPopulation(data, rd.ntype, nattyp, types);
      if (!res)
	{
	  g_error("Can't store the nodes in the VisuData object.");
	}
      DBG_fprintf(stderr, "Cube: there are %d types in this file.\n", rd.ntype);
      if (DEBUG)
	for (i = 0; i < rd.ntype; i++)
	  fprintf(stderr, " | %d atom(s) for type %d.\n", nattyp[i], i);
      g_free(nattyp);
      g_free(types);
         
      /* Store the coordinates */
      visu_data_setBoxGeometryFull(data, rd.box, VISU_DATA_BOX_PERIODIC);
      for(iNodes = 0; iNodes < rd.natom; iNodes++)
	visu_data_addNodeFromElement(data, rd.nodeTypes[iNodes],
                                     rd.coords + 3 * iNodes, FALSE, FALSE);
      visu_data_applyBoxGeometry(data);
      visu_data_setFileCommentary(data, rd.comment, 0);

      flag = g_malloc(sizeof(gboolean));
      *flag = TRUE;
      g_object_set_data_full(G_OBJECT(data), SCALAR_FIELD_DEFINED_IN_STRUCT_FILE,
			     (gpointer)flag, g_free);
    }

  /* Free the local data. */
  reader_free(&rd);  

  if (res < 0)
    /* The file is not a Cube file. */
    return FALSE;
  else if (res > 0)
    /* The file is a Cube file but some errors occured. */
    return TRUE;
  /* Everything is OK. */
  *error = (GError*)0;
  return TRUE;
}

static gboolean cubeDensityLoad(const gchar *filename, GList **fieldList,
				GError **error, GHashTable *table _U_)
{
  int res;
  struct cube_reader rd;
  ScalarField *field;

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

  rd.filename = filename;
  rd.flux = g_io_channel_new_file(filename, "r", error);
  if (!rd.flux)
    return FALSE;

  res = read_cube_file(&rd, error, TRUE);

  if (res == 0)
    {
      field = scalarFieldNew(rd.filename);
      if (!field)
	g_warning("impossible to create a ScalarField object.");
      else
	{
	  scalarFieldSet_periodic(field, TRUE);
	  scalarFieldSet_commentary(field, rd.comment);
	  scalarFieldSet_boxFull(field, rd.box);
	  scalarFieldSet_gridSize(field, rd.mesh);
	  DBG_fprintf(stderr, "Cube: transfer density into field object.\n");
	  scalarFieldSet_data(field, rd.density, FALSE);
	  *fieldList = g_list_append(*fieldList, (gpointer)field);
	}
    }

  /* Free the local data. */
  reader_free(&rd);  

  if (res < 0)
    /* The file is not a Cube file. */
    return FALSE;
  else if (res > 0)
    /* The file is a Cube file but some errors occured. */
    return TRUE;
  /* Everything is OK. */
  *error = (GError*)0;
  return TRUE;
}

/* static gboolean xvAtomicLoad(VisuData *data, const gchar* filename, */
/* 			     ToolFileFormat *format, int nSet, GError **error) */
/* { */
/*   GIOChannel *flux; */
/*   GIOStatus status; */
/*   GString *line; */
/*   int i, n; */
/*   double cart[3][3]; */
/*   float box[6]; */

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

/*   flux = g_io_channel_new_file(filename, "r", error); */
/*   if (!flux) */
/*     return FALSE; */

/*   for (i = 0; i < 3; i++) */
/*     { */
/*       status = g_io_channel_read_line_string(flux, line, NULL, error); */
/*       if (status != G_IO_STATUS_NORMAL) */
/* 	{ */
/* 	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT, */
/* 			       _("Wrong Siesta format, missing the %d box line.\n"), */
/* 			       i + 1); */
/* 	  return FALSE; */
/* 	} */
/*       g_strstrip(line->str); */

/*       if (sscanf(line->str, "%lf %lf %lf\n", */
/* 		 cart[i], cart[i] + 1, cart[i] + 2) != 4) */
/* 	{ */
/* 	  *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT, */
/* 			       _("Wrong Siesta format, missing float values" */
/* 				 " for box definition.\n")); */
/* 	  return FALSE; */
/* 	} */
/*     } */
/*   DBG_fprintf(stderr, "Siesta: read box %8g %8g %8g\n" */
/* 	      "              %8g %8g %8g\n" */
/* 	      "              %8g %8g %8g\n", */
/* 	      cart[0][0], cart[0][1], cart[0][2], */
/* 	      cart[1][0], cart[1][1], cart[1][2], */
/* 	      cart[2][0], cart[2][1], cart[2][2]); */
/*   status = g_io_channel_read_line_string(flux, line, NULL, error); */
/*   if (status != G_IO_STATUS_NORMAL) */
/*     { */
/*       *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT, */
/* 			   _("Wrong Siesta format, missing the number" */
/* 			     " of elements line.\n")); */
/*       return FALSE; */
/*     } */
/*   g_strstrip(line->str); */
/*   if (sscanf(line->str, "%d\n", &n)) */
/*     { */
/*       *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT, */
/* 			   _("Wrong Siesta format, missing the number" */
/* 			     " of elements.\n")); */
/*       return FALSE; */
/*     } */


/*   if (!tool_matrix_reducePrimitiveVectors(box, cart)) */
/*     { */
/*       *error = g_error_new(VISU_ERROR_RENDERING, RENDERING_ERROR_FORMAT, */
/* 			   _("Wrong Siesta format, primitive vectors are not 3D.\n")); */
/*       return FALSE; */
/*     } */
  


/*   if (g_io_channel_shutdown(flux, FALSE, (GError**)0) != G_IO_STATUS_NORMAL) */
/*     g_warning("Cube: can't close file."); */
/* } */