/*******************************************************************************
 * code reads data dumped from the vulcan neutrino-hydrodynamics
 *  simulation code for further purpose:
 *  - graphics
 *  - analysis
 *******************************************************************************/

#define PRIME						    /* setup to read omportSnData.h file */

#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/mman.h>
#include <netinet/in.h>
#include <dx/dx.h>
#include <dx/modflags.h>
#include <importSnData.h>

extern Error DXAddModule(char *, ...);
extern Error m_ImportVulcanData3D(Object *, Object *);

#define NDIM_READ 2
#define SET 1
#define TWO_DIM 2
#define THREE_DIM 3

#define NUM_IN_PORTS 7
#define NUM_OUT_PORTS 2

/****************************************************************************/

void DXEntry()
{
  DXAddModule("ImportVilcanData3D", m_ImportVulcanData3D, 0,
	      NUM_IN_PORTS, "file_name", "band_select", "band_select",
	      "view_dimesnion", "segments", "PiFactor", "velocity_dimension",
	      NUM_OUT_PORTS, "sn_data", "time");
}

/****************************************************************************/

Error m_ImportVulcanData3D(Object *in, Object *out)
{
  Field f=NULL;
  Array parray=NULL;
  Array darray=NULL;
  Array varray=NULL;
  Array carray=NULL;
  
  u32 i, j;
  FILE *fp;
  char *filename;

  i32 offset;
  i32 which_scalar, which_vector;
  char InString[255];
  i32 StreamLength;
  i32 StreamLengthCheck;
  char jobLabel[255];
  i64 ntime;
  i32 numelm,npoints,nus_sn,nfreq;
  f64 time; 
  i32 NumbScalarVariables,NumbVectorVariables,NumbTensorVariables;	
  i32 nzmax;
  i32 iii;
  i32 dimension,segments;
  f32 pi_factor;
  i32 NUM_DIM,cnt;
  f64 deg;
  i32 *nodes, **nump;
  f32 **pData, *pptr;
  f32 *pDataR;
  f64 *ScalarVariables,*VarScalarMinMax;
  f64 *ScalarVariablesDummy,*VarScalarMinMaxDummy;
  i32 *NodeCount;
  f64 *VectorVariables,*VarVectorMinMax;
  f64 **VectorVariablesDummy,**VarVectorMinMaxDummy; 
  f64 *vptr;
  i32 nc, nv, inode, ibreak, vel_vec_dim;
  Cube cb;
  Quadrilateral qd;
  
  for(i = 0; i < NUM_OUT_PORTS; i++){			    /* Initialize all outputs to NULL */
    out[i] = NULL;
  }
  if(in[0] == NULL) {					    /* Error checks: required inputs are verified. */
    DXSetError(ERROR_MISSING_DATA, "\"DATA_*\" must be specified");
    return ERROR;
  }
  else if(!DXExtractString(in[0], &filename)) {
    DXSetError(ERROR_BAD_PARAMETER, "filename must be a string");
    return ERROR;
  }
  if(!(fp = fopen(filename, "r"))) {			    /* Check to see that the file is accessible */
    DXSetError(ERROR_BAD_PARAMETER, "file \"%s\" is not accessible", filename);
    return ERROR;
  }
  if(in[1] == NULL){
    which_scalar = 0;
  }
  else{
    DXExtractInteger(in[1], &which_scalar);
  }
  if(in[2] == NULL){
    which_vector = 0;
  }
  else{
    DXExtractInteger(in[2], &which_vector);
  }
  if(in[3] == NULL){
    dimension = TWO_DIM;
    segments = 1;
    NUM_DIM = 2;
  }
  else {
    DXExtractInteger(in[3], &dimension);
  }
  if(dimension == TWO_DIM) {
    segments = 1;
    NUM_DIM = 2;
  }
  else if(dimension == THREE_DIM) {
    NUM_DIM = 3;
    if(in[4] == NULL){
      segments = 36;
    }
    else{
      DXExtractInteger(in[4], &segments);
    }
    if(in[5] == NULL){
      pi_factor = 2;
    }
    else{
      DXExtractFloat(in[5], &pi_factor);
    }
  }
  if (in[6] == NULL){
    vel_vec_dim = 3;
  }
  else{ 
    DXExtractInteger(in[6], &vel_vec_dim);
  }
  
  /* 
   * Since output "object" is structure Field/Group, it initially
   * is a copy of input "View_data".
   */
  out[0] = DXCopy(in[0], COPY_STRUCTURE);
  if(!out[0]){
    goto error;
  }
  /* 
   * If in[0] was an array, then no copy is actually made - Copy
   * returns a pointer to the input object.  Since this can't be written to
   * we postpone explicitly copying it until the leaf level, when we'll need
   * to be creating writable arrays anyway.
   */
  if(out[0] == in[0]){
    out[0] = NULL;
  }
  /* 
   * Initialize all outputs to NULL
   */
  for(i = 0; i < NUM_OUT_PORTS; i++){
    out[i] = NULL;
  }
  /* 
   * START actually read the file and FILL opendx-structure 
   * Check header and binary format 
   * Switch to change endiness should be introduced later
   */
  fread(&StreamLength, sizeof(i32), 1, fp);
  fread(&InString, sizeof(char), StreamLength, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("HEADER START string read incorrectly\n");
    printf("StreamLength = %d \n", StreamLength);
    printf("StreamLength = %d \n", StreamLengthCheck);
  }

  fread(&StreamLength, sizeof(i32), 1, fp);		    /* Check version */
  fread(&InString, sizeof(char), StreamLength, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("VERSION read incorrectly\n");
    printf("StreamLength = %d \n", StreamLength);
    printf("StreamLength = %d \n", StreamLengthCheck);
  }
  /* 
   * Read the Job Label 
   */
  fread(&StreamLength, sizeof(i32), 1, fp);
  fread(jobLabel, sizeof(char), StreamLength, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("Filename read incorrectly\n");
    printf("StreamLength = %d \n", StreamLength);
    printf("StreamLength = %d \n", StreamLengthCheck);
  }
  /* 
   * Read simulation time steps: ntime
   * Number of quads: numelm
   * Number of points: npoints
   * Don't know what: nus_sn
   * Number of Frequencies: nfreq 
   */
  fread(&StreamLength, sizeof(i32), 1, fp);
  fread(&ntime, sizeof(ntime), 1, fp);
  fread(&numelm, sizeof(i32), 1, fp);
  fread(&npoints, sizeof(i32), 1, fp);
  fread(&nfreq, sizeof(i32), 1, fp);
  fread(&nus_sn, sizeof(i32), 1, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("Integer parameter read incorrectly\n");
    printf("StreamLength = %d \n", StreamLength);
    printf("StreamLength = %d \n", StreamLengthCheck);
  }

  /* 
   * Read Time: time
   */
  fread(&StreamLength, sizeof(i32), 1, fp);
  fread(&time, sizeof(f64), 1, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("TIME read incorrectly \n");
    printf("StreamLength = %d \n", StreamLength);
    printf("StreamLength = %d \n", StreamLengthCheck);
  }

  /* 
   * Output "time" is a value; set up an appropriate array
   */
  out[1] = (Object) DXNewArray(TYPE_DOUBLE, CATEGORY_REAL, 0, 0);
  if(!out[1])
    goto error;
  if(!DXAddArrayData((Array) out[1], 0, 1, &time))
    goto error;

  /* 
   * Read Number of Scalar Variables: NumbScalarVariables
   * Number of Flux Variables: NumbVectorVariables 
   * Number of Tensor Variables: NumbTensorVariables 
   */
  fread(&StreamLength, sizeof(i32), 1, fp);
  fread(&NumbScalarVariables, sizeof(i32), 1, fp);
  fread(&NumbVectorVariables, sizeof(i32), 1, fp);
  fread(&NumbTensorVariables, sizeof(i32), 1, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("Number of variables read incorrectly\n");
    printf("StreamLength=%d\n", StreamLength);
    printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
  }

  /* Create a new field. */
  f = DXNewField();
  if(!f)
    goto error;

  /* 
   * Read nzmax
   */
  fread(&StreamLength, sizeof(i32), 1, fp);
  fread(&nzmax, sizeof(i32), 1, fp);
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("NZMAX incorrectly read\n");
    printf("StreamLength=%d\n", StreamLength);
    printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
  }

  /* 
   * Read number of corners for each cell: nodes
   */
  fread(&StreamLength, sizeof(i32), 1, fp);

  nodes = (i32 *)calloc(numelm, sizeof(i32));
  for(i = 0; i < numelm; i++) {
    fread(nodes + i, sizeof(i32), 1, fp);
  }
  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("\nNODES read incorrectly\n");
    printf("StreamLength=%d\n", StreamLength);
    printf("StreamLengthCheck=%d\n", StreamLengthCheck);
  }

  /* 
   * Read Connections
   */

  fread(&StreamLength, sizeof(i32), 1, fp);

  nump = (i32 **)calloc(numelm, sizeof(i32 *));
  for(i = 0; i < numelm; i++) {
    nump[i] = (i32 *)calloc(nodes[i], sizeof(i32));
    for(j = 0; j < nodes[i]; j++)
      fread(nump[i] + j, sizeof(i32), 1, fp);
  }

  fread(&StreamLengthCheck, sizeof(i32), 1, fp);

  if(StreamLength != StreamLengthCheck) {
    printf("\nCONNECTIONS read incorrectly\n");
    printf("StreamLength=%d\n", StreamLength);
    printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
  }

  /* 
   * Read co-ordinates: Xp, Yp (Same for 2D and 3D)
   */
  pData = (f32 **)calloc(NDIM_READ, sizeof(f32 *));

  for(i = 0; i < NDIM_READ; i++) {

    fread(&StreamLength, sizeof(i32), 1, fp);

    pData[i] = (f32 *)calloc(npoints, sizeof(f32));
    for(j = 0; j < npoints; j++)
      fread(pData[i] + j, sizeof(f32), 1, fp);

    fread(&StreamLengthCheck, sizeof(i32), 1, fp);
  }

  if(StreamLength != StreamLengthCheck) {
    printf("\nCO-ORDINATES read incorrectly\n");
    printf("StreamLength=%d\n", StreamLength);
    printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
  }
  fflush(stderr);
  fflush(stdout);

  pDataR = (f32 *)calloc(3 * segments * npoints, sizeof(f32));
  pptr = pDataR;
  inode = 0;

  if(dimension == THREE_DIM) {
    for(j = 0; j < npoints; j++) {
      *pptr++ = pData[0][j];
      inode = inode + 1;
      *pptr++ = 0.;
      inode = inode + 1;
      
      *pptr++ = pData[1][j];
      inode = inode + 1;
    }
  }
  else {
    for(j = 0; j < npoints; j++) {
      for(i = 0; i < NDIM_READ; i++) {
	*pptr++ = pData[i][j];
	inode = inode + 1;
      }
      *pptr++ = 0.;
      inode = inode + 1;
    }
  }

  /* Prepare DX Array for 3D View */

  if(dimension == THREE_DIM) {
    deg = (pi_factor * M_PI) / segments;

    for(j = 1; j < segments; j++) {
      for(i = 0; i < npoints; i++) {
	*pptr++ = cos(j * deg) * pDataR[3 * i];
	*pptr++ = -sin(j * deg) * pDataR[3 * i];
	*pptr++ = pDataR[3 * i + 2];
      }
    }

  }

  /* assign to dx-data-structure */
  parray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
  if(!parray)
    goto error;

  if(!DXAddArrayData(parray, 0, npoints * segments, NULL))
    goto error;

  if(!DXAddArrayData(parray, 0, npoints * segments, pDataR))
    goto error;

  if(!DXSetComponentValue(f, "positions", (Object) parray))
    goto error;


  for(i = 0; i < NDIM_READ; i++)
    free(pData[i]);

  free(pData);
  free(pDataR);

  /* assign quads */

  if(NUM_DIM == 2) {					    /* If 2D, its a quad */
    for(i = 0; i < numelm; i++) {
      qd =
	DXQuad(*(nump[i]) - 1, *(nump[i] + 1) - 1, *(nump[i] + 3) - 1,
	*(nump[i] + 2) - 1);
      f = DXAddQuad(f, i, qd);
    }
  }
  else {						    /* If 3D, its a cube */
    cnt = 0;
    carray = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 8);
    if(!carray)
      goto error;
    DXSetComponentValue(f, "connections", (Object) carray);
    DXSetAttribute((Object) carray, "element type",
      (Object) DXNewString("cubes"));
    for(i = 0; i < numelm; i++) {
      for(j = 0; j < segments - 1; j++) {
	cb.p = *(nump[i]) + (j * npoints) - 1;		    /* Current Frame */
	cb.q = *(nump[i] + 1) + (j * npoints) - 1;
	cb.r = *(nump[i] + 3) + (j * npoints) - 1;
	cb.s = *(nump[i] + 2) + (j * npoints) - 1;
	cb.t = *(nump[i]) + ((j + 1) * npoints) - 1;	    /* Next Frame */
	cb.u = *(nump[i] + 1) + ((j + 1) * npoints) - 1;
	cb.v = *(nump[i] + 3) + ((j + 1) * npoints) - 1;
	cb.w = *(nump[i] + 2) + ((j + 1) * npoints) - 1;
	if(!DXAddArrayData(carray, cnt++, 1, &cb))
	  goto error;
      }
      if(pi_factor == 2) {				    /* Last frame joined to first frame only if pi_factor = 2 */
	cb.p = *(nump[i]) + ((segments - 1) * npoints) - 1; /* Current Frame */
	cb.q = *(nump[i] + 1) + ((segments - 1) * npoints) - 1;
	cb.r = *(nump[i] + 3) + ((segments - 1) * npoints) - 1;
	cb.s = *(nump[i] + 2) + ((segments - 1) * npoints) - 1;
	cb.t = *(nump[i]) - 1;				    /* Next Frame */
	cb.u = *(nump[i] + 1) - 1;
	cb.v = *(nump[i] + 3) - 1;
	cb.w = *(nump[i] + 2) - 1;
	if(!DXAddArrayData(carray, cnt++, 1, &cb))
	  goto error;
      }
    }
  }

  /************************************************************************************
   *                                                                                  *
   * Read and assign ScalarVariables                                                  *
   *                                                                                  *
   ************************************************************************************/

  ScalarVariables = (f64 *)calloc(npoints, sizeof(f64));
  VarScalarMinMax = (f64 *)calloc(2, sizeof(f64));
  NodeCount = (i32 *)calloc(npoints, sizeof(i32));
  ScalarVariablesDummy = (f64 *)calloc(numelm, sizeof(f64));
  VarScalarMinMaxDummy = (f64 *)calloc(2, sizeof(f64));

  for(i = 0; i < NumbScalarVariables; i++) {
    fread(&StreamLength, sizeof(i32), 1, fp);		    /* Read Variable */
    for(j = 0; j < 2; j++){
      fread(VarScalarMinMaxDummy + j, sizeof(f64), 1, fp);
    }
    for(j = 0; j < numelm; j++){
      fread(ScalarVariablesDummy + j, sizeof(f64), 1, fp);
    }
    fread(&StreamLengthCheck, sizeof(i32), 1, fp);

    if(StreamLength != StreamLengthCheck) {		    /* Check correct reading */
      printf("\nScalar variable = %d read incorrectly\n", i);
      printf("StreamLength=%d\n", StreamLength);
      printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
    }
    if(i == which_scalar) {				    /* variable selector */
      for(j = 0; j < 2; j++){
	VarScalarMinMax[j] = VarScalarMinMaxDummy[j];
      }
      for(j = 0; j < npoints; j++) {
	ScalarVariables[j] = 0.;
	NodeCount[j] = 0;
      }
      for(j = 0; j < numelm; j++) {			    /* distribute cell-centered values to nodes */
	for(nc = 0; nc < nodes[j]; nc++) {
	  inode = nump[j][nc] - 1;
	  ScalarVariables[inode] += ScalarVariablesDummy[j];
	  NodeCount[inode]++;
	}
      }
      for(j = 0; j < npoints; j++){
	ScalarVariables[j] = ScalarVariables[j] / NodeCount[j];
      }
      darray = DXNewArray(TYPE_DOUBLE, CATEGORY_REAL, 1, 1); /* assign variable to dx-array */
      if(!darray){
	goto error;
      }
      if(dimension == TWO_DIM) {
	if(!DXAddArrayData(darray, 0, npoints, ScalarVariables)){
	  goto error;
	}
	if(!DXSetComponentValue(f, "scalar_data", (Object) darray)){
	  goto error;
	}
      }
      else {
	for(i = 0; i < segments; i++){
	  if(!DXAddArrayData(darray, i * npoints, npoints, ScalarVariables)){
	    goto error;
	  }
	}
	if(!DXSetComponentValue(f, "scalar_data", (Object) darray)){
	  goto error;
	}
      }
      if(!DXSetAttribute((Object) darray, "dep", (Object) DXNewString("positions"))){
	goto error;
      }
    }
  }
  free(ScalarVariables);
  free(VarScalarMinMax);
  free(ScalarVariablesDummy);
  free(VarScalarMinMaxDummy);
  free(NodeCount);

  /************************************************************************************
   *                                                                                  *
   * Read and assign VectorVariables                                                  *
   *                                                                                  *
   ************************************************************************************/

  VectorVariables = (f64 *)calloc(3 * segments * npoints, sizeof(f64));
  VarVectorMinMax = (f64 *)calloc(2, sizeof(f64));
  VectorVariablesDummy = (f64 **)calloc(NDIM_READ, sizeof(f64 *));
  VarVectorMinMaxDummy = (f64 **)calloc(NDIM_READ, sizeof(f64 *));

  for(nv = 0; nv < NumbVectorVariables; nv++) {		    /* Read the two components of the vector */
    ibreak = 0;
    iii = NDIM_READ;
    if(nv == 0){
      iii = vel_vec_dim;
    }
    for(i = 0; i < iii; i++) {
      if(ibreak != 0){
	break;
      }
      fread(&StreamLength, sizeof(i32), 1, fp);
      VarVectorMinMaxDummy[i] = (f64 *)calloc(2, sizeof(f64));
      for(j = 0; j < 2; j++){
	fread(VarVectorMinMaxDummy[i] + j, sizeof(f64), 1, fp);
      }
      VectorVariablesDummy[i] = (f64 *)calloc(npoints, sizeof(f64));
      for(j = 0; j < npoints; j++){
	fread(VectorVariablesDummy[i] + j, sizeof(f64), 1, fp);
      }
      fread(&StreamLengthCheck, sizeof(i32), 1, fp);
      if(StreamLength != StreamLengthCheck) {		    /* Check correct reading */
	printf("\nVector variable = %d read incorrectly\n", i);
	printf("StreamLength=%d\n", StreamLength);
	printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
      }
    }
    if(nv == which_vector) {				    /* variable selector */
      for(j = 0; j < 2; j++){
	VarVectorMinMax[j] = *(VarVectorMinMaxDummy[0] + j);
      }
      vptr = VectorVariables;				    /* re-arrange vectors */
      inode = 0;
      if(NUM_DIM == 2) {
	for(j = 0; j < npoints; j++) {
	  for(i = 0; i < NDIM_READ; i++) {
	    *vptr++ = VectorVariablesDummy[i][j];
	    inode = inode + 1;
	  }
	  *vptr++ = 0.;
	  inode = inode + 1;
	}
      }
      else {
	for(j = 0; j < npoints; j++) {
	  *vptr++ = VectorVariablesDummy[0][j];
	  inode = inode + 1;
	  *vptr++ = 0.;
	  inode = inode + 1;
	  *vptr++ = VectorVariablesDummy[1][j];
	  inode = inode + 1;
	}
	deg = (pi_factor * M_PI) / segments;
	for(j = 1; j < segments; j++) {
	  for(i = 0; i < npoints; i++) {
	    *vptr++ = cos(j * deg) * VectorVariables[3 * i];
	    *vptr++ = -sin(j * deg) * VectorVariables[3 * i];
	    *vptr++ = VectorVariables[3 * i + 2];
	  }
	}
      }
    }
    inode = 0;
    for(j = 0; j < npoints; j++) {
      for(i = 0; i < 3; i++) {
	inode = inode + 1;
      }
    }
    varray = DXNewArray(TYPE_DOUBLE, CATEGORY_REAL, 1, 3);  /* assign variable to dx-array */
    if(!darray){
      goto error;
    }
    if(!DXAddArrayData(varray, 0, npoints * segments, VectorVariables)){
      goto error;
    }
    if(!DXSetComponentValue(f, "vector_data", (Object) varray)){
      goto error;
    }
    if(!DXSetAttribute((Object) varray, "dep", (Object) DXNewString("positions"))){
      goto error;
    }
    ibreak = 1;
  }

  offset =  (NumbVectorVariables - 1 - which_vector) * NDIM_READ *
    (2 * sizeof(i32) + (2 + npoints) * sizeof(f64));
  fseek(fp, offset, SEEK_CUR);

  free(VectorVariables);
  free(VarVectorMinMax);
  free(VectorVariablesDummy);
  free(VarVectorMinMaxDummy);
  for(i = 0; i < numelm; i++){				    /* deallocate connections */
    free(nump[i]);
  }
  free(nump);
  free(nodes);

  /* finish up */
  
  if(!DXEndField(f)){					    /* Set field default values & bounding box. */
    goto error;
  }
  out[0] = (Object) f;					    /* Set output to be the created field. */

  parray = NULL;
  darray = NULL;

  return OK;

  /* 
   * On error, any successfully-created outputs are deleted.
   */
error:
  for(i = 0; i < 2; i++) {
    if(in[i] != out[i]){
      DXDelete(out[i]);
    }
    out[i] = NULL;
  }
  return ERROR;
}