/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Library: share/read_table-lib.c
*
* %Identification
* Written by: EF
* Date: Aug 28, 2002
* Origin: ILL
* Release: McStas 1.6
* Version: $Revision: 1.29 $
*
* This file is to be imported by components that may read data from table files
* It handles some shared functions. Embedded within instrument in runtime mode.
* Variable names have prefix 'mc_rt_' for 'McStas Read Table' to avoid conflicts
*
* Usage: within SHARE
* %include "read_table-lib"
*
* $Id: read_table-lib.c,v 1.29 2005/10/14 11:38:28 farhi Exp $
*
* $Log: read_table-lib.c,v $
* Revision 1.29  2005/10/14 11:38:28  farhi
* Corrected missing #define
*
* Revision 1.28  2005/10/12 14:04:29  farhi
* Added function to parse header, Table_ParseHeader(header, "symbol1", ... , NULL)
* Useful for complex sample components, as well as mcformat/mcconvert stuff.
*
* Revision 1.27  2005/10/05 08:50:53  farhi
* Extended buffer for Table line read (fgetl) to 64 ko instead of 4 ko.
* It failed with large matrices (e.g. more than 100 columns)
*
* Revision 1.26  2005/09/16 14:19:03  farhi
* Correct bug #56: SEGV when opening a non existent file with read-table. Was reported on usage of PowderN
*
* Revision 1.25  2005/08/31 14:50:29  farhi
* Fixed bugs when handling vectors with swapped (i,j) indexes
*
* Revision 1.24  2005/07/25 14:55:08  farhi
* DOC update:
* checked all parameter [unit] + text to be OK
* set all versions to CVS Revision
*
* Revision 1.23  2005/07/20 13:08:43  farhi
* Changed Table_Init calling sequence (overrides Table_Alloc)
*
* Revision 1.22  2005/07/12 14:46:26  farhi
* Added Table_Alloc to create a user empty Table
* and Table_SetElement
*
* Revision 1.21  2005/07/08 13:15:43  farhi
* Mismatch in argument swap
*
* Revision 1.20  2005/07/07 14:16:57  farhi
* Min and Max are computed for both row and column vectors
*
* Revision 1.19  2005/07/06 09:50:45  farhi
* Added check for non existing data file in Table_*_Array functions
*
* Revision 1.18  2005/07/06 08:44:28  farhi
* Display headers in Table_Info
* Also works with Arrays
*
* Revision 1.17  2005/07/05 14:30:27  farhi
* misprint
*
* Revision 1.16  2005/07/05 14:25:42  farhi
* added file size in t_Table structure
*
* Revision 1.15  2005/07/05 12:06:40  farhi
* added new functions for table Array handling
* to be used in Isotropic_sqw and mcformat
*
* Revision 1.13  2005/01/20 14:16:43  farhi
* New functions to read separately all numerical bmocks in a text data file
* Will be used for Data conversion from PGPLOT/McStas (mcformat tool)
*
* Revision 1.12  2004/09/10 15:12:02  farhi
* Make these libs easier to externalize (lower dependencies) and add comment about how to make these independent for external linkage.
*
* Revision 1.11  2004/09/09 13:48:02  farhi
* Code clean-up
*
* Revision 1.10  2004/09/03 13:46:50  farhi
* Correct misprint in comment
*
* Revision 1.9  2003/05/20 15:12:33  farhi
* malloc size for read table binary now needs less memory
*
* Revision 1.8  2003/02/11 12:28:46  farhi
* Variouxs bug fixes after tests in the lib directory
* mcstas_r  : disable output with --no-out.. flag. Fix 1D McStas output
* read_table:corrected MC_SYS_DIR -> MCSTAS define
* monitor_nd-lib: fix Log(signal) log(coord)
* HOPG.trm: reduce 4000 points -> 400 which is enough and faster to resample
* Progress_bar: precent -> percent parameter
* CS: ----------------------------------------------------------------------
*
* Revision 1.8  2003/02/06 14:14:41  farhi
* Corrected MC_SYS_DIR into MCSTAS definition of default lib location
*
* Revision 1.2 2002/12/19 12:48:07 ef
* Added binary import. Fixed Rebin. Added Stat.
*
* Revision 1.1 2002/08/29 11:39:00 ef
* Initial revision extracted from lib/optics/Monochromators...
*******************************************************************************/

#ifndef READ_TABLE_LIB_H
#include "read_table-lib.h"
#endif

/*******************************************************************************
* long Read_Table(t_Table *Table, char *name, int block_number)
*   ACTION: read a single Table from a text file
*   input   Table: pointer to a t_Table structure
*           name:  file name from which table should be extracted
*           block_number: if the file does contain more than one
*                 data block, then indicates which one to get (from index 1)
*                 a 0 value means append/catenate all
*   return  initialized single Table t_Table structure containing data, header, ...
*           number of read elements (-1: error, 0:header only)
* The routine stores any line starting with '#', '%' and ';' into the header
* File is opened, read and closed
* Other lines are interpreted as numerical data, and stored.
* Data block should be a rectangular matrix or vector.
* Data block may be rebined with Table_Rebin (also sort in ascending order)
*******************************************************************************/
  long Table_Read(t_Table *mc_rt_Table, char *mc_rt_File, long mc_rt_block_number)
  { /* reads all or a single data block from 'file' and returns a Table structure  */
    return(Table_Read_Offset(mc_rt_Table, mc_rt_File, mc_rt_block_number, NULL, 0));
  } /* end Table_Read */

/*******************************************************************************
* long Table_Read_Offset(t_Table *Table, char *name, int block_number, long *offset
*                        long max_lines)
*   ACTION: read a single Table from a text file, starting at offset
*     Same as Table_Read(..) except:
*   input   offset:    pointer to an offset (*offset should be 0 at start)
*           max_lines: max number of data rows to read from file (0 means all)
*   return  initialized single Table t_Table structure containing data, header, ...
*           number of read elements (-1: error, 0:header only)
*           updated *offset position (where end of reading occured)
*******************************************************************************/
  long Table_Read_Offset(t_Table *mc_rt_Table, char *mc_rt_File,
                         long mc_rt_block_number, long *mc_rt_offset,
                         long mc_rt_max_lines)
  { /* reads all/a data block in 'file' and returns a Table structure  */
    FILE *mc_rt_hfile;
    long  mc_rt_nelements;
    long  mc_rt_begin;
    long  mc_rt_filesize=0;
    struct stat mc_rt_stfile;

    if (!mc_rt_Table) return(-1);
    Table_Init(mc_rt_Table, 0, 0);
    if (!mc_rt_File)  return(-1);
    if (strlen(mc_rt_File) == 0) return (-1);
    mc_rt_hfile = fopen(mc_rt_File, "r");
    if(!mc_rt_hfile)
    {
      char mc_rt_path[256];
      char mc_rt_dir[256];

      if (!strchr(mc_rt_File, MC_PATHSEP_C))
      {
        strcpy(mc_rt_dir, getenv("MCSTAS") ? getenv("MCSTAS") : MCSTAS);
        sprintf(mc_rt_path, "%s%c%s%c%s", mc_rt_dir, MC_PATHSEP_C, "data", MC_PATHSEP_C, mc_rt_File);
        mc_rt_hfile = fopen(mc_rt_path, "r");
      }
      if(!mc_rt_hfile)
      {
        fprintf(stderr, "Error: Could not open input file '%s' (Table_Read)\n", mc_rt_File);
        return (-1);
      }
    }
    if (mc_rt_offset && *mc_rt_offset) fseek(mc_rt_hfile, *mc_rt_offset, SEEK_SET);
    else { stat(mc_rt_File,&mc_rt_stfile); mc_rt_filesize = mc_rt_stfile.st_size; }
    mc_rt_begin     = ftell(mc_rt_hfile);
    mc_rt_nelements = Table_Read_Handle(mc_rt_Table, mc_rt_hfile, mc_rt_block_number, mc_rt_max_lines);
    mc_rt_Table->begin = mc_rt_begin;
    mc_rt_Table->end   = ftell(mc_rt_hfile);
    mc_rt_Table->filesize = (mc_rt_filesize>0 ? mc_rt_filesize : 0);
    strncpy(mc_rt_Table->filename, mc_rt_File, 128);
    if (mc_rt_offset) *mc_rt_offset=mc_rt_Table->end;
    fclose(mc_rt_hfile);
    return(mc_rt_nelements);

  } /* end Table_Read_Offset */

/*******************************************************************************
* long Table_Read_Offset_Binary(t_Table *Table, char *File, char *type,
*                               long *offset, long rows, long columns)
*   ACTION: read a single Table from a binary file, starting at offset
*     Same as Table_Read_Offset(..) except that it handles binary files.
*   input   type: may be "float"/NULL or "double"
*           offset: pointer to an mc_rt_offset (*offset should be 0 at start)
*           rows   : number of rows (0 means read all)
*           columns: number of columns
*   return  initialized single Table t_Table structure containing data, header, ...
*           number of read elements (-1: error, 0:header only)
*           updated *offset position (where end of reading occured)
*******************************************************************************/
  long Table_Read_Offset_Binary(t_Table *mc_rt_Table, char *mc_rt_File, char *mc_rt_type,
                                long *mc_rt_offset, long mc_rt_rows, long mc_rt_columns)
  { /* reads all/a data block in binary 'file' and returns a Table structure  */
    long    mc_rt_nelements, mc_rt_sizeofelement;
    long    mc_rt_filesize;
    FILE   *mc_rt_hfile;
    struct stat mc_rt_stfile;
    double *mc_rt_data;
    long    mc_rt_i;
    long    mc_rt_begin;

    if (!mc_rt_Table) return(-1);

    Table_Init(mc_rt_Table, 0, 0);
    if (!mc_rt_File)  return(-1);

    stat(mc_rt_File,&mc_rt_stfile);
    mc_rt_filesize = mc_rt_stfile.st_size;
    mc_rt_hfile = fopen(mc_rt_File, "r");
    if(!mc_rt_hfile)
    {
      char mc_rt_path[256];
      char mc_rt_dir[256];

      if (!strchr(mc_rt_File, MC_PATHSEP_C))
      {
        strcpy(mc_rt_dir, getenv("MCSTAS") ? getenv("MCSTAS") : MCSTAS);
        sprintf(mc_rt_path, "%s%c%s%c%s", mc_rt_dir, MC_PATHSEP_C, "data", MC_PATHSEP_C, mc_rt_File);
        mc_rt_hfile = fopen(mc_rt_path, "r");
      }
      if(!mc_rt_hfile)
      {
        fprintf(stderr, "Error: Could not open input file '%s' (Table_Read_Offset_Binary)\n", mc_rt_File);
        return (-1);
      }
    }
    if (mc_rt_type && !strcmp(mc_rt_type,"double")) mc_rt_sizeofelement = sizeof(double);
    else  mc_rt_sizeofelement = sizeof(float);
    if (mc_rt_offset && *mc_rt_offset) fseek(mc_rt_hfile, *mc_rt_offset, SEEK_SET);
    mc_rt_begin     = ftell(mc_rt_hfile);
    if (mc_rt_rows && mc_rt_filesize > mc_rt_sizeofelement*mc_rt_columns*mc_rt_rows)
      mc_rt_nelements = mc_rt_columns*mc_rt_rows;
    else mc_rt_nelements = (long)(mc_rt_filesize/mc_rt_sizeofelement);
    if (!mc_rt_nelements || mc_rt_filesize <= *mc_rt_offset) return(0);
    mc_rt_data    = (double*)malloc(mc_rt_nelements*mc_rt_sizeofelement);
    if (!mc_rt_data) {
      fprintf(stderr,"Error: allocating %d elements for %s file '%s'. Too big (Table_Read_Offset_Binary).\n", mc_rt_nelements, mc_rt_type, mc_rt_File);
      exit(-1);
    }
    mc_rt_nelements = fread(mc_rt_data, mc_rt_sizeofelement, mc_rt_nelements, mc_rt_hfile);

    if (!mc_rt_data || !mc_rt_nelements)
    {
      fprintf(stderr,"Error: reading %d elements from %s file '%s' (Table_Read_Offset_Binary)\n", mc_rt_nelements, mc_rt_type, mc_rt_File);
      exit(-1);
    }
    mc_rt_Table->begin   = mc_rt_begin;
    mc_rt_Table->end     = ftell(mc_rt_hfile);
    if (mc_rt_offset) *mc_rt_offset=mc_rt_Table->end;
    fclose(mc_rt_hfile);
    mc_rt_data = (double*)realloc(mc_rt_data, (double)mc_rt_nelements*mc_rt_sizeofelement);
    /* copy file data into Table */
    if (mc_rt_type && !strcmp(mc_rt_type,"double")) mc_rt_Table->data = mc_rt_data;
    else {
      float  *mc_rt_s;
      double *mc_rt_dataf;
      mc_rt_s     = (float*)mc_rt_data;
      mc_rt_dataf = (double*)malloc(sizeof(double)*mc_rt_nelements);
      for (mc_rt_i=0; mc_rt_i<mc_rt_nelements; mc_rt_i++)
        mc_rt_dataf[mc_rt_i]=mc_rt_s[mc_rt_i];
      free(mc_rt_data);
      mc_rt_Table->data = mc_rt_dataf;
    }
    strcpy(mc_rt_Table->filename, mc_rt_File);
    mc_rt_Table->rows    = mc_rt_nelements/mc_rt_columns;
    mc_rt_Table->columns = mc_rt_columns;
    mc_rt_Table->array_length = 1;
    mc_rt_Table->block_number = 1;
    mc_rt_Table->filesize=mc_rt_filesize;

    Table_Stat(mc_rt_Table);

    return(mc_rt_nelements);
  } /* end Table_Read_Offset_Binary */

/*******************************************************************************
* long Read_Table_Handle(t_Table *Table, FILE *fid, int block_number, long max_lines)
*   ACTION: read a single Table from a text file handle (private)
*   input   Table:pointer to a t_Table structure
*           fid:  pointer to FILE handle
*           block_number: if the file does contain more than one
*                 data block, then indicates which one to get (from index 1)
*                 a 0 value means append/catenate all
*           max_lines: if non 0, only reads that number of lines
*   return  initialized single Table t_Table structure containing data, header, ...
*           modified Table t_Table structure containing data, header, ...
*           number of read elements (-1: error, 0:header only)
* The routine stores any line starting with '#', '%' and ';' into the header
* Other lines are interpreted as numerical data, and stored.
* Data block should be a rectangular matrix or vector.
* Data block may be rebined with Table_Rebin (also sort in ascending order)
*******************************************************************************/
  long Table_Read_Handle(t_Table *mc_rt_Table, FILE *mc_rt_hfile,
                         long mc_rt_block_number, long mc_rt_max_lines)
  { /* reads all/a data block from 'file' handle and returns a Table structure  */
    double *mc_rt_Data;
    char *mc_rt_Header;
    long  mc_rt_malloc_size         = 1024;
    long  mc_rt_malloc_size_h       = 4096;
    long  mc_rt_Rows = 0,   mc_rt_Columns = 0;
    long  mc_rt_count_in_array      = 0;
    long  mc_rt_count_in_header     = 0;
    long  mc_rt_block_Current_index = 0;
    char  mc_rt_flag_In_array       = 0;
    char  mc_rt_flag_End_row_loop   = 0;

    if (!mc_rt_Table) return(-1);
    Table_Init(mc_rt_Table, 0, 0);

    if(!mc_rt_hfile) {
       fprintf(stderr, "Error: File handle is NULL (Table_Read_Handle).\n");
       return (-1);
    }
    mc_rt_Header = (char*)  malloc(mc_rt_malloc_size_h*sizeof(char));
    mc_rt_Data   = (double*)malloc(mc_rt_malloc_size  *sizeof(double));
    if ((mc_rt_Header == NULL) || (mc_rt_Data == NULL)) {
       fprintf(stderr, "Error: Could not allocate Table and Header (Table_Read_Handle).\n");
       return (-1);
    }
    mc_rt_Header[0] = '\0';

    do { /* while (!mc_rt_flag_End_row_loop) */
      char  mc_rt_line[64*1024];
      long  mc_rt_back_pos=0;   /* ftell start of line */

      mc_rt_back_pos = ftell(mc_rt_hfile);
      if (fgets(mc_rt_line, 64*1024, mc_rt_hfile) != NULL) { /* analyse line */
        int mc_rt_i=0;
        char  mc_rt_flag_Store_into_header=0;
        /* first skip blank and tabulation characters */
        while (mc_rt_line[mc_rt_i] == ' ' || mc_rt_line[mc_rt_i] == '\t') mc_rt_i++;
        /* handle comments: stored in header */
        if ((mc_rt_line[mc_rt_i] == '#') || (mc_rt_line[mc_rt_i] == '%')
        || (mc_rt_line[mc_rt_i] == ';') || (mc_rt_line[mc_rt_i] == '/'))
        { /* line is a comment */
          mc_rt_flag_Store_into_header=1;
        } else {
          double mc_rt_X;

          /* get the number of columns splitting line with strtok */
          if (sscanf(mc_rt_line,"%lg ",&mc_rt_X) == 1)
          { /* line begins at least with one num */
            char  *mc_rt_InputTokens, *mc_rt_lexeme;
            char   mc_rt_flag_End_Line= 0;
            long   mc_rt_block_Num_Columns     = 0;

            mc_rt_InputTokens            = mc_rt_line;

            do { /* while (!mc_rt_flag_End_Line) */
              mc_rt_lexeme      = (char *)strtok(mc_rt_InputTokens, " ,;\t\n\r");
              mc_rt_InputTokens = NULL;
              if ((mc_rt_lexeme != NULL) && (strlen(mc_rt_lexeme) != 0))
              { /* reading line: the token is not empty */
                if (sscanf(mc_rt_lexeme,"%lg ",&mc_rt_X) == 1)
                { /* reading line: the token is a number in the line */
                  if (!mc_rt_flag_In_array)
                  { /* reading num: not already in a block: starts a new data block */
                    mc_rt_block_Current_index++;
                    mc_rt_flag_In_array    = 1;
                    mc_rt_block_Num_Columns= 0;
                    if (mc_rt_block_number)
                    { /* initialise a new data block */
                      mc_rt_Rows = 0;
                      mc_rt_count_in_array = 0;
                    } /* else append */
                  }
                  /* reading num: all blocks or selected block */
                  if ( mc_rt_flag_In_array &&  ((mc_rt_block_number == 0) || (mc_rt_block_number == mc_rt_block_Current_index)) )
                  {
                    /* starting block: already the desired number of rows ? */
                    if (mc_rt_block_Num_Columns == 0
                      && mc_rt_max_lines && mc_rt_Rows >= mc_rt_max_lines) {
                      mc_rt_flag_End_Line      = 1;
                      mc_rt_flag_End_row_loop  = 1;
                      mc_rt_flag_In_array      = 0;
                      /* reposition to begining of line (ignore line) */
                      fseek(mc_rt_hfile, mc_rt_back_pos, SEEK_SET);
                    } else { /* store into data array */
                      if (mc_rt_count_in_array >= mc_rt_malloc_size)
                      { /* realloc data buffer if necessary */
                        mc_rt_malloc_size = mc_rt_count_in_array+1024;
                        mc_rt_Data     = (double*)realloc(mc_rt_Data, mc_rt_malloc_size*sizeof(double));
                        if (mc_rt_Data == NULL)
                        {
                          fprintf(stderr, "Error: Can not re-allocate memory %i (Table_Read_Handle).\n", mc_rt_malloc_size*sizeof(double));
                          return (-1);
                        }
                      }
                      if (mc_rt_block_Num_Columns == 0) mc_rt_Rows++;
                      mc_rt_Data[mc_rt_count_in_array] = mc_rt_X;
                      mc_rt_count_in_array++;
                      mc_rt_block_Num_Columns++;
                    }
                  } /* reading num: end if mc_rt_flag_In_array */
                  else
                  { /* reading num: passed selected block */
                    if (mc_rt_block_number > mc_rt_block_Current_index)
                    { /* we finished to extract block -> force end of file reading */
                      mc_rt_flag_End_Line      = 1;
                      mc_rt_flag_End_row_loop  = 1;
                      mc_rt_flag_In_array  = 0;
                    }
                    /* else (if read all blocks) continue */
                  }
                } /* end reading num: end if sscanf mc_rt_lexeme -> numerical */
                else
                { /* reading line: the token is not numerical in that line. end block */
                  mc_rt_flag_End_Line = 1;
                  mc_rt_flag_In_array = 0;
                }
              }
              else
              { /* no more tokens in mc_rt_line */
                mc_rt_flag_End_Line = 1;
                if (mc_rt_block_Num_Columns) mc_rt_Columns = mc_rt_block_Num_Columns;
              }
            } while (!mc_rt_flag_End_Line); /* end while mc_rt_flag_End_Line */
          }
          else
          { /* ascii line: does not begin with a number: ignore line */
            mc_rt_flag_In_array          = 0;
            mc_rt_flag_Store_into_header = 1;
          }
        } /* end: if not line comment else numerical */
        if (mc_rt_flag_Store_into_header) { /* add line into header */
          mc_rt_count_in_header += strlen(mc_rt_line);
          if (mc_rt_count_in_header+4096 > mc_rt_malloc_size_h)
          { /* if succeed and in array : add (and realloc if necessary) */
            mc_rt_malloc_size_h = mc_rt_count_in_header+4096;
            mc_rt_Header     = (char*)realloc(mc_rt_Header, mc_rt_malloc_size_h*sizeof(char));
          }
          strncat(mc_rt_Header, mc_rt_line, 4096);
          mc_rt_flag_In_array  = 0; /* will start a new data block */
          /* exit line and file if passed desired block */
          if (mc_rt_block_number && mc_rt_block_number == mc_rt_block_Current_index) {
            mc_rt_flag_End_row_loop  = 1;
          }
        }
      } /* end: if fgets */
      else mc_rt_flag_End_row_loop = 1; /* else fgets : end of file */
    } while (!mc_rt_flag_End_row_loop); /* end while mc_rt_flag_End_row_loop */

    mc_rt_Table->block_number = mc_rt_block_number;
    mc_rt_Table->array_length = 1;
    if (mc_rt_count_in_header) mc_rt_Header    = (char*)realloc(mc_rt_Header, mc_rt_count_in_header*sizeof(char));
    mc_rt_Table->header       = mc_rt_Header;
    if (mc_rt_count_in_array*mc_rt_Rows*mc_rt_Columns == 0)
    {
      mc_rt_Table->rows         = 0;
      mc_rt_Table->columns      = 0;
      free(mc_rt_Data);
      return (0);
    }
    if (mc_rt_Rows * mc_rt_Columns != mc_rt_count_in_array)
    {
      fprintf(stderr, "Warning: Read_Table :%s Data has %li values that should be %li x %li\n", (!mc_rt_block_number ? " catenated" : ""), mc_rt_count_in_array, mc_rt_Rows, mc_rt_Columns);
      mc_rt_Columns = mc_rt_count_in_array; mc_rt_Rows = 1;
    }
    mc_rt_Data     = (double*)realloc(mc_rt_Data, mc_rt_count_in_array*sizeof(double));
    mc_rt_Table->data         = mc_rt_Data;
    mc_rt_Table->rows         = mc_rt_Rows;
    mc_rt_Table->columns      = mc_rt_Columns;

    Table_Stat(mc_rt_Table);
    return (mc_rt_count_in_array);

  } /* end Table_Read_Handle */

/*******************************************************************************
* long Rebin_Table(t_Table *Table)
*   ACTION: rebin a single Table, sorting 1st column in ascending order
*   input   Table: single table containing data.
*                  The data block is reallocated in this process
*   return  updated Table with increasing, evenly spaced first column (index 0)
*           number of data elements (-1: error, 0:header only)
*******************************************************************************/
  long Table_Rebin(t_Table *mc_rt_Table)
  {
    double mc_rt_new_step=0;
    long   mc_rt_i;
    long   mc_rt_tmp;
    char   mc_rt_monotonic = 1;
    /* performs linear interpolation on X axis (0-th column) */

    if (!mc_rt_Table) return(-1);
    if (!mc_rt_Table->data
    || mc_rt_Table->rows*mc_rt_Table->columns == 0 || !mc_rt_Table->step_x)
      return(0);
    mc_rt_tmp   = mc_rt_Table->rows;
    mc_rt_new_step = mc_rt_Table->step_x;
    for (mc_rt_i=0; mc_rt_i < mc_rt_Table->rows - 1; mc_rt_i++)
    {
      double mc_rt_current_step;
      double mc_rt_X, mc_rt_diff;
      mc_rt_X            = Table_Index(*mc_rt_Table,mc_rt_i  ,0);
      mc_rt_diff         = Table_Index(*mc_rt_Table,mc_rt_i+1,0) - mc_rt_X;
      mc_rt_current_step = fabs(mc_rt_diff);
      if ((mc_rt_Table->max_x - mc_rt_Table->min_x)*mc_rt_diff < 0 && mc_rt_monotonic && mc_rt_Table->columns > 1)
      {
        char mc_rt_buffer[256];
        if (!mc_rt_Table->block_number) strcpy(mc_rt_buffer, "catenated");
        else sprintf(mc_rt_buffer, "block %i", mc_rt_Table->block_number);
        fprintf(stderr, "Warning: Rebin_Table :%s Data from file '%s' (%li x %li) is not monotonic (at row %li)\n", mc_rt_buffer, mc_rt_Table->filename,
          mc_rt_Table->rows, mc_rt_Table->columns, mc_rt_i);
        mc_rt_monotonic = 0;
      }
      if (mc_rt_current_step > 0 && mc_rt_current_step < mc_rt_new_step) mc_rt_new_step = mc_rt_current_step;
      else mc_rt_tmp--;
    } /* for */
    if (fabs(mc_rt_new_step/mc_rt_Table->step_x) >= 0.98)
      return (mc_rt_Table->rows*mc_rt_Table->columns);
    if (mc_rt_tmp > 0 && mc_rt_new_step > 0 && mc_rt_Table->columns > 1)  /* table was not already evenly sampled */
    {
      long mc_rt_Length_Table;
      double *mc_rt_New_Table;
      /* modify step if leads to too many points */
      if  (mc_rt_Table->rows > 2000)
        if (mc_rt_new_step < mc_rt_Table->step_x)
          mc_rt_new_step = mc_rt_Table->step_x;
      if (mc_rt_new_step*10 < mc_rt_Table->step_x)
        mc_rt_new_step = mc_rt_Table->step_x/10;
      mc_rt_Length_Table = ceil(fabs(mc_rt_Table->max_x - mc_rt_Table->min_x)/mc_rt_new_step);
      mc_rt_New_Table    = (double*)malloc(mc_rt_Length_Table*mc_rt_Table->columns*sizeof(double));

      for (mc_rt_i=0; mc_rt_i < mc_rt_Length_Table; mc_rt_i++)
      {
        long   mc_rt_j;
        long   mc_rt_old_i;
        double mc_rt_X;
        double mc_rt_X1, mc_rt_X2, mc_rt_Y1, mc_rt_Y2;
        char   mc_rt_test=0;
        mc_rt_X = mc_rt_Table->min_x + mc_rt_i*mc_rt_new_step;
        mc_rt_New_Table[mc_rt_i*mc_rt_Table->columns] = mc_rt_X;
        /* look for index surrounding X in the old table -> index old_i, old-1 */
        for (mc_rt_old_i=1; mc_rt_old_i < mc_rt_Table->rows-1; mc_rt_old_i++)
        {
          mc_rt_X2 = Table_Index(*mc_rt_Table,mc_rt_old_i  ,0);
          mc_rt_X1 = Table_Index(*mc_rt_Table,mc_rt_old_i-1,0);
          if (mc_rt_Table->min_x < mc_rt_Table->max_x)
            mc_rt_test = ((mc_rt_X1 <= mc_rt_X) && (mc_rt_X < mc_rt_X2));
          else
            mc_rt_test = ((mc_rt_X2 <= mc_rt_X) && (mc_rt_X < mc_rt_X1));
          if (mc_rt_test) break;
        }

        for (mc_rt_j=1; mc_rt_j < mc_rt_Table->columns; mc_rt_j++)
        {
          mc_rt_Y2 = Table_Index(*mc_rt_Table,mc_rt_old_i  ,mc_rt_j);
          mc_rt_Y1 = Table_Index(*mc_rt_Table,mc_rt_old_i-1,mc_rt_j);
          if (mc_rt_X2-mc_rt_X1)
          {
          /* linear interpolation */
            double mc_rt_slope = (mc_rt_Y2-mc_rt_Y1)/(mc_rt_X2-mc_rt_X1);
            mc_rt_New_Table[mc_rt_i*mc_rt_Table->columns+mc_rt_j] = mc_rt_Y1+mc_rt_slope*(mc_rt_X-mc_rt_X1);
          }
          else
            mc_rt_New_Table[mc_rt_i*mc_rt_Table->columns+mc_rt_j] = mc_rt_Y2;
        }

      } /* end for  i */
      mc_rt_Table->rows = mc_rt_Length_Table;
      mc_rt_Table->step_x = mc_rt_new_step;
      free(mc_rt_Table->data);
      mc_rt_Table->data = mc_rt_New_Table;
    } /* end if tmp */
    return (mc_rt_Table->rows*mc_rt_Table->columns);
  } /* end Rebin_Table */

/*******************************************************************************
* double Table_Index(t_Table Table, long i, long j)
*   ACTION: read an element [i,j] of a single Table
*   input   Table: table containing data
*           i : index of row      (0:mc_rt_Rows-1)
*           j : index of column   (0:Columns-1)
*   return  Value = data[i][j]
* Returns Value from the i-th row, j-th column of Table
* Tests are performed on indexes i,j to avoid errors
*******************************************************************************/
  double Table_Index(t_Table mc_rt_Table, long mc_rt_i, long mc_rt_j)
  {
    long mc_rt_AbsIndex;

    if (mc_rt_i < 0)        mc_rt_i = 0;
    if (mc_rt_j < 0)        mc_rt_j = 0;
    /* handle vectors specifically */
    if (mc_rt_Table.columns == 1 || mc_rt_Table.rows == 1) {
      mc_rt_AbsIndex = mc_rt_i+mc_rt_j;
    } else {
      if (mc_rt_i >= mc_rt_Table.rows)    mc_rt_i = mc_rt_Table.rows-1;
      if (mc_rt_j >= mc_rt_Table.columns) mc_rt_j = mc_rt_Table.columns-1;
      mc_rt_AbsIndex = mc_rt_i*(mc_rt_Table.columns)+mc_rt_j;
    }
    if (mc_rt_AbsIndex >= mc_rt_Table.rows*mc_rt_Table.columns)
      mc_rt_AbsIndex = mc_rt_Table.rows*mc_rt_Table.columns-1;
    else if (mc_rt_AbsIndex < 0) mc_rt_AbsIndex=0;
    if (mc_rt_Table.data != NULL)
      return(mc_rt_Table.data[mc_rt_AbsIndex]);
    else
      return(0);
  } /* end Table_Index */

/*******************************************************************************
* void Table_SetElement(t_Table *Table, long i, long j, double value)
*   ACTION: set an element [i,j] of a single Table
*   input   Table: table containing data
*           i : index of row      (0:mc_rt_Rows-1)
*           j : index of column   (0:Columns-1)
*           value = data[i][j]
* Returns 0 in case of error
* Tests are performed on indexes i,j to avoid errors
*******************************************************************************/
  int Table_SetElement(t_Table *mc_rt_Table, long mc_rt_i, long mc_rt_j,
                        double mc_rt_value)
  {
    long mc_rt_AbsIndex;

    if (mc_rt_i < 0)        mc_rt_i = 0;
    if (mc_rt_i >= mc_rt_Table->rows)    mc_rt_i = mc_rt_Table->rows-1;
    if (mc_rt_j < 0)        mc_rt_j = 0;
    if (mc_rt_j >= mc_rt_Table->columns) mc_rt_j = mc_rt_Table->columns-1;
    mc_rt_AbsIndex = mc_rt_i*(mc_rt_Table->columns)+mc_rt_j;
    if (mc_rt_Table->data != NULL)
      mc_rt_Table->data[mc_rt_AbsIndex] = mc_rt_value;
    else return(0);
    return(1);
  } /* end Table_SetElement */

/*******************************************************************************
* double Table_Value(t_Table Table, double X, long j)
*   ACTION: read column [j] of a single Table at row which 1st column is X
*   input   Table: table containing data. Must be monotonic and evenly sampled
*                  use table_Rebin for that
*           X : data value in the first column (index 0)
*           j : index of column from which is extracted the Value (0:Columns-1)
*   return  Value = data[index for X][j]
* Returns Value from the j-th column of Table corresponding to the
* X value for the 1st column (index 0)
* Tests are performed (within Table_Index) on indexes i,j to avoid errors
* NOTE: data should rather be monotonic, and evenly sampled.
*******************************************************************************/
  double Table_Value(t_Table mc_rt_Table, double X, long j)
  {
    long   mc_rt_Index;
    double mc_rt_Value;

    if (mc_rt_Table.step_x != 0)
      mc_rt_Index = floor((X - mc_rt_Table.min_x)/mc_rt_Table.step_x);
    else mc_rt_Index=0;
    mc_rt_Value = Table_Index(mc_rt_Table, mc_rt_Index, j);

    return(mc_rt_Value);
  } /* end Table_Value */

/*******************************************************************************
* void Table_Free(t_Table *Table)
*   ACTION: free a single Table
*   return: empty Table
*******************************************************************************/
  void Table_Free(t_Table *mc_rt_Table)
  {
    if (!mc_rt_Table) return;
    if (mc_rt_Table->data   != NULL) free(mc_rt_Table->data);
    if (mc_rt_Table->header != NULL) free(mc_rt_Table->header);
    mc_rt_Table->data   = NULL;
    mc_rt_Table->header = NULL;
  } /* end Table_Free */

/******************************************************************************
* void Table_Info(t_Table Table)
*    ACTION: print informations about a single Table
*******************************************************************************/
  long Table_Info(t_Table mc_rt_Table)
  {
    char mc_rt_buffer[256];
    long ret=0;

    if (!mc_rt_Table.block_number) strcpy(mc_rt_buffer, "catenated");
    else sprintf(mc_rt_buffer, "block %i", mc_rt_Table.block_number);
    printf("Table from file '%s' (%s)", mc_rt_Table.filename, mc_rt_buffer);
    if (mc_rt_Table.filesize>0) printf(" of size %li", mc_rt_Table.filesize);
    if ((mc_rt_Table.data   != NULL) && (mc_rt_Table.rows*mc_rt_Table.columns))
    {
      printf(" is %li x %li ", mc_rt_Table.rows, mc_rt_Table.columns);
      if (mc_rt_Table.rows*mc_rt_Table.columns > 1)
        printf("(x=%g:%g).\n", mc_rt_Table.min_x, mc_rt_Table.max_x);
      else printf("(x=%g).\n", mc_rt_Table.min_x);
      ret = mc_rt_Table.rows*mc_rt_Table.columns;
    }
    else printf(" is empty.\n");
    if (mc_rt_Table.header && strlen(mc_rt_Table.header)) {
      char *header;
      int  i;
      header = malloc(80);
      if (!header) return(ret);
      for (i=0; i<80; header[i++]=0);
      if (strlen(mc_rt_Table.header) > 75) {
        strncpy(header, mc_rt_Table.header, 75-5);
        strcat( header, " ...");
      } else strcpy(header, mc_rt_Table.header);
      for (i=0; i<strlen(header); i++)
        if (header[i] == '\n' || header[i] == '\r') header[i] = ';';
      printf("  '%s'\n", header);
      free(header);
    }
    return(ret);
  } /* end Table_Info */

/******************************************************************************
* long Table_Init(t_Table *Table, m, n)
*   ACTION: initialise a Table to empty m by n table
*   return: empty Table
******************************************************************************/
long Table_Init(t_Table *mc_rt_Table, long rows, long columns)
{
  double *mc_rt_data=NULL;
  long   i;

  if (!mc_rt_Table) return(0);

  mc_rt_Table->header  = NULL;
  mc_rt_Table->filename[0]= '\0';
  mc_rt_Table->filesize= 0;
  mc_rt_Table->min_x   = 0;
  mc_rt_Table->max_x   = 0;
  mc_rt_Table->step_x  = 0;
  mc_rt_Table->block_number = 0;
  mc_rt_Table->array_length = 0;
  mc_rt_Table->begin   = 0;
  mc_rt_Table->end     = 0;

  if (rows*columns >= 1) {
    mc_rt_data    = (double*)malloc(rows*columns*sizeof(double));
    if (mc_rt_data) for (i=0; i < rows*columns; mc_rt_data[i++]=0);
    else {
      fprintf(stderr,"Error: allocating %d double elements."
                     "Too big (Table_Init).\n", rows*columns);
      rows = columns = 0;
    }
  }
  mc_rt_Table->rows    = (rows >= 1 ? rows : 0);
  mc_rt_Table->columns = (columns >= 1 ? columns : 0);
  mc_rt_Table->data    = mc_rt_data;
  return(mc_rt_Table->rows*mc_rt_Table->columns);
} /* end Table_Init */


/******************************************************************************
* void Table_Stat(t_Table *Table)
*   ACTION: computes min/max/mean step of 1st column for a single table (private)
*   return: updated Table
*******************************************************************************/
  static void Table_Stat(t_Table *mc_rt_Table)
  {
    long   mc_rt_i;
    double mc_rt_max_x, mc_rt_min_x;
    double mc_rt_row=1;
    double mc_rt_n;

    if (!mc_rt_Table) return;
    if (!mc_rt_Table->rows || !mc_rt_Table->columns) return;
    if (mc_rt_Table->rows == 1) mc_rt_row=0;
    mc_rt_max_x = mc_rt_Table->data[0];
    mc_rt_min_x = mc_rt_Table->data[(mc_rt_Table->rows-1)*mc_rt_Table->columns];
    mc_rt_n     = (mc_rt_row ? mc_rt_Table->rows : mc_rt_Table->columns);

    for (mc_rt_i=0; mc_rt_i < mc_rt_n; mc_rt_i++)
    {
      double mc_rt_X;
      mc_rt_X = (mc_rt_row ? Table_Index(*mc_rt_Table,mc_rt_i  ,0)
                           : Table_Index(*mc_rt_Table,0, mc_rt_i));
      if (mc_rt_X < mc_rt_min_x) mc_rt_min_x = mc_rt_X;
      if (mc_rt_X > mc_rt_max_x) mc_rt_max_x = mc_rt_X;
    } /* for */
    mc_rt_Table->max_x = mc_rt_max_x;
    mc_rt_Table->min_x = mc_rt_min_x;
    mc_rt_Table->step_x = (mc_rt_Table->max_x - mc_rt_Table->min_x)/mc_rt_n;
  } /* end Table_Stat */

/******************************************************************************
* t_Table *Table_Read_Array(char *File, long *blocks)
*   ACTION: read as many data blocks as available, iteratively from file
*   return: initialized t_Table array, last element is an empty Table.
*           the number of extracted blocks in non NULL pointer *blocks
*******************************************************************************/
  t_Table *Table_Read_Array(char *mc_rt_File, long *mc_rt_blocks)
  {
    t_Table *mc_rt_Table_Array=NULL;
    long mc_rt_offset=0;
    long mc_rt_block_number=0;
    long mc_rt_allocated=256;
    long mc_rt_nelements=1;

    /* fisrt allocate an initial empty t_Table array */
    mc_rt_Table_Array = (t_Table *)malloc(mc_rt_allocated*sizeof(t_Table));
    if (!mc_rt_Table_Array) {
      fprintf(stderr, "Error: Can not allocate memory %i (Table_Read_Array).\n",
         mc_rt_allocated*sizeof(t_Table));
      *mc_rt_blocks = 0;
      return (NULL);
    }

    while (mc_rt_nelements > 0)
    {
      t_Table mc_rt_Table;

      /* access file at mc_rt_offset and get following block */
      mc_rt_nelements = Table_Read_Offset(&mc_rt_Table, mc_rt_File, 1,
      &mc_rt_offset,0);
      /* if ok, set t_Table block number else exit loop */
      mc_rt_block_number++;
      mc_rt_Table.block_number = mc_rt_block_number;
      /* if t_Table array is not long enough, expand and realocate */
      if (mc_rt_block_number >= mc_rt_allocated-1) {
        mc_rt_allocated += 256;
        mc_rt_Table_Array = (t_Table *)realloc(mc_rt_Table_Array,
           mc_rt_allocated*sizeof(t_Table));
        if (!mc_rt_Table_Array) {
          fprintf(stderr, "Error: Can not re-allocate memory %i (Table_Read_Array).\n",
              mc_rt_allocated*sizeof(t_Table));
          *mc_rt_blocks = 0;
          return (NULL);
        }
      }
      /* store it into t_Table array */
      mc_rt_Table_Array[mc_rt_block_number-1] = mc_rt_Table;
      /* continues until we find an empty block */
    }
    /* send back number of extracted blocks */
    if (mc_rt_blocks) *mc_rt_blocks = mc_rt_block_number-1;

    /* now store total number of elements in Table array */
    for (mc_rt_offset=0; mc_rt_offset < mc_rt_block_number;
      mc_rt_Table_Array[mc_rt_offset++].array_length = mc_rt_block_number);

    return(mc_rt_Table_Array);
  } /* end Table_Read_Array */
/*******************************************************************************
* void Table_Free_Array(t_Table *Table)
*   ACTION: free a Table array
*******************************************************************************/
  void Table_Free_Array(t_Table *mc_rt_Table)
  {
    long mc_rt_index=0;
    if (!mc_rt_Table) return;
    do {
        if (mc_rt_Table[mc_rt_index].data || mc_rt_Table[mc_rt_index].header)
          Table_Free(&mc_rt_Table[mc_rt_index]);
        else mc_rt_index=-1;
    } while (mc_rt_index>= 0);
    free(mc_rt_Table);
  } /* end Table_Free_Array */

/******************************************************************************
* void Table_Info_Array(t_Table *Table)
*    ACTION: print informations about a Table array
*    return: number of elements in the Table array
*******************************************************************************/
  long Table_Info_Array(t_Table *mc_rt_Table)
  {
    long mc_rt_index=0;

    if (!mc_rt_Table) return(-1);
    while (mc_rt_index < mc_rt_Table[mc_rt_index].array_length
       && (mc_rt_Table[mc_rt_index].data || mc_rt_Table[mc_rt_index].header)
       && (mc_rt_Table[mc_rt_index].rows*mc_rt_Table[mc_rt_index].columns) ) {
      Table_Info(mc_rt_Table[mc_rt_index]);
      mc_rt_index++;
    }
    printf("This Table array contains %i elements\n", mc_rt_index);
    return(mc_rt_index);
  } /* end Table_Info_Array */

/******************************************************************************
* char **Table_ParseHeaderchar *header, symbol1, symbol2, ..., NULL)
*    ACTION: search for char* symbols in header and return their value or NULL
*            Last argument MUST be NULL
*    return: array of char* with line following each symbol, or NULL if not found
*******************************************************************************/
#ifndef MyNL_ARGMAX
#define MyNL_ARGMAX 50
#endif

char **Table_ParseHeader(char *header, ...){
  va_list ap;
  char exit_flag=0;
  char counter=0;
  char **ret;

  if (!header) return(NULL);

  ret = (char**)calloc(MyNL_ARGMAX, sizeof(char*));
  if (!ret) {
    printf("Table_ParseHeader: %s: Cannot allocate %i values array for Parser (Table_ParseHeader).\n",
      MyNL_ARGMAX);
    return(NULL);
  }

  va_start(ap, header);
  while(!exit_flag && counter < MyNL_ARGMAX-1)
  {
    char *arg_char;
    char *pos;
    /* get variable argument value as a char */
    arg_char = va_arg(ap, char *);
    ret[counter] = NULL;
    if (!arg_char){
      exit_flag = 1; break;
    }
    /* search for the symbol in the header */
    pos = strstr(header, arg_char);
    if (pos) {
      char *eol_pos;
      eol_pos = strchr(pos+strlen(arg_char), '\n');
      if (!eol_pos)
        eol_pos = strchr(pos+strlen(arg_char), '\r');
      if (!eol_pos) eol_pos = pos+strlen(pos)-1;
      ret[counter] = (char*)malloc(eol_pos - pos);
      if (!ret[counter]) {
        printf("Table_ParseHeader: %s: Cannot allocate value[%i] array for Parser (Table_ParseHeader).\n",
          counter);
        exit_flag = 1; break;
      }
      strncpy(ret[counter], pos+strlen(arg_char), eol_pos - pos - strlen(arg_char));
    }
    counter++;
  }
  va_end(ap);
  return(ret);
}

/* end of read_table-lib.c */
