#include "ldstab.h"

//*********************************************************************************************************************
long int find_nearest_larger_nele(scan_table_single_def *scan_table_all_poi,long int nele,long int *nmax,long int ntables){

  bool shorter_instead_of_larger=false;
  return find_nearest_larger_nele_flip(scan_table_all_poi,nele,nmax,ntables,shorter_instead_of_larger);
}

//*********************************************************************************************************************
long int find_nearest_shorter_nele(scan_table_single_def *scan_table_all_poi,long int nele,long int *nmax,long int ntables){

  bool shorter_instead_of_larger=true;
  return find_nearest_larger_nele_flip(scan_table_all_poi,nele,nmax,ntables,shorter_instead_of_larger);
}

//*********************************************************************************************************************
long int find_nearest_larger_nele_flip(scan_table_single_def *scan_table_all_poi,long int nele,long int *nmax,long int ntables,bool shorter_instead_of_larger){
  
  //ricerca binaria

  //printf("finding nele %li\n",nele);


  long int itable_max=0;
  long int itable_min=ntables-1;
  long int itable_this=itable_max;    //ALL'INIZIO MI METTO QUI

  scan_table_single_def *table_min=scan_table_all_poi+itable_min;
  scan_table_single_def *table_max=scan_table_all_poi+itable_max;
  scan_table_single_def *table_this=scan_table_all_poi+itable_this;

  long int nele_min=table_min->nele;
  long int nele_max=table_max->nele;
  long int nele_this=table_this->nele;

  //controllo se ho nele è tra i valori tabulati
  if(nele_min > nele){
    return -1;
  }
  if(nele_max < nele){
    return -1;
  }

  //controllo se nele è agli estremi calcolati
  if(nele_min == nele){
    *nmax=nele_min;
    return itable_min;
  }
  if(nele_max == nele){
    *nmax=nele_max;
    return itable_max;
  }


  long int delta_table,semidelta_table;
  
 search_max_start:
  table_min=scan_table_all_poi+itable_min;
  table_max=scan_table_all_poi+itable_max;
  nele_min=table_min->nele;
  nele_max=table_max->nele;
  //printf("nele min max during scan:%li %li (%li) %li\n",nele_min,nele_max,nele,nele_this);

  //itable_min ha l'indice più alto
  //itable_max ha l'indice più basso
  delta_table=(itable_min-itable_max);
  semidelta_table=floor((double)(itable_min-itable_max+1)/2);

  //se il delta indici è 1, allora il valore cercato è itable_min oppure itable_max 
  if(delta_table <= 1){

    if(shorter_instead_of_larger == false){

      if(nele_min == nele){
	*nmax=nele_min;
	//printf("found 2.1: %li\n",nele_min);

	return itable_min;
      }
      *nmax=nele_max;       //quì c'è la scelta del nearest larger
      //printf("found 2.2: %li\n",nele_max);
      return itable_max;

    }else{

      if(nele_max == nele){
	*nmax=nele_max;
	//printf("found 3.1: %li %li itab: %li %li\n",nele_max,delta_table,itable_min,itable_max);
	return itable_max;
      }
      *nmax=nele_min;       //quì c'è la scelta del nearest **shorter**
      //printf("found 3.2: %li %li itab: %li %li\n",nele_min,delta_table,itable_min,itable_max);

      return itable_min;

    }

  }  //fine di if delta_table <= 1



  
  //  itable_this=itable_this+delta_table;
  itable_this=itable_max+semidelta_table;
  //printf("new table index: %li\n",itable_this);
  table_this=scan_table_all_poi+itable_this;
  nele_this=table_this->nele;
  if(nele == nele_this){    
    *nmax=nele_this;
    //printf("found 1: %li\n",nele_this);
    return itable_this;
  }
  
  if(nele_this > nele){    
    itable_max=itable_this;
  }else{
    itable_min=itable_this;

  }

  //printf("new tables limits: %li %li\n",itable_min,itable_max);
  goto search_max_start;
  
}







//*********************************************************************************************************************
long int find_nearest_larger_mele(scan_table_single_def *scan_table_this_poi,double mele,long int *mmax){

  bool shorter_instead_of_larger=false;
  return find_nearest_larger_mele_flip(scan_table_this_poi,mele,mmax,shorter_instead_of_larger);
}

//*********************************************************************************************************************
long int find_nearest_shorter_mele(scan_table_single_def *scan_table_this_poi,double mele,long int *mmax){

  bool shorter_instead_of_larger=true;
  return find_nearest_larger_mele_flip(scan_table_this_poi,mele,mmax,shorter_instead_of_larger);
}

//*********************************************************************************************************************
long int find_nearest_larger_mele_flip(scan_table_single_def *scan_table_this_poi,double mele,long int *mmax,bool shorter_instead_of_larger){
  
  //ricerca binaria

  //printf("searching for true mele %lf\n",mele);

  long int nmele=scan_table_this_poi->nmele;

  long int imele_min=0;
  long int imele_max=nmele-1;
  long int imele_this=imele_min;    //ALL'INIZIO MI METTO QUI

  double mele_min=*(scan_table_this_poi->mele_v+imele_min);
  double mele_max=*(scan_table_this_poi->mele_v+imele_max);


  //controllo se mele è tra i valori tabulati
  if(mele_min > mele){
    return -1;
  }
  if(mele_max < mele){
    return -1;
  }

  //controllo se mele è agli estremi calcolati
  if(mele_min == mele){
    *mmax=mele_min;
    return imele_min;
  }
  if(mele_max == mele){
    *mmax=mele_max;
    return imele_max;
  }


  long int delta_imele;
  long int semidelta_imele;
  double mele_this;
  
 search_mmax_start:
  mele_min=*(scan_table_this_poi->mele_v+imele_min);
  mele_max=*(scan_table_this_poi->mele_v+imele_max);



  //imele_max ha l'indice più alto
  //imele_min ha l'indice più basso
  delta_imele=imele_max-imele_min;
  semidelta_imele=floor((double)(imele_max-imele_min+1)/2);

  //se il delta indici è 1, allora il valore cercato è imele_min oppure imele_max 
  if(delta_imele <= 1){

    if(shorter_instead_of_larger == false){

      if(mele_min == mele){
	*mmax=mele_min;
	return imele_min;
      }
      *mmax=mele_max;       //quì c'è la scelta del nearest larger
      return imele_max;

    }else{

      if(mele_max == mele){
	*mmax=mele_max;
	return imele_max;
      }
      *mmax=mele_min;       //quì c'è la scelta del nearest **shorter**
      return imele_min;

    }

  }  //fine di if delta_imele <= 1



  
  imele_this=imele_min+semidelta_imele;
  mele_this=*(scan_table_this_poi->mele_v+imele_this);
  if(mele == mele_this){    
    *mmax=mele_this;
    return imele_this;
  }
  
  if(mele_this > mele){    
    imele_max=imele_this;
  }else{
    imele_min=imele_this;
  }

  goto search_mmax_start;
  
}



//#######################################################
scan_table_single_def *load_all_scan_table(long int *ntables){


  char datadir[STRINGLENMAX];
  strcpy(datadir,"./");
  scan_table_single_def *scan_table_all_poi=load_all_scan_table_fromdir(ntables,datadir);

  return scan_table_all_poi;
}

//#######################################################
scan_table_single_def *load_all_scan_table_fromdir(long int *ntables,char *datadir){

  char scan_stat_filename[STRINGLENMAX];
  FILE *scan_stat_file_poi;
  strcpy(scan_stat_filename,datadir);
  strcat(scan_stat_filename,"stab.txt");
  scan_stat_file_poi=fopen(scan_stat_filename,"r"); 

  char dummys[STRINGLENMAX];
  char release_date[STRINGLENMAX];  
  long int version,format_version;
  fscanf(scan_stat_file_poi,"%s %s %s %s %li\n",dummys,dummys,dummys,dummys,&format_version);
  printf("scan tables format version %li\n",format_version);
  if(format_version != 1){
    fprintf(stderr,"this version of \"ldstab\" read scan table file with format version 1 only\n");
    exit(-1);
  }
  fscanf(scan_stat_file_poi,"%s %s %s %li\n",dummys,dummys,dummys,&version);
  printf("scan tables version %li\n",version);
  fscanf(scan_stat_file_poi,"%s %s %s\n",dummys,dummys,release_date);
  printf("release date: %s\n",release_date);
  long int ncomments;
  fscanf(scan_stat_file_poi,"%s %s %li\n",dummys,dummys,&ncomments);
  long int icomment;
  for(icomment=0;icomment<ncomments;icomment++){
    fgets(dummys,STRINGLENMAX,scan_stat_file_poi);
    printf("COMMENT: %s",dummys);
  }

  fscanf(scan_stat_file_poi,"%8s %6li\n",dummys,ntables);
  //fgets(dummys,100,scan_stat_file_poi);
  printf("found %li tables\n",*ntables);
  
  scan_table_single_def *scan_table_all_poi=(scan_table_single_def *)malloc(*ntables*sizeof(scan_table_single_def));
  if(scan_table_all_poi == NULL){
    fprintf(stderr,"asked %li B\n",*ntables*sizeof(scan_table_single_def));
    fprintf(stderr,"unable to get memory for scan data\n");
    exit(-1);
  }

  long int nele;
  long int mele;
  double mdist_2p0sigma;
  double mdist_2p5sigma;
  double mdist_3p0sigma;
  double mdist_3p5sigma;
  double mdist_4p0sigma;
  double mdist_4p5sigma;
  double mdist_5p0sigma;
  long int imele;
  long int nmele;
  long int inele;

  scan_table_single_def *scan_table_single_poi;
  for(inele=0;inele<*ntables;inele++){
    
    fscanf(scan_stat_file_poi,"%5s %10li\n",dummys,&nele);    
    fscanf(scan_stat_file_poi,"%16c\n",dummys);
    fscanf(scan_stat_file_poi,"%26c\n",dummys);
    fscanf(scan_stat_file_poi,"%6s %6li\n",dummys,&nmele);
    //printf("table %li contains %li meles\n",nele,nmele);

    //scan_table_single_poi=scan_table_all_poi+inele;
    scan_table_single_poi=scan_table_all_poi+*ntables-1-inele;
    scan_table_single_poi->nele=nele;      
    scan_table_single_poi->nmele=nmele;

    scan_table_single_poi->mele_v=(long int *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_0p5sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_1p5sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_2p5sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_3p5sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_4p5sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_1p0sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_2p0sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_3p0sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_4p0sigma_v=(double *)malloc(sizeof(long int)*nmele);
    scan_table_single_poi->mdist_5p0sigma_v=(double *)malloc(sizeof(long int)*nmele);


   for(imele=0;imele<nmele;imele++){
      fscanf(scan_stat_file_poi,"%li %lf %lf %lf %lf %lf %lf %lf\n",&mele,&mdist_2p0sigma,&mdist_2p5sigma,&mdist_3p0sigma,&mdist_3p5sigma,&mdist_4p0sigma,&mdist_4p5sigma,&mdist_5p0sigma);
      //printf("%li %lf %lf %lf %lf %lf %lf %lf\n",mele,mdist_2p0sigma,mdist_2p5sigma,mdist_3p0sigma,mdist_3p5sigma,mdist_4p0sigma,mdist_4p5sigma,mdist_5p0sigma);
      *(scan_table_single_poi->mele_v+imele)=mele;
      *(scan_table_single_poi->mdist_0p5sigma_v+imele)=1.*0;
      *(scan_table_single_poi->mdist_1p0sigma_v+imele)=1.*0;
      *(scan_table_single_poi->mdist_1p5sigma_v+imele)=1.*0;
      *(scan_table_single_poi->mdist_2p0sigma_v+imele)=mdist_2p0sigma;
      *(scan_table_single_poi->mdist_2p5sigma_v+imele)=mdist_2p5sigma;
      *(scan_table_single_poi->mdist_3p0sigma_v+imele)=mdist_3p0sigma;
      *(scan_table_single_poi->mdist_3p5sigma_v+imele)=mdist_3p5sigma;
      *(scan_table_single_poi->mdist_4p0sigma_v+imele)=mdist_4p0sigma;
      *(scan_table_single_poi->mdist_4p5sigma_v+imele)=mdist_4p5sigma;
      *(scan_table_single_poi->mdist_5p0sigma_v+imele)=mdist_5p0sigma;
    }
  }

  fclose(scan_stat_file_poi);

  return scan_table_all_poi;
}

//**************************************************+
int interpol_mdist_scan_stat(long int nele,long int mele,scan_table_single_def *scan_table_all_poi, double *mdist_sigma_v,long int *ntables){

  //RITORNA negativo in caso di problemi
  long int iele;
  
  if(mele < 2){
    for(iele=0;iele<10;iele++){ 
      *(mdist_sigma_v+iele)=-1;
    }
    return -1;
  }
  
  if(mele > nele){
    for(iele=0;iele<10;iele++){ 
      *(mdist_sigma_v+iele)=-1;
    }
    return -1;
  }



  long int nmax;

  long int itabhigh=find_nearest_larger_nele(scan_table_all_poi,nele,&nmax,*ntables);
  if(itabhigh < 0){
    for(iele=0;iele<10;iele++){ 
      *(mdist_sigma_v+iele)=-1;
    }
    return -1;
  }


  long int nmin;
  long int itablow=find_nearest_shorter_nele(scan_table_all_poi,nele,&nmin,*ntables);
  if(itablow < 0){
    return -1;
  }
  

  //printf("found itab: %li %li and nmin/max %li %li\n",itablow,itabhigh,nmin,nmax);

  double itablev[2]={itablow,itabhigh};
  

   
  //adesso devo cercare sia per ilow, sia per ihigh i valori di mele più vicini a
  //quelli da interpolare
  
  
  double mdist_1sigma_lowhigh_v[2];
  double mdist_2sigma_lowhigh_v[2];
  double mdist_3sigma_lowhigh_v[2];
  double mdist_4sigma_lowhigh_v[2];
  double mdist_5sigma_lowhigh_v[2];
  double mdist_0p5sigma_lowhigh_v[2];
  double mdist_1p5sigma_lowhigh_v[2];
  double mdist_2p5sigma_lowhigh_v[2];
  double mdist_3p5sigma_lowhigh_v[2];
  double mdist_4p5sigma_lowhigh_v[2];

  
  long int itab;
  long int nmele;
  double mele_ref;
  scan_table_single_def *scan_table_this_poi;
  long int mmin,mmax;
  double mele_pre,mele_post;
  
  double mdist_1sigma_pre;
  double mdist_2sigma_pre;
  double mdist_3sigma_pre;
  double mdist_4sigma_pre;
  double mdist_5sigma_pre;
  
  double mdist_0p5sigma_pre;
  double mdist_1p5sigma_pre;
  double mdist_2p5sigma_pre;
  double mdist_3p5sigma_pre;
  double mdist_4p5sigma_pre;
  
  double mdist_1sigma_post;
  double mdist_2sigma_post;
  double mdist_3sigma_post;
  double mdist_4sigma_post;
  double mdist_5sigma_post;
  
  double mdist_0p5sigma_post;
  double mdist_1p5sigma_post;
  double mdist_2p5sigma_post;
  double mdist_3p5sigma_post;
  double mdist_4p5sigma_post;
  
  long imelehigh,imelelow;

  for(iele=0;iele<2;iele++){
    itab=itablev[iele];
    scan_table_this_poi=scan_table_all_poi+itab;
    nmele=scan_table_this_poi->nmele;
    
    mele_ref=(double)(mele-2)/nele*scan_table_this_poi->nele+2;
    
    imelehigh=find_nearest_larger_mele(scan_table_this_poi,mele_ref,&mmax);
    //da rivedere il seguente
    if(imelehigh == -1){
      imelehigh=nmele-1;
    }			 
    
    imelelow=find_nearest_shorter_mele(scan_table_this_poi,mele_ref,&mmin);
    //da rivedere il seguente
    if(imelelow == -1){
      imelelow=0;
    } 
    //printf("found imele: %li %li, meles: %li %li\n",imelelow,imelehigh,mmin,mmax);
    

    mele_pre=*(scan_table_this_poi->mele_v+imelelow);
    mele_post=*(scan_table_this_poi->mele_v+imelehigh);
    
    
    
    mdist_1sigma_pre=*(scan_table_this_poi->mdist_1p0sigma_v+imelelow);
    mdist_2sigma_pre=*(scan_table_this_poi->mdist_2p0sigma_v+imelelow);
    mdist_3sigma_pre=*(scan_table_this_poi->mdist_3p0sigma_v+imelelow);
    mdist_4sigma_pre=*(scan_table_this_poi->mdist_4p0sigma_v+imelelow);
    mdist_5sigma_pre=*(scan_table_this_poi->mdist_5p0sigma_v+imelelow);
    
    mdist_0p5sigma_pre=*(scan_table_this_poi->mdist_0p5sigma_v+imelelow);
    mdist_1p5sigma_pre=*(scan_table_this_poi->mdist_1p5sigma_v+imelelow);
    mdist_2p5sigma_pre=*(scan_table_this_poi->mdist_2p5sigma_v+imelelow);
    mdist_3p5sigma_pre=*(scan_table_this_poi->mdist_3p5sigma_v+imelelow);
    mdist_4p5sigma_pre=*(scan_table_this_poi->mdist_4p5sigma_v+imelelow);
    
    mdist_1sigma_post=*(scan_table_this_poi->mdist_1p0sigma_v+imelehigh);
    mdist_2sigma_post=*(scan_table_this_poi->mdist_2p0sigma_v+imelehigh);
    mdist_3sigma_post=*(scan_table_this_poi->mdist_3p0sigma_v+imelehigh);
    mdist_4sigma_post=*(scan_table_this_poi->mdist_4p0sigma_v+imelehigh);
    mdist_5sigma_post=*(scan_table_this_poi->mdist_5p0sigma_v+imelehigh);
    
    mdist_0p5sigma_post=*(scan_table_this_poi->mdist_0p5sigma_v+imelehigh);
    mdist_1p5sigma_post=*(scan_table_this_poi->mdist_1p5sigma_v+imelehigh);
    mdist_2p5sigma_post=*(scan_table_this_poi->mdist_2p5sigma_v+imelehigh);
    mdist_3p5sigma_post=*(scan_table_this_poi->mdist_3p5sigma_v+imelehigh);
    mdist_4p5sigma_post=*(scan_table_this_poi->mdist_4p5sigma_v+imelehigh);

    //printf("mele at 2 sigma pre/post: %le %le\n",mdist_2sigma_pre,mdist_2sigma_post);

    
    if(mele_post != mele_pre){
     mdist_1sigma_lowhigh_v[iele]=mdist_1sigma_pre+(mdist_1sigma_post-mdist_1sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_2sigma_lowhigh_v[iele]=mdist_2sigma_pre+(mdist_2sigma_post-mdist_2sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_3sigma_lowhigh_v[iele]=mdist_3sigma_pre+(mdist_3sigma_post-mdist_3sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_4sigma_lowhigh_v[iele]=mdist_4sigma_pre+(mdist_4sigma_post-mdist_4sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_5sigma_lowhigh_v[iele]=mdist_5sigma_pre+(mdist_5sigma_post-mdist_5sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     
     mdist_0p5sigma_lowhigh_v[iele]=mdist_0p5sigma_pre+(mdist_0p5sigma_post-mdist_0p5sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_1p5sigma_lowhigh_v[iele]=mdist_1p5sigma_pre+(mdist_1p5sigma_post-mdist_1p5sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_2p5sigma_lowhigh_v[iele]=mdist_2p5sigma_pre+(mdist_2p5sigma_post-mdist_2p5sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_3p5sigma_lowhigh_v[iele]=mdist_3p5sigma_pre+(mdist_3p5sigma_post-mdist_3p5sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
     mdist_4p5sigma_lowhigh_v[iele]=mdist_4p5sigma_pre+(mdist_4p5sigma_post-mdist_4p5sigma_pre)*(double)(mele_ref-mele_pre)/(mele_post-mele_pre);
   }else{
      mdist_1sigma_lowhigh_v[iele]=mdist_1sigma_pre;
      mdist_2sigma_lowhigh_v[iele]=mdist_2sigma_pre;
      mdist_3sigma_lowhigh_v[iele]=mdist_3sigma_pre;
      mdist_4sigma_lowhigh_v[iele]=mdist_4sigma_pre;
      mdist_5sigma_lowhigh_v[iele]=mdist_5sigma_pre;
      //printf("mdist pre 5sigma: %lf %lf %i\n",mdist_5sigma_pre,mdist_5sigma_post,iele);
      
      mdist_0p5sigma_lowhigh_v[iele]=mdist_0p5sigma_pre;
      mdist_1p5sigma_lowhigh_v[iele]=mdist_1p5sigma_pre;
      mdist_2p5sigma_lowhigh_v[iele]=mdist_2p5sigma_pre;
      mdist_3p5sigma_lowhigh_v[iele]=mdist_3p5sigma_pre;
      mdist_4p5sigma_lowhigh_v[iele]=mdist_4p5sigma_pre;
    }
    
  }// fine del ciclo for
  
  scan_table_single_def *scan_table_this_low_poi=scan_table_all_poi+itablow;
  scan_table_single_def *scan_table_this_high_poi=scan_table_all_poi+itabhigh;
  
  long int nele_pre=scan_table_this_low_poi->nele;
  long int nele_post=scan_table_this_high_poi->nele;
  mdist_1sigma_pre=mdist_1sigma_lowhigh_v[0];
  mdist_1sigma_post=mdist_1sigma_lowhigh_v[1];
  mdist_2sigma_pre=mdist_2sigma_lowhigh_v[0];
  mdist_2sigma_post=mdist_2sigma_lowhigh_v[1];
  mdist_3sigma_pre=mdist_3sigma_lowhigh_v[0];
  mdist_3sigma_post=mdist_3sigma_lowhigh_v[1];
  mdist_4sigma_pre=mdist_4sigma_lowhigh_v[0];
  mdist_4sigma_post=mdist_4sigma_lowhigh_v[1];
  mdist_5sigma_pre=mdist_5sigma_lowhigh_v[0];
  mdist_5sigma_post=mdist_5sigma_lowhigh_v[1];
  
  mdist_0p5sigma_pre=mdist_0p5sigma_lowhigh_v[0];
  mdist_0p5sigma_post=mdist_0p5sigma_lowhigh_v[1];
  mdist_1p5sigma_pre=mdist_1p5sigma_lowhigh_v[0];
  mdist_1p5sigma_post=mdist_1p5sigma_lowhigh_v[1];
  mdist_2p5sigma_pre=mdist_2p5sigma_lowhigh_v[0];
  mdist_2p5sigma_post=mdist_2p5sigma_lowhigh_v[1];
  mdist_3p5sigma_pre=mdist_3p5sigma_lowhigh_v[0];
  mdist_3p5sigma_post=mdist_3p5sigma_lowhigh_v[1];
  mdist_4p5sigma_pre=mdist_4p5sigma_lowhigh_v[0];
  mdist_4p5sigma_post=mdist_4p5sigma_lowhigh_v[1];
  

  double mdist_1sigma;
  double mdist_2sigma;
  double mdist_3sigma;
  double mdist_4sigma;
  double mdist_5sigma;

  double mdist_0p5sigma;
  double mdist_1p5sigma;
  double mdist_2p5sigma;
  double mdist_3p5sigma;
  double mdist_4p5sigma;



  if(nele_post != nele_pre){
    mdist_1sigma=mdist_1sigma_pre+(mdist_1sigma_post-mdist_1sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_2sigma=mdist_2sigma_pre+(mdist_2sigma_post-mdist_2sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_3sigma=mdist_3sigma_pre+(mdist_3sigma_post-mdist_3sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_4sigma=mdist_4sigma_pre+(mdist_4sigma_post-mdist_4sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_5sigma=mdist_5sigma_pre+(mdist_5sigma_post-mdist_5sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    
    mdist_0p5sigma=mdist_0p5sigma_pre+(mdist_0p5sigma_post-mdist_0p5sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_1p5sigma=mdist_1p5sigma_pre+(mdist_1p5sigma_post-mdist_1p5sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_2p5sigma=mdist_2p5sigma_pre+(mdist_2p5sigma_post-mdist_2p5sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_3p5sigma=mdist_3p5sigma_pre+(mdist_3p5sigma_post-mdist_3p5sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
    mdist_4p5sigma=mdist_4p5sigma_pre+(mdist_4p5sigma_post-mdist_4p5sigma_pre)*(double)(nele-nele_pre)/(nele_post-nele_pre);
  }else{
    mdist_1sigma=mdist_1sigma_pre;
    mdist_2sigma=mdist_2sigma_pre;
    mdist_3sigma=mdist_3sigma_pre;
    mdist_4sigma=mdist_4sigma_pre;
    mdist_5sigma=mdist_5sigma_pre;
    
    mdist_0p5sigma=mdist_0p5sigma_pre;
    mdist_1p5sigma=mdist_1p5sigma_pre;
    mdist_2p5sigma=mdist_2p5sigma_pre;
    mdist_3p5sigma=mdist_3p5sigma_pre;
    mdist_4p5sigma=mdist_4p5sigma_pre;
    //printf("mdist 5s: %le\n", mdist_5sigma_post);
  }
  
  
  
  *(mdist_sigma_v+0)=mdist_1sigma;
  *(mdist_sigma_v+1)=mdist_2sigma;
  *(mdist_sigma_v+2)=mdist_3sigma;
  *(mdist_sigma_v+3)=mdist_4sigma;
  *(mdist_sigma_v+4)=mdist_5sigma;
  
  *(mdist_sigma_v+5)=mdist_0p5sigma;
  *(mdist_sigma_v+6)=mdist_1p5sigma;
  *(mdist_sigma_v+7)=mdist_2p5sigma;
  *(mdist_sigma_v+8)=mdist_3p5sigma;
  *(mdist_sigma_v+9)=mdist_4p5sigma;
  
//  for(int iele=0;iele<10;iele++){
//    printf("mdist: %le\n", *(mdist_sigma_v+iele));
//  }
  
  return 0;
}


//******************************************************
void load_corr_table(corr_table_def *corr_table){

  char datadir[STRINGLENMAX];
  strcpy(datadir,"./");
  load_corr_table_fromdir(corr_table,datadir);
  return;
}
//******************************************************
void load_corr_table_fromdir(corr_table_def *corr_table,char *datadir){

  char correction_file[STRINGLENMAX];
  strcpy(correction_file,datadir);
  strcat(correction_file,"/ctab.txt");
  
  double nsigma;
  
  //determina gli elementi effettivi di nelev
  //determina gli elementi effettivi di ntscan_max
  long int nele_appo=0;
  double tscan_appo=0.;
  double pcoinc_appo=0.;
  double errm_appo=0.;
  double errp_appo=0.;

  FILE *lunc=fopen(correction_file,"r"); 
  if(lunc == NULL){
    fprintf(stderr,"opening corr file\n");
    exit(-1);
  }

  char dummys[STRINGLENMAX];
  char release_date[STRINGLENMAX];  
  long int version,format_version;
  fscanf(lunc,"%s %s %s %s %li\n",dummys,dummys,dummys,dummys,&format_version);
  printf("corr tables format version %li\n",format_version);
  if(format_version != 1){
    fprintf(stderr,"this version of \"ldstab\" read corr table file with format version 1 only\n");
    exit(-1);
  }
  fscanf(lunc,"%s %s %s %li\n",dummys,dummys,dummys,&version);
  printf("corr tables version %li\n",version);
  fscanf(lunc,"%s %s %s\n",dummys,dummys,release_date);
  printf("release date: %s\n",release_date);
  long int ncomments;
  fscanf(lunc,"%s %s %li\n",dummys,dummys,&ncomments);
  long int icomment;
  for(icomment=0;icomment<ncomments;icomment++){
    fgets(dummys,STRINGLENMAX,lunc);
    printf("COMMENT: %s",dummys);
  }





  long int nele_appo_old=-1;
  double tscan_appo_old=-1.;

  long int iele=-1;
  long int itscan=-1;
  long int itscan_max=-1;
  
  while(fscanf(lunc,"%li %lf %lf %lf %lf\n",&nele_appo,&tscan_appo,&pcoinc_appo,&errm_appo,&errp_appo) == 5 ){
    //printf("%li %lf %lf %lf %lf\n",nele_appo,tscan_appo,pcoinc_appo,errm_appo,errp_appo);
    if(nele_appo != nele_appo_old){
      iele=iele+1;
      nele_appo_old=nele_appo;
      corr_table->nelev[iele]=nele_appo;
      itscan=0;
    }else{
      itscan=itscan+1;
    }
    if(itscan > itscan_max){
      itscan_max=itscan;
    }
    //printf("   iele:%li, itscan:%li\n",iele,itscan);
    corr_table->tscanv[itscan]=tscan_appo;
    corr_table->pcoincv[iele][itscan]=pcoinc_appo/100;
    nsigma=calculate_normal_quantile(pcoinc_appo/100);
    corr_table->nsigmav[iele][itscan]=nsigma;

  }
  fclose(lunc);
  corr_table->nnele_max=iele+1;
  corr_table->ntscan_max=itscan+1;

  return;
}

//***********************************
double calculate_normal_quantile(double p){

  long int nele=1000;
  long int il=-1;
  long int ih=-1;
  long int iele;
  for(iele=0;iele<nele;iele++){
    double nsigma=(double)iele/nele*7.5;
    double prob=1.-(erf(nsigma/sqrt(2))/2+0.5);
    //printf("nsigma: %lf, prob: %le, p:%le\n",nsigma,prob,p);
    if(prob <= p){
      il=iele;
      goto after_calc_normal_for;
    }
  }
  after_calc_normal_for:

  //exit (0);
  if(il == -1){
    return 1.;
  }

  if(il == 0){
    return 0.;
  }

  ih=il-1;
  double nsigma_l=(double)il/1000*7.5;
  double nsigma_h=(double)ih/1000*7.5;
  double p_l=1.-(erf(nsigma_l/sqrt(2))/2+0.5);
  double p_h=1.-(erf(nsigma_h/sqrt(2))/2+0.5);
  double nsigma=nsigma_l+(nsigma_h-nsigma_l)*(p-p_l)/(p_h-p_l);
  
  double prob_recalc=1.-(erf(nsigma/sqrt(2))/2+0.5);
  //printf("nsigma %lf, prob %lf, Precalc %lf\n",nsigma,p,prob_recalc);
  //printf(" nsigma_h %lf, nsigma_l %lf, ih %li, il %li\n",nsigma_h,nsigma_l,ih,il);

  //char aaa[100];
  //gets(aaa);
  return nsigma;
}

//*****************************************
double interpol_correction_table(long int nele_tot,double tscan_found,corr_table_def *corr_table){

  long int nnele_max=corr_table->nnele_max;
  long int ntscan_max=corr_table->ntscan_max;

  long int iele_up=-1;
  long int iele_h;
  long int iele_l;
  long int iele;
  for(iele=0;iele<nnele_max;iele++){
    if(corr_table->nelev[iele] >= nele_tot){
      iele_up=iele;
      goto after_interpol_corr_for;
    }    
  }
 after_interpol_corr_for:
  if(iele_up == -1){
    iele_h=nnele_max-1;
    iele_l=iele_h;
  }else if(iele_up== 0){
    printf("nele_tot outside range, min is %li\n",corr_table->nelev[0]);
    iele_h=0;
    iele_l=iele_h;
  }else{
    iele_h=iele_up;
    iele_l=iele_h-1;
  }

  long int itscan_up=-1;
  long int itscan_h;
  long int itscan_l;
  long int itscan;
  for(itscan=0;itscan<ntscan_max;itscan++){
      if(corr_table->tscanv[itscan] >= tscan_found){
	itscan_up=itscan;
	goto after_interpol_corr_for2;
      }
    }
 after_interpol_corr_for2:
  if(itscan_up == -1){
    itscan_h=ntscan_max-1;
    itscan_l=itscan_h;
  }else if(iele_up == 0){
    printf("tscan outside range, using min: %lf\n",corr_table->tscanv[0]);
    itscan_h=0;
    itscan_l=itscan_h;
  }else{
    itscan_h=itscan_up;
    itscan_l=itscan_h-1;
  }

  double tscan_l=corr_table->tscanv[itscan_l];
  double tscan_h=corr_table->tscanv[itscan_h];

  long int nele_l=corr_table->nelev[iele_l];
  long int nele_h=corr_table->nelev[iele_h];

  double nsigma_at_nele_low;
  double nsigma_l_at_nele_low=corr_table->nsigmav[iele_l][itscan_l];
  double nsigma_h_at_nele_low=corr_table->nsigmav[iele_l][itscan_h];
  if(tscan_h != tscan_l){
    nsigma_at_nele_low=nsigma_l_at_nele_low+(nsigma_h_at_nele_low-nsigma_l_at_nele_low)*(tscan_found-tscan_l)/(tscan_h-tscan_l);
  }else{
    nsigma_at_nele_low=nsigma_l_at_nele_low;
  }

  double nsigma_at_nele_high;
  double nsigma_l_at_nele_high=corr_table->nsigmav[iele_h][itscan_l];
  double nsigma_h_at_nele_high=corr_table->nsigmav[iele_h][itscan_h];
  if(tscan_h != tscan_l){
    nsigma_at_nele_high=nsigma_l_at_nele_high+(nsigma_h_at_nele_high-nsigma_l_at_nele_high)*(tscan_found-tscan_l)/(tscan_h-tscan_l);
  }else{
    nsigma_at_nele_high=nsigma_l_at_nele_high;
  }

  double nsigma;
  if(nele_h != nele_l){
    nsigma=nsigma_at_nele_low+(nsigma_at_nele_high-nsigma_at_nele_low)*(1.*nele_tot-nele_l)/(nele_h-nele_l);
  }else{
    nsigma=nsigma_at_nele_low;
  }

   double prob=1.-(erf(nsigma/sqrt(2))/2+0.50);

  return prob;
}