/******************************************************************************* * * 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_idata = 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; iheader = 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 */