/*
''LyMAS: Predicting Large-Scale Lyman-alpha; Forest Statistics from the  Dark Matter Density Field''
  S. Peirani, D. H. Weinberg, S. Colombi, J. Blaizot & Y. Dubois
  Special thanks to Thierry Sousbie
*/

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


int swapI(int val)
{
    int out;
    const char *i=(const char *)&val;
    char *o=(char *)&out;

    o[3]=i[0];
    o[2]=i[1];
    o[1]=i[2];
    o[0]=i[3];

    return out; 
}

float swapF(float val)
{
    float out;
    const char *i=(const char *)&val;
    char *o=(char *)&out;

    o[3]=i[0];
    o[2]=i[1];
    o[1]=i[2];
    o[0]=i[3];

    return out; 
}

double swapD(double val)
{
    double out;
    const char *i=(const char *)&val;
    char *o=(char *)&out;

    o[7]=i[0];
    o[6]=i[1];
    o[5]=i[2];
    o[4]=i[3];
    o[3]=i[4];
    o[2]=i[5];
    o[1]=i[6];
    o[0]=i[7];

    return out; 
}




int main(int argc, char **argv)
{
   char *data_file;
   FILE *fd;
   int dummy;

   int F_space, DM_space;
   float F_smoothing, DM_smoothing;
   int F_dim, DM_dim;
   float log_rho_min, log_rho_max;
   double *P; 
   float log_rho, bin;
   int i, K;
   int swap_endians=0;

   if ( (argc<2) || (argc>2) ) 
     { 
       printf("\n\tUsage : read_prob data_file  \n\n");
       exit(1);  
     }

   data_file=argv[1];

   if (!(fd=fopen(data_file,"r")))
     {
       printf("can't open file `%s`\n",data_file);
       exit(0);
     }

   printf(" \n reading `%s' ...\n",data_file); fflush(stdout);

   //.....
   fread(&dummy, sizeof(dummy), 1, fd); 
   if (dummy!=8)
     {
       printf("  seems we need to swap indians...\n");
       swap_endians=1;
     }
   fread(&F_space, sizeof(int), 1, fd);   
   fread(&DM_space, sizeof(int), 1, fd);  
   fread(&dummy, sizeof(dummy), 1, fd);  
   //.....
   fread(&dummy, sizeof(dummy), 1, fd);
   fread(&F_smoothing, sizeof(float), 1, fd);   
   fread(&DM_smoothing, sizeof(float), 1, fd);  
   fread(&dummy, sizeof(dummy), 1, fd);
   //..... 
   fread(&dummy, sizeof(dummy), 1, fd);
   fread(&F_dim, sizeof(int), 1, fd);  
   fread(&DM_dim, sizeof(int), 1, fd);  
   fread(&dummy, sizeof(dummy), 1, fd);
   //..... 
   fread(&dummy, sizeof(dummy), 1, fd);
   fread(&log_rho_min, sizeof(float), 1, fd);   
   fread(&log_rho_max, sizeof(float), 1, fd);  
   fread(&dummy, sizeof(dummy), 1, fd);
   //..... 
   if(!(P=malloc(F_dim*DM_dim*sizeof(double))))
    {
      fprintf(stderr,"failed to allocate memory.\n");
      exit(0);
    }
   fread(&dummy, sizeof(dummy), 1, fd);  
   fread(P, sizeof(*P), F_dim*DM_dim , fd);
   fread(&dummy, sizeof(dummy), 1, fd); 
   //..... 
 
   fclose(fd);
   printf(" reading done\n\n");

   if(swap_endians)
     {
       F_space=swapI(F_space);
       DM_space=swapI(DM_space);
       F_smoothing=swapF(F_smoothing);
       DM_smoothing=swapF(DM_smoothing);
       F_dim=swapI(F_dim);
       DM_dim=swapI(DM_dim);
       log_rho_min=swapF(log_rho_min);
       log_rho_max=swapF(log_rho_max);

       for(i=0;i< F_dim*DM_dim;i++)
	 P[i]=swapD(P[i]);
      }

   if(F_space==0) printf(" Flux : real-space\n");
   else if (F_space==1) printf(" Flux : redshift-space\n");

   if(DM_space==0) printf(" DM density field : real-space\n");
   else if (DM_space==1) printf(" DM density field : redshift-space\n\n");

   printf(" Smoothing scales:\n");
   printf(" Flux : %.3f Mpc/h\n",F_smoothing);
   printf(" DM density field : %.3f Mpc/h\n\n",DM_smoothing);

   printf(" F_dim=%d\n",F_dim);
   printf(" DM_dim=%d\n",DM_dim);
   printf(" min(log(1+d))=%f\n",log_rho_min);
   printf(" max(log(1+d))=%f\n\n",log_rho_max);
     
   printf(" Probability distribution at which log(1+d)? (min=%.4f  max=%.4f): ",log_rho_min, log_rho_max);
   do
     {
       K=0;
       scanf("%f",&log_rho);
       if (log_rho<log_rho_min || log_rho>log_rho_max)   
         {
           printf(" You must choose  log(1+d) between %.4f and %.4f: ",log_rho_min, log_rho_max);
           K=1;
	 }
     }
   while (K!=0);

   bin=(log_rho_max-log_rho_min)/DM_dim;
   K=(int)((log_rho-log_rho_min)/bin) +1;
   //printf(" K=%d\n\n",K);
    
   printf(" write file: 'prob_test.dat'\n");
   fd = fopen("prob_test.dat", "w");
   bin=1./F_dim;
   for (i=0;i<F_dim;i++)
     fprintf (fd,"%f\t%f\n", bin*(i+0.5), P[F_dim*(K-1) + i]);

   fclose(fd);

   free(P);

   exit(0);    

}