00001
00002
00004 #include <stdlib.h>
00005 #include <stdio.h>
00006 #include <iostream>
00007 #include <fstream>
00008 #include <ctype.h>
00009 #include <math.h>
00010 #include <string.h>
00011
00012 #include "../nr.h"
00013 #include "../def.h"
00014
00015 using namespace std;
00016
00017 void readData(GlobExp *ex,Glob *globs,long expNr)
00018 {
00019 int lineNr=0;
00020 long k,l,index;
00021 long nLines=0,maxLeng=0;
00022 long nobsmax=ex[expNr].nobs;
00023 char c,*line,*arg;
00024 double *x;
00025 double **y;
00026 double **sig;
00027
00028 ifstream in;
00029
00030 FILE *fp;
00031
00032 fp=fopen(ex[expNr].fileName,"r");
00033 if (fp==NULL) {
00034 cerr << "Error opening data file " << ex[expNr].fileName << "\n";
00035 exit(1);
00036 }
00037
00038
00039 while((c=fgetc(fp))!=EOF)
00040 {
00041 for(k=1;c!='\n';k++)
00042 {
00043 if(k>maxLeng)
00044 maxLeng=k;
00045 c=fgetc(fp);
00046 }
00047 nLines++;
00048 }
00049 #ifdef DEBUGREADDATA
00050 *dbg << "Total number of lines in " << ex[expNr].fileName << " : " << nLines << endl;
00051 *dbg << "Maximal length of line in " << ex[expNr].fileName << " : " << maxLeng << endl;
00052 #endif
00053 fclose(fp);
00054
00055 in.open(ex[expNr].fileName);
00056
00057
00058 line=new char[maxLeng+2];
00059 arg=new char[maxLeng+2];
00060 x=dvector(1,nLines);
00061 y=dmatrix(1,nLines,1,ex[expNr].nobs);
00062 sig=dmatrix(1,nLines,1,ex[expNr].nobs);
00063
00064
00065 while(!in.eof())
00066 {
00067 in.getline(line,maxLeng+2,'\n');
00068 if(isspace(line[0])==FALSE && line[0]!= '#' && line[0] != '\0' )
00069 {
00070 lineNr++;
00071 #ifdef DEBUGREADDATA
00072 *dbg << "Extracting data line no. : " << lineNr << endl;
00073 #endif
00074 l=0;
00075 index=0;
00076 for(k=0;k<=strlen(line);k++)
00077 {
00078 if(isspace(line[k])!=FALSE)
00079 {
00080 arg[l]='\0';
00081 l=0;
00082 if(index==0)
00083 x[lineNr]=atof(arg);
00084 else if(index <=2*ex[expNr].nobs && (index % 2)==1)
00085 y[lineNr][(index+1)/2]=atof(arg);
00086 else if(index <=2*ex[expNr].nobs && (index % 2)==0)
00087 sig[lineNr][index/2]=atof(arg);
00088 index++;
00089 }
00090 else
00091 {
00092 arg[l]=line[k];
00093 l++;
00094 }
00095 }
00096 arg[l]='\0';
00097 if(index==0)
00098 x[lineNr]=atof(arg);
00099 else if(index <=2*ex[expNr].nobs && (index % 2)==1)
00100 y[lineNr][(index+1)/2]=atof(arg);
00101 else if(index <=2*ex[expNr].nobs && (index % 2)==0)
00102 sig[lineNr][index/2]=atof(arg);
00103
00104
00105 if((index % 2)==1)
00106 index--;
00107 if(index==0)
00108 {
00109 cerr << "Data line " << lineNr << " contains no data." << endl;
00110 exit(1);
00111 }
00112 if(index/2 <= nobsmax)
00113 nobsmax=index/2;
00114 }
00115
00116 }
00117 in.close();
00118
00119 ex[expNr].nobs=nobsmax;
00120 ex[expNr].nMeasure=lineNr;
00121 if(ex[expNr].nMeasure==0)
00122 {
00123 cerr << "File " << ex[expNr].fileName << " contains no data." << endl;
00124 exit(1);
00125 }
00126
00127
00128 for(k=2;k<=ex[expNr].nMeasure;k++)
00129 {
00130 if(x[k-1] >= x[k])
00131 {
00132 cerr << "Data not in temporal ascending order at dataline " << k << "\n";
00133 exit(1);
00134 }
00135 }
00136
00137 if(ex[expNr].fitstart == -1e99)
00138 ex[expNr].fitstart=x[1];
00139 if(ex[expNr].fitend == 1e99)
00140 ex[expNr].fitend=x[ex[expNr].nMeasure];
00141
00142
00143 index=1;
00144 for(k=1;k<=ex[expNr].nMeasure;k++)
00145 {
00146 if((x[k] >= ex[expNr].fitstart) && (x[k] <= ex[expNr].fitend))
00147 {
00148 index++;
00149 }
00150 }
00151 index--;
00152
00153
00154 ex[expNr].nMeasure=index;
00155
00156
00157 ex[expNr].xMeasure=dvector(1,ex[expNr].nMeasure);
00158 ex[expNr].yMeasure=dmatrix(1,ex[expNr].nMeasure,1,ex[expNr].nobs);
00159 ex[expNr].sigma=dmatrix(1,ex[expNr].nMeasure,1,ex[expNr].nobs);
00160
00161
00162 index=1;
00163 for(k=1;k<=lineNr;k++)
00164 {
00165 if((x[k] >= ex[expNr].fitstart) && (x[k] <= ex[expNr].fitend))
00166 {
00167 ex[expNr].xMeasure[index]=x[k];
00168 for(l=1;l<=ex[expNr].nobs;l++)
00169 {
00170 ex[expNr].yMeasure[index][l]=y[k][l];
00171 ex[expNr].sigma[index][l]=sig[k][l];
00172 }
00173 index++;
00174 }
00175 }
00176 index--;
00177
00178
00179 ex[expNr].firstMeasure=1;
00180 ex[expNr].lastMeasure=ex[expNr].nMeasure;
00181
00182
00183 free_dvector(x,1,nLines);
00184 free_dmatrix(y,1,nLines,1,ex[expNr].nobs);
00185 free_dmatrix(sig,1,nLines,1,ex[expNr].nobs);
00186 }