/* AAK:
 * Utility to merge sdf files generated by PAMR
 * Assumes that the sdf files will cover a box region
 */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <bbhutil.h>
#include <sdfmerge.h>
#include <iv.h>
#include <dv.h>
#include <math.h>

#define max_sdf_file 1024
#define max_dim 3
#define max_fname_size 256
#define max_name_size 64
#define m_inf -1.0E25
#define p_inf 1.0E25

double round(double x);


struct my_SDF {
   char file[max_fname_size];
   char name[max_name_size];
   char **cnames;
   int dim;
   int size;
   int coord_size;
   int shape[max_dim];
   double bbox[max_dim*2];
   double *data;
   double *coords;
   int position[max_dim];
} ;

typedef struct my_SDF SDF;

int main(int argc, char *argv[]) {

   SDF gf[max_sdf_file];
   int num_sdf_files;
   int i, j, k, l;
   int lev;
   double time;
   int crd_pty;
   int debug = 0;
   int check_compatibility = 1;
   double bbox[max_dim*2];
   int shape[max_dim];
   double dx[max_dim];
   int size;
   double *data;
   int dim;
   int rc;
   int count;
   int trace = 1;


   num_sdf_files = argc;
   lev = 1;

   for(i=0;i<max_dim;i++) {
      bbox[2*i] = p_inf;
      bbox[2*i+1] = m_inf;
      shape[i] = 0;
   }

   for(i=1;i<num_sdf_files;i++) {

      strcpy(gf[i].file,argv[i]);

      printf("reading rank\n");
      gft_read_rank(gf[i].file,lev,&gf[i].dim);
      printf("done!\n");

      ivls(gf[i].shape,1,max_dim);

      printf("reading shape\n");
      gft_read_shape(gf[i].file,lev,gf[i].shape);

      printf("reading name\n");
      gft_read_name(gf[i].file,1,gf[i].name);

      gf[i].coord_size = 0;
      gf[i].size = 1;

      printf("some calc...\n");
      for( j=0; j<gf[i].dim; j++ ) {
         gf[i].coord_size += gf[i].shape[j];
         gf[i].size *= gf[i].shape[j];
      }

      printf("allocating memory\n");
      gf[i].coords = vec_alloc(gf[i].coord_size);
      gf[i].data = vec_alloc(gf[i].size);

      gf[i].cnames=(char **)malloc(1*sizeof(char *));
      gf[i].cnames[0] = (char *)malloc(6*sizeof(char));


      if (trace == 1) {
        printf("reading %s file...\n",gf[i].file);
      }
      gft_read_full( gf[i].file,lev,gf[i].shape,gf[i].cnames[0],gf[i].dim,
            &time, gf[i].coords, gf[i].data);
      if (trace == 1) {
        printf("done!\n");
      }

      gf[i].bbox[0] = gf[i].coords[0];
      gf[i].bbox[1] = gf[i].coords[gf[i].shape[0]-1];

      crd_pty = 0;

      for(j=0;j<gf[i].dim;j++) {
         gf[i].bbox[j*2] =  gf[i].coords[crd_pty];
         gf[i].bbox[j*2+1] = gf[i].coords[crd_pty+gf[i].shape[j]-1];
         crd_pty += gf[i].shape[j];
      }

      for(j=0;j<gf[1].dim;j++) {
         if (gf[i].bbox[2*j] < bbox[2*j]) { bbox[2*j] = gf[i].bbox[2*j]; }
         if (gf[i].bbox[2*j+1] > bbox[2*j+1] ) { bbox[2*j+1] = gf[i].bbox[2*j+1]; }
      }


   if (trace == 1) {
   }
   } //End of loop over all grid functions


   size = 1;
   dim = gf[1].dim;

   ivls(shape,1,max_dim);
   for(i=0;i<dim;i++) {
      dx[i] = (gf[1].bbox[2*i+1]-gf[1].bbox[2*i] ) / (gf[1].shape[i] - 1);
      shape[i] = (int)round(( bbox[2*i+1] - bbox[2*i] ) / dx[i] ) + 1;
      size *= shape[i];
   }

   data = vec_alloc(size);
   dvls(data,0.0,size);

   for (i=1;i<num_sdf_files;i++) {
      ivls(gf[i].position,0,max_dim);
      for(j=0;j<dim;j++) {
         gf[i].position[j] = (int)round( (gf[i].bbox[2*j] - bbox[2*j]) / dx[j] );
      }
   } //End of loop over all grid functions

   lev = 1;
   rc = 1;
   while (rc == 1) {

      count = 0;
      for (i=1;i<num_sdf_files;i++){
           rc = gft_read_full( gf[i].file,lev,gf[i].shape,gf[i].cnames[0],gf[i].dim,
            &time, gf[i].coords, gf[i].data);
           if (rc == 1) {
             count++;
             inj_grid_func_(data,gf[i].data,shape,gf[i].shape,gf[i].position,&dim);
           }
      }

      if ((rc == 1) && (count !=num_sdf_files-1)) {
         printf("Could not read one of the sdf files at lev:%d\n",lev);
         exit (1);
      }
      else if (rc == 1) {
         if (trace == 1) {
           printf("Output at time %f, lev %d:\n",time,lev);
         }
         gft_out_bbox("merged.sdf",time,shape,dim,bbox,data);
         lev++;
      }

   } //END OF WHILE POSSIBLE TO READ

   printf("sdf files merged to > merged.sdf upto level: %d\n",lev-1);

   gft_close_all();

} //END OF MAIN