/*  lanl_ic2gadget  -- convert binary format initial conditions particle files from the 
                       Heitmann et al. 2004/2007 Code Comparison Project 
                       (so-called "LANL" format) to Gadget2 format

Usage:  lanl_ic2gadget lanl_ic_file > ic_gadget_particles.dat 

Example: lanl_ic2gadget particles_ic_lg > ic_lg_gadget.dat

Note:   lanl_ic2gadget is clever enough to figure out if you are converting the 
        large box (L = 256 Mpc/h) IC file or the small box (L = 64 Mpc/h)
        file. 

Warning:   Do not use this to convert the z = 0 particle distributions to Gadget
           format!!! Use lanl_out2gadget instead!!!

More info:  http://t8web.lanl.gov/people/heitmann/arxiv/data.html

Developed from Gadget source C code snippets by Chris Orban (Ohio State University), 10/08    

*/

#include <stdio.h>
#include <string.h>
#include <math.h>

#define NumPart 256*256*256



struct io_header_1
{
  int      npart[6];
  double   mass[6];
  double   time;
  double   redshift;
  int      flag_sfr;
  int      flag_feedback;
  int      npartTotal[6];
  int      flag_cooling;
  int      num_files;
  double   BoxSize;
  double   Omega0;
  double   OmegaLambda;
  double   HubbleParam;
  char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
} header1;


struct particle_data
{
  float  x,vx;
  float  y,vy;
  float  z,vz;
  float  pmass;
  int    ptag;                                                                              
} P[NumPart];


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

   float h = 0.71; // The Heitmann et al. 05 comparison uses Ho = 71 km/s/Mpc
   float sqrta;
   int i;
   size_t count;

   if (argc== 1)
   {
      printf("Tool for converting Heitmann et al. 2004 initial conditions files to Gadget format!!! \n\n");  
      printf("Usage: \n \n");
      printf("           lanl_ic2gadget lanl_ic_file > ic_gadget_particles.dat \n \n");
      printf("Example:   lanl_to_gadget particles_ic_lg > ic_lg_gadget.dat \n\n");
      printf("See lanl_ic2gadget.c header for more details. \n");
      exit(0); 
   }


   FILE *fp = fopen(argv[1], "rb");

   for(i=0;i<NumPart;i++){
     fread(&P[i],sizeof(struct particle_data),1,fp);
     /*    printf("P[%d].x = %5.10f \n",i,P[i].x);
    printf("P[%d].vx = %5.10f \n",i,P[i].vx);
    printf("P[%d].y = %5.10f \n",i,P[i].y);
    printf("P[%d].vy = %5.10f \n",i,P[i].vy);
    printf("P[%d].z = %5.10f \n",i,P[i].z);
    printf("P[%d].vz = %5.10f \n",i,P[i].vz);
    printf("P[%d].pmass = %g \n",i,P[i].pmass); //Diff by factor of h from Gadget units
    printf("P[%d].ptag = %d \n",i,P[i].ptag); */
   }

    if(fp == NULL) {
	printf("Something has gone horribly wrong! \n");
        return 0;
    }

  fclose(fp);

  header1.npart[0]= header1.npartTotal[0]= 0.;
  header1.npart[1]= header1.npartTotal[1]= 256*256*256;
  header1.npart[2]= header1.npartTotal[2]= 0.;
  header1.npart[3]= header1.npartTotal[3]= 0.;
  header1.npart[4]= header1.npartTotal[4]= 0.;
  header1.npart[5]= header1.npartTotal[5]= 0.;

  header1.mass[0]=  0.;
  header1.mass[1]=  P[0].pmass*h/1.e10; //Now the units are in h^-1 10^10 Msun
  header1.mass[2]=  0.;                  
  header1.mass[3]=  0.;
  header1.mass[4]=  0.;
  header1.mass[5]=  0.;

  header1.redshift   = 50.; // i.e. z_ini = 50
  header1.time   = 1./(1+header1.redshift); // i.e. a_start = 1/(1+z)
  header1.flag_sfr   = 0;
  header1.flag_feedback   = 0;
  header1.flag_cooling   = 0;
  header1.num_files   = 1; 

  int counter=0;
  float boxsize;
  for(i=0;i<5000;i++) {
    if ( P[i].x > 91. )  { 
      counter++ ; 
       } 
  }

  if (counter > 0) 
    { 
      header1.BoxSize = 256000;
    } else {
      header1.BoxSize = 64000;
    }

  //  header1.BoxSize   = boxsize*1000.; // units now in h^-1 kpc
  header1.Omega0   = 0.314;
  header1.OmegaLambda   = 0.686;
  header1.HubbleParam   = 0.71;

  // The following is meant to mimic the write commands in io.c from gadget-1.1

  int blklen;
  float dummy[3];

//  FILE *fd = fopen("ic_lg_gadget.dat", "w");

  #define BLKLEN fwrite(&blklen, sizeof(blklen), 1, stdout);

  blklen=sizeof(header1);
  BLKLEN;
  fwrite(&header1, sizeof(header1), 1, stdout);
  BLKLEN;

  blklen=NumPart*3*sizeof(float);

  BLKLEN;
  for(i=0;i<NumPart;i++)
    {

      dummy[0] = 1000*h*P[i].x; // Now the units are in h^-1 kpc
      dummy[1] = 1000*h*P[i].y;
      dummy[2] = 1000*h*P[i].z;
      fwrite(dummy,sizeof(float),3,stdout);
    }
  BLKLEN;

  BLKLEN;
  sqrta = (float)sqrt(header1.time);
  for(i=0;i<NumPart;i++)
    {
      dummy[0] = sqrta*P[i].vx;
      dummy[1] = sqrta*P[i].vy;
      dummy[2] = sqrta*P[i].vz;
      fwrite(dummy,sizeof(float),3,stdout);
    }
  BLKLEN;

  blklen=NumPart*sizeof(int);
  BLKLEN;
  for(i=0;i<NumPart;i++)
    {
      P[i].ptag += 1;
      fwrite(&P[i].ptag,sizeof(int),1,stdout);
    }
  BLKLEN;

//  fclose(fd);   
}