// ########################################################################
// ##### read_data.c
// ##### - Please defien the reading file name "fname_in".
// ##### - Please search "OUTPUT SAMPLE". You can find the readed data and their units.
// ########################################################################



// ####################################
// ##### Includes

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



// ####################################
// ##### Options

#define  N_DATABLOCKS  (24)


// ## Units in GADGET
#define  UNIT_GADGET_LEN   (1.e03)  //>! kpc    [pc]
#define  UNIT_GADGET_VEL   (1.e05)  //>! km/sec [cm/sec]
#define  UNIT_GADGET_MASS  (1.e10)  //>! 1E10   [M_sun]


// ## MKS system of units: metre(m), kilogram(kg), and seconds(s)
#define  CONST_PERSEC    (3.0856776e+16)  //>! [m]
#define  CONST_BOLTZMANN (1.3806503e-23)  //>! [m^2 kg s^{-2} K^{-1}]
#define  CONST_PROTON_M  (1.67262158e-27) //>! [kg]
#define  CONST_SUN_M     (1.9891e+30)     //>! [kg]



// ####################################
// ##### Macros

#define  FLOAT_IN  float

#define  OpenFileAndCheck( File, Fname, Type )  if(!( File = fopen(Fname, Type))){printf( "####ERROR(%s/%s:L%d)####\n", __FILE__, __func__, __LINE__ );printf( "## Failed to open file '%s'.\n", Fname);exit( __LINE__ );}



// ####################################
// ##### Grobal variables

struct io_header
{
  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[96];
} header;

int NumPart, Ngas, ntot_withmasses;

struct particle_data
{

  FLOAT_IN pos[3];
  FLOAT_IN vel[3];
  int      id;
  int      splv;
  FLOAT_IN mass;
  FLOAT_IN ener;
  FLOAT_IN dens;
  FLOAT_IN abu[15];
  FLOAT_IN hsml;

  FLOAT_IN m_weight;
} P;


enum ind_chem_species {
    i_elec   , i_HI    , i_HII  , i_HeI , i_HeII
  , i_HeIII  , i_H2I   , i_H2II , i_HM  , i_DI
  , i_DII    , i_HDI   , i_HDII , i_DM  , i_GAMMA
  , N_ABU
};  //>! N_ABU = 16



// ####################################
// ##### Prototype declarations

long long int nseek( int num );



// ####################################
// ##### Main-programs

int main ( int argc, char **argv )
{
  // ##### Local Variables
  FILE   *file_R[N_DATABLOCKS];
  char   fname_in[512];
  int    i, k, n;
  double atime, a3inv, conv_len, conv_vel, conv_mass, conv_temp, conv_dens_Msuncc, conv_dens_gcc;


  //<! ##[ Please define the input filename ]##
  sprintf( fname_in, "./Run-Ref_N1E12"  ); //<! Run-Ref,  Nmax=1e12[/cc], 500[pc] around the newborn protostar
  //sprintf( fname_in, "./Run-B_N1E12"    ); //<! Run-B,    Nmax=1e12[/cc], 500[pc] around the newborn protostar
  //sprintf( fname_in, "./Run-J001_N1E12" ); //<! Run-J001, Nmax=1e12[/cc], 500[pc] around the newborn protostar
  //sprintf( fname_in, "./Run-J300_N1E08" ); //<! Run-J300, Nmax=1e08[/cc], 500[pc] around the newborn protostar


  //>! ## Header Data
  OpenFileAndCheck( file_R[0], fname_in, "r" );
  fseek( file_R[0], nseek(0), SEEK_CUR );
  fread( &header, sizeof(header), 1, file_R[0] );
  fclose( file_R[0] );

  printf( "\n## ----- Header Data\n");
  for( k=0; k<6; k++ ){
    if( header.npartTotal[k] > 0 ){
      printf( "## Type %1d : N_particle = %d", k, header.npartTotal[k] );
      if( header.mass[k] > 0 ) printf( " M_particle = %e [Msun]", header.mass[k]*1.e10 );
      printf( "\n" );
    }
  }
  printf( "## -----\n\n");


  //>! ## Summary of data size
  for(k=0, NumPart=0, ntot_withmasses=0; k<6; k++){
    NumPart += header.npartTotal[k];
    if(header.mass[k] == 0 && header.npartTotal[k] > 0){
      ntot_withmasses += header.npartTotal[k];
    }
  }
  Ngas = header.npartTotal[0];


  // ## Unit Conversion
  atime = a3inv = 1.;
  conv_len         = UNIT_GADGET_LEN*atime/header.HubbleParam;  //>! [(c)kpc/h]    -> [pc]
  conv_vel         = sqrt( atime );                             //>! [(c)km/sec]   -> [km/sec]
  conv_mass        = UNIT_GADGET_MASS/header.HubbleParam;       //>! [1e10 Msun/h] -> [Msun]
  conv_temp        = pow(UNIT_GADGET_VEL, 2.) / (CONST_BOLTZMANN*1.e7) * (CONST_PROTON_M*1.e3);
  conv_dens_Msuncc = conv_mass*CONST_SUN_M/CONST_PROTON_M * pow(conv_len*CONST_PERSEC*1.e2, -3);  // [1e10 Msun/h (kpc/h)^{-3} (comoving)] -> [g cm^{-3}]
  conv_dens_gcc    = conv_mass*CONST_SUN_M*1.e3 * pow(conv_len*CONST_PERSEC*1.e2, -3);  // [1e10 Msun/h (kpc/h)^{-3} (comoving)] -> [g cm^{-3}]


  //>! ########################################
  //>! ##### Open Files & Move Cursol

  for( n=1; n<N_DATABLOCKS; n++ ){
    OpenFileAndCheck( file_R[n], fname_in, "r" );
    fseek( file_R[n], nseek(n), SEEK_CUR );
  }


  //>! ########################################
  //>! ##### Read Data

  for( k=0; k<6; k++ ){
    if( header.npartTotal[k] > 0 ){
      for( n=0; n<header.npart[k]; n++ ){

        // ## Reading Data
        fread(&P.pos[0],   sizeof(FLOAT_IN),  3, file_R[ 1]);
        fread(&P.vel[0],   sizeof(FLOAT_IN),  3, file_R[ 2]);
        fread(&P.id,       sizeof(int),       1, file_R[ 3]);
        fread(&P.splv,     sizeof(int),       1, file_R[ 4]);
        if( header.mass[k] == 0 ) fread(&P.mass, sizeof(FLOAT_IN), 1, file_R[ 5]);
        else                      P.mass = header.mass[k];
        if( k == 0 ){
          fread(&P.ener,   sizeof(FLOAT_IN),  1, file_R[ 6]);
          fread(&P.dens,   sizeof(FLOAT_IN),  1, file_R[ 7]);
          for( i=0; i<N_ABU; i++ )
            fread(&P.abu[i], sizeof(FLOAT_IN), 1, file_R[8+i]);
          fread(&P.hsml,   sizeof(FLOAT_IN),  1, file_R[23]);

          for( i=0,P.m_weight=0; i<N_ABU-1; i++ ) P.m_weight += P.abu[i];
          P.m_weight = 1. / P.m_weight; //<! molecular weight
        }



        // ##[ OUTPUT SAMPLE ]##
        // - particle type (k): 0(gas), 1(dark matter), 4(star)

        if( n%(header.npart[k]/3)==0 ){
          printf( "## Type=%d, %d-th particle (id=%d, mass=%1.2e[Msun])\n", k, n, P.id, P.mass*conv_mass );
          printf( "## pos  = (%9.2e,%9.2e,%9.2e) [pc]\n"
                  , P.pos[0]*conv_len, P.pos[1]*conv_len, P.pos[2]*conv_len ); //<! Coordinate
          printf( "## vel  = (%9.2e,%9.2e,%9.2e) [km/s]\n"
                  , P.vel[0]*conv_vel, P.vel[1]*conv_vel, P.vel[2]*conv_vel ); //<! Velocity
          if( k==0 ){
            //<! for Gas Data
            printf( "## temp = %1.2e [K]\n", conv_temp*P.ener*(P.abu[i_GAMMA]-1.)*P.m_weight ); //<! Gas temperature
            printf( "## dens = %1.2e [gram/cm^3]\n", P.dens*conv_dens_gcc ); //<! Gas density
            printf( "## elec = %1.2e\n", P.abu[i_elec] ); //<! Chemical abundance
            printf( "## HI   = %1.2e\n", P.abu[i_HI] );
            printf( "## HII  = %1.2e\n", P.abu[i_HII] );
            printf( "## HeI  = %1.2e\n", P.abu[i_HeI] );
            printf( "## HeII = %1.2e\n", P.abu[i_HeII] );
            printf( "## HeIII= %1.2e\n", P.abu[i_HeIII] );
            printf( "## H2I  = %1.2e\n", P.abu[i_H2I] );
            printf( "## H2II = %1.2e\n", P.abu[i_H2II] );
            printf( "## HM   = %1.2e\n", P.abu[i_HM] );
            printf( "## DI   = %1.2e\n", P.abu[i_DI] );
            printf( "## DII  = %1.2e\n", P.abu[i_DII] );
            printf( "## HDI  = %1.2e\n", P.abu[i_HDI] );
            printf( "## HDII = %1.2e\n", P.abu[i_HDII] );
            printf( "## DM   = %1.2e\n", P.abu[i_DM] );
            printf( "## gamma= %1.2e\n", P.abu[i_GAMMA] ); //<! adiabatic index
            printf( "## hsml = %1.2e [pc]\n", P.hsml*conv_len ); //<! SPH smoothing length
          }
          printf( "\n" );
        }

        // ##[ OUTPUT SAMPLE ]##



      }  //>! Loop for Particle[Type]
    }
  }  //>! Loop for Type


  //>! ########################################
  //>! ##### Close Files

  for( n=1; n<N_DATABLOCKS; n++ ) fclose( file_R[n] );


  return 0;
}


long long int nseek( int num )
{
  long long int ldummy, ldummy_gap, ldummy_h, ldummy_f3, ldummy_i1, ldummy_fm, ldummy_fg;

  ldummy_gap = sizeof(int) * 2;  //>! Gap between data block
  ldummy_h   = sizeof(header);
  ldummy_f3  = sizeof(FLOAT_IN) * NumPart * 3;
  ldummy_i1  = sizeof(int)      * NumPart;
  ldummy_fm  = sizeof(FLOAT_IN) * ntot_withmasses;
  ldummy_fg  = sizeof(FLOAT_IN) * Ngas;

               ldummy  = sizeof(int);             //>! [ 0] Header
  if(num >  0) ldummy += ldummy_gap + ldummy_h;   //>! [ 1] Position
  if(num >  1) ldummy += ldummy_gap + ldummy_f3;  //>! [ 2] Velocity
  if(num >  2) ldummy += ldummy_gap + ldummy_f3;  //>! [ 3] ID
  if(num >  3) ldummy += ldummy_gap + ldummy_i1;  //>! [ 4] SpLv
  if(num >  4) ldummy += ldummy_gap + ldummy_i1;  //>! [ 5] Mass
  if(num >  5)
    if( ntot_withmasses > 0 ) ldummy += ldummy_gap + ldummy_fm;  //>! [ 6] InternalEnergy
  if(num >  6) ldummy += (ldummy_gap + ldummy_fg)*(num-6);  //>! [ 7] Density - [23] Hsml

  return ldummy;
}

#undef FLOAT_IN