external_packages/hydro_from_external_file/arsenal.cpp

0001 // Ver 1.7.0
0002 // Zhi Qiu
0003 /*==========================================================================================
0004 Change logs: see arsenal.h
0005 ==========================================================================================*/
0006
0007 #include <vector>
0008 #include <iostream>
0009 #include <string>
0010 #include <sstream>
0011 #include <stdlib.h>
0012 #include <cmath>
0013 #include <iomanip>
0014 #include <cstdarg>
0015
0016 #include "arsenal.h"
0017
0018 #define OUTPUT_PRECISION 10
0019
0020 // for log_gamma function
0021 #define D1      -0.5772156649015328605195174
0022 #define D2      0.4227843350984671393993777
0023 #define D4      1.791759469228055000094023
0024 #define SQRTPI  0.9189385332046727417803297
0025 #define FRTBIG  1.42E+09
0026 #define PNT68   0.6796875
0027 #define XBIG    4.08E+36
0028 #define MACHINE_EPSILON 2.22044604925031e-016
0029
0030 using namespace std;
0031
0032 //**********************************************************************
0033 double sixPoint2dInterp(double x, double y,
0034     double v00, double v01, double v02, double v10, double v11, double v20)
0035 {
0036   /* Assume a 2d function f(x,y) has a quadratic form:
0037     f(x,y) = axx*x^2 + axy*x*y + ayy*y^2 + bx*x + by*y + c
0038     Also assume that:
0039     f(0,0)=v00, f(0,1)=v01, f(0,2)=v02, f(1,0)=v10, f(1,1)=v11, f(2,0)=v20
0040     Then its value at (x,y) can be solved.
0041     The best result is produced if x and y are between 0 and 1.
0042   */
0043
0044   // Calculate coefficients:
0045   double axx = 1.0/2.0*(v00 - 2*v10 + v20);
0046   double axy = v00 - v01 - v10 + v11;
0047   double ayy = 1.0/2.0*(v00 - 2*v01 + v02);
0048   double bx = 1.0/2.0*(-3.0*v00 + 4*v10 - v20);
0049   double by = 1.0/2.0*(-3.0*v00 + 4*v01 - v02);
0050   double c = v00;
0051
0052   // Calcualte f(x,y):
0053   return axx*x*x + axy*x*y + ayy*y*y + bx*x + by*y + c;
0054 }
0055
0056
0057 //**********************************************************************
0058 double interpCubicDirect(vector<double>* x, vector<double>* y, double x0)
0059 // Returns the interpreted value of y=y(x) at x=x0 using cubic polynomial interpolation method.
0060 // -- x,y: the independent and dependent double x0ables; x is assumed to be equal spaced and increasing
0061 // -- x0: where the interpolation should be performed
0062 {
0063   long size = x->size();
0064   if (size==1) {cout<<"interpCubicDirect warning: table size = 1"; return (*y)[0];}
0065   double dx = (*x)[1]-(*x)[0]; // increment in x
0066
0067   // if close to left end:
0068   if (abs(x0-(*x)[0])<dx*1e-30) return (*y)[0];
0069
0070   // find x's integer index
0071   long idx = floor((x0-(*x)[0])/dx);
0072
0073   if (idx<0 || idx>=size-1)
0074   {
0075     cout    << "interpCubicDirect: x0 out of bounds." << endl
0076             << "x ranges from " << (*x)[0] << " to " << (*x)[size-1] << ", "
0077             << "x0=" << x0 << ", " << "dx=" << dx << ", " << "idx=" << idx << endl;
0078     exit(-1);
0079   }
0080
0081   if (idx==0)
0082   {
0083     // use quadratic interpolation at left end
0084     double A0 = (*y)[0], A1 = (*y)[1], A2 = (*y)[2], deltaX = x0 - (*x)[0]; // deltaX is the increment of x0 compared to the closest lattice point
0085     return (A0-2.0*A1+A2)/(2.0*dx*dx)*deltaX*deltaX - (3.0*A0-4.0*A1+A2)/(2.0*dx)*deltaX + A0;
0086   }
0087   else if (idx==size-2)
0088   {
0089     // use quadratic interpolation at right end
0090     double A0 = (*y)[size-3], A1 = (*y)[size-2], A2 = (*y)[size-1], deltaX = x0 - ((*x)[0] + (idx-1)*dx);
0091     return (A0-2.0*A1+A2)/(2.0*dx*dx)*deltaX*deltaX - (3.0*A0-4.0*A1+A2)/(2.0*dx)*deltaX + A0;
0092   }
0093   else
0094   {
0095     // use cubic interpolation
0096     double A0 = (*y)[idx-1], A1 = (*y)[idx], A2 = (*y)[idx+1], A3 = (*y)[idx+2], deltaX = x0 - ((*x)[0] + idx*dx);
0097     //cout << A0 << "  " << A1 << "  " << A2 << "  " << A3 << endl;
0098     return (-A0+3.0*A1-3.0*A2+A3)/(6.0*dx*dx*dx)*deltaX*deltaX*deltaX
0099             + (A0-2.0*A1+A2)/(2.0*dx*dx)*deltaX*deltaX
0100             - (2.0*A0+3.0*A1-6.0*A2+A3)/(6.0*dx)*deltaX
0101             + A1;
0102   }
0103
0104 }
0105
0106
0107
0108
0109 //**********************************************************************
0110 double interpLinearDirect(vector<double>* x, vector<double>* y, double x0)
0111 // Returns the interpreted value of y=y(x) at x=x0 using linear interpolation method.
0112 // -- x,y: the independent and dependent double x0ables; x is assumed to be equal spaced and increasing
0113 // -- x0: where the interpolation should be performed
0114 {
0115   long size = x->size();
0116   if (size==1) {cout<<"interpLinearDirect warning: table size = 1"<<endl; return (*y)[0];}
0117   double dx = (*x)[1]-(*x)[0]; // increment in x
0118
0119   // if close to left end:
0120   if (abs(x0-(*x)[0])<dx*1e-30) return (*y)[0];
0121
0122   // find x's integer index
0123   long idx = floor((x0-(*x)[0])/dx);
0124
0125   if (idx<0 || idx>=size-1)
0126   {
0127     cout    << "interpLinearDirect: x0 out of bounds." << endl
0128             << "x ranges from " << (*x)[0] << " to " << (*x)[size-1] << ", "
0129             << "x0=" << x0 << ", " << "dx=" << dx << ", " << "idx=" << idx << endl;
0130     exit(-1);
0131   }
0132
0133   return (*y)[idx] + ((*y)[idx+1]-(*y)[idx])/dx*(x0-(*x)[idx]);
0134
0135 }
0136
0137
0138
0139
0140 //**********************************************************************
0141 double interpNearestDirect(vector<double>* x, vector<double>* y, double x0)
0142 // Returns the interpreted value of y=y(x) at x=x0 using nearest interpolation method.
0143 // -- x,y: the independent and dependent double x0ables; x is assumed to be equal spaced and increasing
0144 // -- x0: where the interpolation should be performed
0145 {
0146   long size = x->size();
0147   if (size==1) {cout<<"interpNearestDirect warning: table size = 1"<<endl; return (*y)[0];}
0148   double dx = (*x)[1]-(*x)[0]; // increment in x
0149
0150   // if close to left end:
0151   if (abs(x0-(*x)[0])<dx*1e-30) return (*y)[0];
0152
0153   // find x's integer index
0154   long idx = floor((x0-(*x)[0])/dx);
0155
0156   if (idx<0 || idx>=size-1)
0157   {
0158     cout    << "interpNearestDirect: x0 out of bounds." << endl
0159             << "x ranges from " << (*x)[0] << " to " << (*x)[size-1] << ", "
0160             << "x0=" << x0 << ", " << "dx=" << dx << ", " << "idx=" << idx << endl;
0161     exit(-1);
0162   }
0163
0164   return x0-(*x)[idx]>dx/2 ? (*y)[idx+1] : (*y)[idx];
0165
0166 }
0167
0168
0169
0170
0171 //**********************************************************************
0172 double interpCubicMono(vector<double>* x, vector<double>* y, double xx)
0173 // Note that this function does NOT perform well with small x and y table spacing; in which case use "direct" version instead.
0174 // Returns the interpreted value of y=y(x) at x=x0 using cubic polynomial interpolation method.
0175 // -- x,y: the independent and dependent double x0ables; x is *NOT* assumed to be equal spaced but it has to be increasing
0176 // -- xx: where the interpolation should be performed
0177 {
0178   long size = x->size();
0179   if (size==1) {cout<<"interpCubicMono warning: table size = 1"<<endl; return (*y)[0];}
0180
0181   // if close to left end:
0182   if (abs(xx-(*x)[0])<((*x)[1]-(*x)[0])*1e-30) return (*y)[0];
0183
0184   // find x's integer index
0185   long idx = binarySearch(x, xx);
0186
0187   if (idx<0 || idx>=size-1)
0188   {
0189     cout    << "interpCubicMono: x0 out of bounds." << endl
0190             << "x ranges from " << (*x)[0] << " to " << (*x)[size-1] << ", "
0191             << "xx=" << xx << ", " << "idx=" << idx << endl;
0192     exit(-1);
0193   }
0194
0195   if (idx==0)
0196   {
0197     // use linear interpolation at the left end
0198     return (*y)[0] + ( (*y)[1]-(*y)[0] )/( (*x)[1]-(*x)[0] )*( xx-(*x)[0] );
0199   }
0200   else if (idx==size-2)
0201   {
0202     // use linear interpolation at the right end
0203     return (*y)[size-2] + ( (*y)[size-1]-(*y)[size-2] )/( (*x)[size-1]-(*x)[size-2] )*( xx-(*x)[size-2] );
0204   }
0205   else
0206   {
0207     // use cubic interpolation
0208     long double y0 = (*y)[idx-1], y1 = (*y)[idx], y2 = (*y)[idx+1], y3 = (*y)[idx+2];
0209     long double y01=y0-y1, y02=y0-y2, y03=y0-y3, y12=y1-y2, y13=y1-y3, y23=y2-y3;
0210     long double x0 = (*x)[idx-1], x1 = (*x)[idx], x2 = (*x)[idx+1], x3 = (*x)[idx+2];
0211     long double x01=x0-x1, x02=x0-x2, x03=x0-x3, x12=x1-x2, x13=x1-x3, x23=x2-x3;
0212     long double x0s=x0*x0, x1s=x1*x1, x2s=x2*x2, x3s=x3*x3;
0213     long double denominator = x01*x02*x12*x03*x13*x23;
0214     long double C0, C1, C2, C3;
0215     C0 = (x0*x02*x2*x03*x23*x3*y1
0216           + x1*x1s*(x0*x03*x3*y2+x2s*(-x3*y0+x0*y3)+x2*(x3s*y0-x0s*y3))
0217           + x1*(x0s*x03*x3s*y2+x2*x2s*(-x3s*y0+x0s*y3)+x2s*(x3*x3s*y0-x0*x0s*y3))
0218           + x1s*(x0*x3*(-x0s+x3s)*y2+x2*x2s*(x3*y0-x0*y3)+x2*(-x3*x3s*y0+x0*x0s*y3))
0219           )/denominator;
0220     C1 = (x0s*x03*x3s*y12
0221           + x2*x2s*(x3s*y01+x0s*y13)
0222           + x1s*(x3*x3s*y02+x0*x0s*y23-x2*x2s*y03)
0223           + x2s*(-x3*x3s*y01-x0*x0s*y13)
0224           + x1*x1s*(-x3s*y02+x2s*y03-x0s*y23)
0225           )/denominator;
0226     C2 = (-x0*x3*(x0s-x3s)*y12
0227           + x2*(x3*x3s*y01+x0*x0s*y13)
0228           + x1*x1s*(x3*y02+x0*y23-x2*y03)
0229           + x2*x2s*(-x3*y01-x0*y13)
0230           + x1*(-x3*x3s*y02+x2*x2s*y03-x0*x0s*y23)
0231           )/denominator;
0232     C3 = (x0*x03*x3*y12
0233           + x2s*(x3*y01+x0*y13)
0234           + x1*(x3s*y02+x0s*y23-x2s*y03)
0235           + x2*(-x3s*y01-x0s*y13)
0236           + x1s*(-x3*y02+x2*y03-x0*y23)
0237           )/denominator;
0238 /*    cout  << x0s*x03*x3s*y12 << "  "
0239           <<  x2*x2s*(x3s*y01+x0s*y13) << "   "
0240           <<  x1s*(x3*x3s*y02+x0*x0s*y23-x2*x2s*y03) << "  "
0241           <<  x2s*(-x3*x3s*y01-x0*x0s*y13) << "  "
0242           <<  x1*x1s*(-x3s*y02+x2s*y03-x0s*y23) << endl;
0243     cout << denominator << endl;
0244
0245     cout << x0 << " " << x1 << "  " << x2 << "  " << x3 << endl;
0246     cout << y0 << " " << y1 << "  " << y2 << "  " << y3 << endl;
0247     cout << C0 << "  " << C1 << "  " << C2 << "  " << C3 << endl;*/
0248     return C0 + C1*xx + C2*xx*xx + C3*xx*xx*xx;
0249   }
0250
0251 }
0252
0253
0254
0255
0256
0257 //**********************************************************************
0258 double interpLinearMono(vector<double>* x, vector<double>* y, double xx)
0259 // Returns the interpreted value of y=y(x) at x=x0 using linear interpolation method.
0260 // -- x,y: the independent and dependent double x0ables; x is *NOT* assumed to be equal spaced but it has to be increasing
0261 // -- xx: where the interpolation should be performed
0262 {
0263   long size = x->size();
0264   if (size==1) {cout<<"interpLinearMono warning: table size = 1"<<endl; return (*y)[0];}
0265
0266   // if close to left end:
0267   if (abs(xx-(*x)[0])<((*x)[1]-(*x)[0])*1e-30) return (*y)[0];
0268
0269   // find x's integer index
0270   long idx = binarySearch(x, xx);
0271
0272   if (idx<0 || idx>=size-1)
0273   {
0274     cout    << "interpLinearMono: x0 out of bounds." << endl
0275             << "x ranges from " << (*x)[0] << " to " << (*x)[size-1] << ", "
0276             << "xx=" << xx << ", " << "idx=" << idx << endl;
0277     exit(-1);
0278   }
0279
0280   return (*y)[idx] + ( (*y)[idx+1]-(*y)[idx] )/( (*x)[idx+1]-(*x)[idx] )*( xx-(*x)[idx] );
0281
0282 }
0283
0284
0285
0286
0287 //**********************************************************************
0288 double interpNearestMono(vector<double>* x, vector<double>* y, double xx)
0289 // Returns the interpreted value of y=y(x) at x=x0 using nearest interpolation method.
0290 // -- x,y: the independent and dependent double x0ables; x is *NOT* assumed to be equal spaced but it has to be increasing
0291 // -- xx: where the interpolation should be performed
0292 {
0293   long size = x->size();
0294   if (size==1) {cout<<"interpNearestMono warning: table size = 1"<<endl; return (*y)[0];}
0295
0296   // if close to left end:
0297   if (abs(xx-(*x)[0])<((*x)[1]-(*x)[0])*1e-30) return (*y)[0];
0298
0299   // find x's integer index
0300   long idx = binarySearch(x, xx);
0301
0302   if (idx<0 || idx>=size-1)
0303   {
0304     cout    << "interpNearestMono: x0 out of bounds." << endl
0305             << "x ranges from " << (*x)[0] << " to " << (*x)[size-1] << ", "
0306             << "xx=" << xx << ", " << "idx=" << idx << endl;
0307     exit(-1);
0308   }
0309
0310   return xx-(*x)[idx] > (*x)[idx+1]-xx ? (*y)[idx+1] : (*y)[idx];
0311
0312 }
0313
0314
0315
0316
0317 //**********************************************************************
0318 double invertFunc(double (*func)(double), double y, double xL, double xR, double dx, double x0, double relative_accuracy)
0319 //Purpose:
0320 //  Return x=func^(-1)(y) using Newton method.
0321 //  -- func: double 1-argument function to be inverted
0322 //  -- xL: left boundary (for numeric derivative)
0323 //  -- xR: right boundary (for numeric derivative)
0324 //  -- dx: step (for numeric derivative)
0325 //  -- x0: initial value
0326 //  -- y: the value to be inverted
0327 //  -- Returns inverted value
0328 //Solve: f(x)=0 with f(x)=table(x)-y => f'(x)=table'(x)
0329 {
0330   double accuracy;
0331   int tolerance;
0332
0333   double XX1, XX2; // used in iterations
0334   double F0, F1, F2, F3, X1, X2; // intermedia variables
0335   int impatience; // number of iterations
0336
0337
0338   // initialize parameters
0339   accuracy = dx*relative_accuracy;
0340
0341   tolerance = 60;
0342   impatience = 0;
0343
0344   // initial value, left point and midxle point
0345   XX2 = x0;
0346   XX1 = XX2-10*accuracy; // this value 10*accuracy is meanless, just to make sure the check in the while statement goes through
0347
0348   while (abs(XX2-XX1)>accuracy)
0349   {
0350     XX1 = XX2; // copy values
0351
0352     // value of function at XX
0353     F0 = (*func)(XX1) - y; // the value of the function at this point
0354
0355     // decide X1 and X2 for differentiation
0356     if (XX1>xL+dx)
0357       X1 = XX1 - dx;
0358     else
0359       X1 = xL;
0360
0361     if (XX1<xR-dx)
0362       X2 = XX1 + dx;
0363     else
0364       X2 = xR;
0365
0366     // get values at X1 and X2
0367     F1 = (*func)(X1);
0368     F2 = (*func)(X2);
0369     F3 = (F1-F2)/(X1-X2); // derivative at XX1
0370
0371     XX2 = XX1 - F0/F3; // Newton's mysterious method
0372
0373     impatience = impatience + 1;
0374     //cout << "impatience=" << impatience << endl;
0375     if (impatience>tolerance)
0376     {
0377       cout << "invertFunc: " << "max number of iterations reached." << endl;
0378       exit(-1);
0379     }
0380
0381   } // <=> abs(XX2-XX1)>accuracy
0382
0383   return XX2;
0384 }
0385
0386
0387
0388 //**********************************************************************
0389 vector<double> *zq_x_global, *zq_y_global;
0390 double invertTableDirect_hook(double xx) {return interpCubicDirect(zq_x_global,zq_y_global,xx);}
0391 double invertTableDirect(vector<double>* x, vector<double>* y, double y0, double x0, double relative_accuracy)
0392 // Return x0=y^(-1)(y0) for y=y(x); use interpCubic and invertFunc.
0393 //  -- x,y: the independent and dependent variables. x is assumed to be equal-spaced.
0394 //  -- y0: where the inversion should be performed.
0395 //  -- x0: initial guess
0396 {
0397   long size = x->size();
0398   if (size==1) return (*y)[0];
0399   zq_x_global = x; zq_y_global = y;
0400   return invertFunc(&invertTableDirect_hook, y0, (*x)[0], (*x)[size-1], (*x)[1]-(*x)[0], x0, relative_accuracy);
0401 }
0402
0403
0404
0405
0406 //**********************************************************************
0407 vector<double> stringToDoubles(string str)
0408 // Return a vector of doubles from the string "str". "str" should
0409 // be a string containing a line of data.
0410 {
0411   stringstream sst(str+" "); // add a blank at the end so the last data will be read
0412   vector<double> valueList;
0413   double val;
0414   sst >> val;
0415   while (sst.eof()==false)
0416   {
0417     valueList.push_back(val);
0418     sst >> val;
0419   }
0420   return valueList;
0421 }
0422
0423
0424 //**********************************************************************
0425 double stringToDouble(string str)
0426 // Return the 1st doubles number read from the string "str". "str" should be a string containing a line of data.
0427 {
0428   stringstream sst(str+" "); // add a blank at the end so the last data will be read
0429   double val;
0430   sst >> val;
0431   return val;
0432 }
0433
0434
0435
0436 //**********************************************************************
0437 vector< vector<double>* >* readBlockData(istream &stream_in)
0438 // Return a nested vector of vector<double>* object. Each column of data
0439 // is stored in a vector<double> array and the collection is the returned
0440 // object. Data are read from the input stream "stream_in". Each line
0441 // of data is processed by the stringToDoubles function. Note that the
0442 // data block is dynamicall allocated and is not release within the
0443 // function.
0444 // Note that all "vectors" are "new" so don't forget to delete them.
0445 // Warning that also check if the last line is read correctly. Some files
0446 // are not endded properly and the last line is not read.
0447 {
0448   vector< vector<double>* >* data;
0449   vector<double> valuesInEachLine;
0450   long lineSize;
0451   long i; // temp variable
0452   char buffer[99999]; // each line should be shorter than this
0453
0454   // first line:
0455   stream_in.getline(buffer,99999);
0456   valuesInEachLine = stringToDoubles(buffer);
0457   // see if it is empty:
0458   lineSize = valuesInEachLine.size();
0459   if (lineSize==0)
0460   {
0461     // empty:
0462     cout << "readBlockData warning: input stream has empty first row; no data read" << endl;
0463     return NULL;
0464   }
0465   else
0466   {
0467     // not empty; allocate memory:
0468     data = new vector< vector<double>* >(lineSize);
0469     for (i=0; i<lineSize; i++) (*data)[i] = new vector<double>;
0470   }
0471
0472   // rest of the lines:
0473   while (stream_in.eof()==false)
0474   {
0475     // set values:
0476     for (i=0; i<lineSize; i++) (*(*data)[i]).push_back(valuesInEachLine[i]);
0477     // next line:
0478     stream_in.getline(buffer,99999);
0479     valuesInEachLine = stringToDoubles(buffer);
0480   }
0481
0482   return data;
0483 }
0484
0485
0486 //**********************************************************************
0487 void releaseBlockData(vector< vector<double>* >* data)
0488 // Use to delete the data block allocated by readBlockData function.
0489 {
0490   if (data)
0491   {
0492     for (unsigned long i=0; i<data->size(); i++) delete (*data)[i];
0493     delete data;
0494   }
0495 }
0496
0497
0498 //**********************************************************************
0499 // From Wikipedia --- the free encyclopeida
0500 //
0501 // Recursive auxiliary function for adaptiveSimpsons() function below
0502 //
0503 double adaptiveSimpsonsAux(double (*f)(double), double a, double b, double epsilon,
0504                            double S, double fa, double fb, double fc, int bottom) {
0505   double c = (a + b)/2, h = b - a;
0506   double d = (a + c)/2, e = (c + b)/2;
0507   double fd = f(d), fe = f(e);
0508   double Sleft = (h/12)*(fa + 4*fd + fc);
0509   double Sright = (h/12)*(fc + 4*fe + fb);
0510   double S2 = Sleft + Sright;
0511   if (bottom <= 0 || fabs(S2 - S) <= 15*epsilon)
0512     return S2 + (S2 - S)/15;
0513   return adaptiveSimpsonsAux(f, a, c, epsilon/2, Sleft,  fa, fc, fd, bottom-1) +
0514          adaptiveSimpsonsAux(f, c, b, epsilon/2, Sright, fc, fb, fe, bottom-1);
0515 }
0516 //
0517 // Adaptive Simpson's Rule
0518 //
0519 double adaptiveSimpsons(double (*f)(double),   // ptr to function
0520                         double a, double b,  // interval [a,b]
0521                         double epsilon,  // error tolerance
0522                         int maxRecursionDepth) {   // recursion cap
0523   double c = (a + b)/2, h = b - a;
0524   double fa = f(a), fb = f(b), fc = f(c);
0525   double S = (h/6)*(fa + 4*fc + fb);
0526   return adaptiveSimpsonsAux(f, a, b, epsilon, S, fa, fb, fc, maxRecursionDepth);
0527 }
0528
0529
0530 //**********************************************************************
0531 double qiu_simpsons(double (*f)(double), // ptr to function
0532                     double a, double b, // interval [a,b]
0533                     double epsilon, int maxRecursionDepth) // recursion maximum
0534 // My version of the adaptive simpsons integration method.
0535 {
0536   double f_1=f(a)+f(b), f_2=0., f_4=0.; // sum of values of f(x) that will be weighted by 1, 2, 4 respectively, depending on where x is
0537   double sum_previous=0., sum_current=0.; // previous and current sum (intgrated value)
0538
0539   long count = 1; // how many new mid-points are there
0540   double length = (b-a), // helf length of the interval
0541          step = length/count; // mid-points are located at a+(i+0.5)*step, i=0..count-1
0542
0543   int currentRecursionDepth = 1;
0544
0545   f_4 = f(a+0.5*step); // mid point of [a,b]
0546   sum_current = (length/6)*(f_1 + f_2*2. + f_4*4.); // get the current sum
0547
0548   do
0549   {
0550     sum_previous = sum_current; // record the old sum
0551     f_2 += f_4; // old mid-points with weight 4 will be new mid-points with weight 2
0552
0553     count*=2; // increase number of mid-points
0554     step/=2.0; // decrease jumping step by half
0555     f_4 = 0.; // prepare to sum up f_4
0556     for (int i=0; i<count; i++) f_4 += f(a+step*(i+0.5)); // sum up f_4
0557
0558     sum_current = (length/6/count)*(f_1 + f_2*2. + f_4*4.); // calculate current sum
0559     //cout << sum_current << endl;
0560
0561     if (currentRecursionDepth>maxRecursionDepth)
0562     {
0563       cout << endl << "Warning qiu_simpsons: maximum recursion depth reached!" << endl << endl;
0564       break; // safety treatment
0565     }
0566     else currentRecursionDepth++;
0567
0568   } while (abs(sum_current-sum_previous)>epsilon);
0569
0570   return sum_current;
0571 }
0572
0573 //**********************************************************************
0574 double qiu_simpsonsRel(double (*f)(double), // ptr to function
0575                     double a, double b, // interval [a,b]
0576                     double epsilon, int maxRecursionDepth) // recursion maximum
0577 // My version of the adaptive simpsons integration method.
0578 {
0579   double f_1=f(a)+f(b), f_2=0., f_4=0.; // sum of values of f(x) that will be weighted by 1, 2, 4 respectively, depending on where x is
0580   double sum_previous=0., sum_current=0.; // previous and current sum (intgrated value)
0581
0582   long count = 1; // how many new mid-points are there
0583   double length = (b-a), // helf length of the interval
0584          step = length/count; // mid-points are located at a+(i+0.5)*step, i=0..count-1
0585
0586   int currentRecursionDepth = 1;
0587
0588   f_4 = f(a+0.5*step); // mid point of [a,b]
0589   sum_current = (length/6)*(f_1 + f_2*2. + f_4*4.); // get the current sum
0590
0591   do
0592   {
0593     sum_previous = sum_current; // record the old sum
0594     f_2 += f_4; // old mid-points with weight 4 will be new mid-points with weight 2
0595
0596     count*=2; // increase number of mid-points
0597     step/=2.0; // decrease jumping step by half
0598     f_4 = 0.; // prepare to sum up f_4
0599     for (int i=0; i<count; i++) f_4 += f(a+step*(i+0.5)); // sum up f_4
0600
0601     sum_current = (length/6/count)*(f_1 + f_2*2. + f_4*4.); // calculate current sum
0602     //cout << sum_current << endl;
0603
0604     if (currentRecursionDepth>maxRecursionDepth)
0605     {
0606       cout << endl << "Warning qiu_simpsons: maximum recursion depth reached!" << endl << endl;
0607       break; // safety treatment
0608     }
0609     else currentRecursionDepth++;
0610
0611   } while (abs(sum_current-sum_previous)/(sum_current-sum_previous)>epsilon);
0612
0613   return sum_current;
0614 }
0615
0616
0617 //**********************************************************************
0618 string toLower(string str)
0619 // Convert all character in string to lower case
0620 {
0621   string tmp = str;
0622   for (string::iterator it=tmp.begin(); it<=tmp.end(); it++) *it = tolower(*it);
0623   return tmp;
0624 }
0625
0626 //**********************************************************************
0627 string trim(string str)
0628 // Convert all character in string to lower case
0629 {
0630   string tmp = str;
0631   long number_of_char = 0;
0632   for (size_t ii=0; ii<str.size(); ii++)
0633     if (str[ii]!=' ' && str[ii]!='\t')
0634     {
0635       tmp[number_of_char]=str[ii];
0636       number_of_char++;
0637     }
0638   tmp.resize(number_of_char);
0639   return tmp;
0640 }
0641
0642
0643 //**********************************************************************
0644 long binarySearch(vector<double>* A, double value, bool skip_out_of_range)
0645 // Return the index of the largest number less than value in the list A
0646 // using binary search. Index starts with 0.
0647 // If skip_out_of_range is set to true, then it will return -1 for those
0648 // samples that are out of the table range.
0649 {
0650    int length = A->size();
0651    int idx_i, idx_f, idx;
0652    idx_i = 0;
0653    idx_f = length-1;
0654    if(value > (*A)[idx_f])
0655    {
0656       if (skip_out_of_range) return -1;
0657       cout << "binarySearch: desired value is too large, exceeding the end of the table." << endl;
0658       exit(-1);
0659    }
0660    if(value < (*A)[idx_i])
0661    {
0662       if (skip_out_of_range) return -1;
0663       cout << "binarySearch: desired value is too small, exceeding the beginning of table." << endl;
0664       exit(-1);
0665    }
0666    idx = (int) floor((idx_f+idx_i)/2.);
0667    while((idx_f-idx_i) > 1)
0668    {
0669      if((*A)[idx] < value)
0670         idx_i = idx;
0671      else
0672         idx_f = idx;
0673      idx = (int) floor((idx_f+idx_i)/2.);
0674    }
0675    return(idx_i);
0676 }
0677
0678
0679 //**********************************************************************
0680 long double gamma_function(long double x)
0681 // gamma.cpp -- computation of gamma function.
0682 //      Algorithms and coefficient values from "Computation of Special
0683 //      Functions", Zhang and Jin, John Wiley and Sons, 1996.
0684 // Returns gamma function of argument 'x'.
0685 //
0686 // NOTE: Returns 1e308 if argument is a negative integer or 0,
0687 //      or if argument exceeds 171.
0688 {
0689   int i,k,m;
0690   long double ga,gr,r=0,z;
0691
0692   static long double g[] = {
0693     1.0,
0694     0.5772156649015329,
0695    -0.6558780715202538,
0696    -0.420026350340952e-1,
0697     0.1665386113822915,
0698    -0.421977345555443e-1,
0699    -0.9621971527877e-2,
0700     0.7218943246663e-2,
0701    -0.11651675918591e-2,
0702    -0.2152416741149e-3,
0703     0.1280502823882e-3,
0704    -0.201348547807e-4,
0705    -0.12504934821e-5,
0706     0.1133027232e-5,
0707    -0.2056338417e-6,
0708     0.6116095e-8,
0709     0.50020075e-8,
0710    -0.11812746e-8,
0711     0.1043427e-9,
0712     0.77823e-11,
0713    -0.36968e-11,
0714     0.51e-12,
0715    -0.206e-13,
0716    -0.54e-14,
0717     0.14e-14};
0718
0719   if (x > 171.0) return 1e308;    // This value is an overflow flag.
0720   if (x == (int)x) {
0721     if (x > 0.0) {
0722       ga = 1.0;               // use factorial
0723       for (i=2;i<x;i++) {
0724         ga *= i;
0725       }
0726      }
0727      else
0728       ga = 1e308;
0729    }
0730    else {
0731     if (fabs(x) > 1.0) {
0732       z = fabs(x);
0733       m = (int)z;
0734       r = 1.0;
0735       for (k=1;k<=m;k++) {
0736         r *= (z-k);
0737       }
0738       z -= m;
0739     }
0740     else
0741       z = x;
0742     gr = g[24];
0743     for (k=23;k>=0;k--) {
0744       gr = gr*z+g[k];
0745     }
0746     ga = 1.0/(gr*z);
0747     if (fabs(x) > 1.0) {
0748       ga *= r;
0749       if (x < 0.0) {
0750         ga = -M_PI/(x*ga*sin(M_PI*x));
0751       }
0752     }
0753   }
0754   return ga;
0755 }
0756
0757
0758
0759
0760
0761
0762 //**********************************************************************
0763 long double log_gamma_function(long double x)
0764 // Return ln(Gamma(x)).
0765 // Courtesy to Codecog.org.
0766 // This funcion is under GP license.
0767 {
0768   if (x <= 0 || x > XBIG)
0769     return HUGE_VAL;
0770
0771   if (x <= MACHINE_EPSILON)
0772     return -log(x);
0773
0774   if (x <= 4)
0775   {
0776     double
0777     p1[8] =
0778     {
0779       4.945235359296727046734888E+00, 2.018112620856775083915565E+02, 2.290838373831346393026739E+03,
0780       1.131967205903380828685045E+04, 2.855724635671635335736389E+04, 3.848496228443793359990269E+04,
0781       2.637748787624195437963534E+04, 7.225813979700288197698961E+03
0782     },
0783     q1[8] =
0784     {
0785       6.748212550303777196073036E+01, 1.113332393857199323513008E+03, 7.738757056935398733233834E+03,
0786       2.763987074403340708898585E+04, 5.499310206226157329794414E+04, 6.161122180066002127833352E+04,
0787       3.635127591501940507276287E+04, 8.785536302431013170870835E+03
0788     },
0789     p2[8] =
0790     {
0791       4.974607845568932035012064E+00, 5.424138599891070494101986E+02, 1.550693864978364947665077E+04,
0792       1.847932904445632425417223E+05, 1.088204769468828767498470E+06, 3.338152967987029735917223E+06,
0793       5.106661678927352456275255E+06, 3.074109054850539556250927E+06
0794     },
0795     q2[8] =
0796     {
0797       1.830328399370592604055942E+02, 7.765049321445005871323047E+03, 1.331903827966074194402448E+05,
0798       1.136705821321969608938755E+06, 5.267964117437946917577538E+06, 1.346701454311101692290052E+07,
0799       1.782736530353274213975932E+07, 9.533095591844353613395747E+06
0800     };
0801
0802     double corr = x >= PNT68 ? 0 : -log(x);
0803     double xden = 1, xnum = 0, xm = x <= 1.5 ? (x > 0.5 ? x - 1 : x) : x - 2;
0804     bool flag = false;
0805     if (x <= 1.5 && (x <= 0.5 || x >= PNT68)) flag = true;
0806     double *p = flag ? p1 : p2, *q = flag ? q1 : q2;
0807
0808     xnum = xnum * xm + p[0], xden = xden * xm + q[0];
0809     xnum = xnum * xm + p[1], xden = xden * xm + q[1];
0810     xnum = xnum * xm + p[2], xden = xden * xm + q[2];
0811     xnum = xnum * xm + p[3], xden = xden * xm + q[3];
0812     xnum = xnum * xm + p[4], xden = xden * xm + q[4];
0813     xnum = xnum * xm + p[5], xden = xden * xm + q[5];
0814     xnum = xnum * xm + p[6], xden = xden * xm + q[6];
0815     xnum = xnum * xm + p[7], xden = xden * xm + q[7];
0816
0817     return (x > 1.5 ? 0 : corr) + xm * ((flag ? D1 : D2) + xm * (xnum / xden));
0818   }
0819
0820   if (x <= 12)
0821   {
0822     double xm = x - 4, xden = -1, xnum = 0,
0823     p[8] =
0824     {
0825       1.474502166059939948905062E+04, 2.426813369486704502836312E+06, 1.214755574045093227939592E+08,
0826       2.663432449630976949898078E+09, 2.940378956634553899906876E+010,1.702665737765398868392998E+011,
0827       4.926125793377430887588120E+011, 5.606251856223951465078242E+011
0828     },
0829     q[8] =
0830     {
0831       2.690530175870899333379843E+03, 6.393885654300092398984238E+05, 4.135599930241388052042842E+07,
0832       1.120872109616147941376570E+09, 1.488613728678813811542398E+010, 1.016803586272438228077304E+011,
0833       3.417476345507377132798597E+011, 4.463158187419713286462081E+011
0834     };
0835
0836     xnum = xnum * xm + p[0], xden = xden * xm + q[0];
0837     xnum = xnum * xm + p[1], xden = xden * xm + q[1];
0838     xnum = xnum * xm + p[2], xden = xden * xm + q[2];
0839     xnum = xnum * xm + p[3], xden = xden * xm + q[3];
0840     xnum = xnum * xm + p[4], xden = xden * xm + q[4];
0841     xnum = xnum * xm + p[5], xden = xden * xm + q[5];
0842     xnum = xnum * xm + p[6], xden = xden * xm + q[6];
0843     xnum = xnum * xm + p[7], xden = xden * xm + q[7];
0844
0845     return D4 + xm * (xnum / xden);
0846   }
0847
0848   double res = 0;
0849   if (x <= FRTBIG)
0850   {
0851     double xsq = x * x,
0852     c[7] =
0853     {
0854       -1.910444077728E-03, 8.4171387781295E-04, -5.952379913043012E-04, 7.93650793500350248E-04,
0855       -2.777777777777681622553E-03, 8.333333333333333331554247E-02, 5.7083835261E-03
0856     };
0857     res = c[6];
0858     res = res / xsq + c[0];
0859     res = res / xsq + c[1];
0860     res = res / xsq + c[2];
0861     res = res / xsq + c[3];
0862     res = res / xsq + c[4];
0863     res = res / xsq + c[5];
0864   }
0865   double corr = log(x);
0866   res /= x, res += SQRTPI - corr / 2 + x * (corr - 1);
0867   return res;
0868 }
0869
0870
0871
0872
0873
0874 //**********************************************************************
0875 double beta_function(double x, double y)
0876 // B(x,y):=Gamma(x)Gamma(y)/Gamma(x+y)
0877 {
0878     return exp(log_gamma_function(x)+log_gamma_function(y)-log_gamma_function(x+y));
0879 }
0880
0881
0882
0883
0884
0885 //**********************************************************************
0886 double binomial_coefficient(double n, double k)
0887 // Returns the binomial coefficient
0888 // C_n^k := Gamma(n+1) / (Gamma(k+1)*Gamma(n-k+1))
0889 {
0890     return 1.0/(n+1)/beta_function(k+1,n-k+1);
0891 }
0892
0893
0894
0895
0896
0897 //**********************************************************************
0898 void print_progressbar(double percentage, int length, string symbol)
0899 // Print out a progress bar with the given percentage. Use a negative value to reset the progress bar.
0900 {
0901   static int status=0;
0902   static int previous_stop=0;
0903
0904   if (percentage<0)
0905   {
0906     // reset
0907     status = 0;
0908     previous_stop = 0;
0909   }
0910
0911   // initializing progressbar
0912   if (status==0)
0913   {
0914     cout << "\r";
0915     cout << "[";
0916     for (int i=1; i<=length; i++) cout << " ";
0917     cout << "]";
0918     cout << "\r";
0919     cout << "[";
0920   }
0921
0922   // plot status
0923   int stop;
0924   if (percentage>=0.99) stop=0.99*length;
0925   else stop = percentage*length;
0926   for (int i=previous_stop; i<stop; i++) cout << symbol;
0927   if (previous_stop<stop) previous_stop=stop;
0928
0929   // simulate a rotating bar
0930   if (status==0) cout << "-";
0931   switch (status)
0932   {
0933     case 1: cout << "\\"; break;
0934     case 2: cout << "|"; break;
0935     case 3: cout << "/"; break;
0936     case 4: cout << "-"; break;
0937   }
0938   cout << "\b";
0939   status++;
0940   if (status==5) status=1;
0941   cout.flush();
0942 }
0943
0944
0945 //**********************************************************************
0946 void formatedPrint(ostream& os, int count, ...)
0947 // For easier scientific data outputing.
0948 {
0949   va_list ap;
0950   va_start(ap, count); //Requires the last fixed parameter (to get the address)
0951   for(int j=0; j<count; j++)
0952       os << scientific << setprecision(OUTPUT_PRECISION) << "  " << va_arg(ap, double); //Requires the type to cast to. Increments ap to the next argument.
0953   va_end(ap);
0954   os << endl;
0955 }
0956
0957
0958
0959 //**********************************************************************
0960 void display_logo(int which)
0961 // Personal amusement.
0962 {
0963   switch (which)
0964   {
0965     case 1:
0966     cout << " ____  ____            _                    " << endl;
0967     cout << "|_   ||   _|          (_)                   " << endl;
0968     cout << "  | |__| |    .---.   __    _ .--.    ____  " << endl;
0969     cout << "  |  __  |   / /__\\\\ [  |  [ `.-. |  [_   ] " << endl;
0970     cout << " _| |  | |_  | \\__.,  | |   | | | |   .' /_ " << endl;
0971     cout << "|____||____|  '.__.' [___] [___||__] [_____]" << endl;
0972     cout << "                                            " << endl;
0973     break;
0974
0975     case 2:
0976     cout << ":::    ::: :::::::::: ::::::::::: ::::    ::: :::::::::" << endl;
0977     cout << ":+:    :+: :+:            :+:     :+:+:   :+:      :+: " << endl;
0978     cout << "+:+    +:+ +:+            +:+     :+:+:+  +:+     +:+  " << endl;
0979     cout << "+#++:++#++ +#++:++#       +#+     +#+ +:+ +#+    +#+   " << endl;
0980     cout << "+#+    +#+ +#+            +#+     +#+  +#+#+#   +#+    " << endl;
0981     cout << "#+#    #+# #+#            #+#     #+#   #+#+#  #+#     " << endl;
0982     cout << "###    ### ########## ########### ###    #### #########" << endl;
0983     break;
0984
0985     case 3:
0986     cout << " __  __     ______     __     __   __     _____    " << endl;
0987     cout << "/\\ \\_\\ \\   /\\  ___\\   /\\ \\   /\\ '-.\\ \\   /\\___  \\  " << endl;
0988     cout << "\\ \\  __ \\  \\ \\  __\\   \\ \\ \\  \\ \\ \\-.  \\  \\/_/  /__ " << endl;
0989     cout << " \\ \\_\\ \\_\\  \\ \\_____\\  \\ \\_\\  \\ \\_\\\\'\\_\\   /\\_____\\" << endl;
0990     cout << "  \\/_/\\/_/   \\/_____/   \\/_/   \\/_/ \\/_/   \\/_____/" << endl;
0991     break;
0992
0993   }
0994
0995 }
0996
0997
0998
0999 //**********************************************************************
1000 void GaussLegendre_getWeight(int npts,double* xg,double* wg, double A, double B, int iop)
1001 // Calculate the sampling location and weight for Gauss-Legendre quadrature
1002 // -- From Hirano and Nara's MC-KLN code.
1003 {
1004 //ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1005 //  gauss.f: Points and weights for Gaussian quadrature                 c
1006 //                       c
1007 //  taken from: "Projects in Computational Physics" by Landau and Paez  c
1008 //         copyrighted by John Wiley and Sons, New York            c
1009 //                                                                      c
1010 //  written by: Oregon State University Nuclear Theory Group            c
1011 //         Guangliang He & Rubin H. Landau                         c
1012 //  supported by: US National Science Foundation, Northwest Alliance    c
1013 //                for Computational Science and Engineering (NACSE),    c
1014 //                US Department of Energy                          c
1015 //                       c
1016 //  comment: error message occurs if subroutine called without a main   c
1017 //  comment: this file has to reside in the same directory as integ.c   c
1018 //ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1019
1020   static const double EPS = 3.0e-14;
1021   int m=(npts+1)/2;
1022   for(int i=0;i<m;i++) {
1023     double  t=cos(M_PI*(i+1-0.25)/(npts+0.5));
1024     double t1=t;
1025     double pp;
1026     do {
1027         double p1=1.0;
1028         double p2=0.0;
1029         double aj=0.0;
1030         for(int j=0;j<npts;j++) {
1031             double p3=p2;
1032             p2=p1;
1033             aj += 1.0;
1034             p1=((2.0*aj-1.0)*t*p2-(aj-1.0)*p3)/aj;
1035         }
1036         pp=npts*(t*p1-p2)/(t*t-1.0);
1037         t1=t;
1038         t=t1-p1/pp;
1039     } while(abs(t-t1)>EPS);
1040     xg[i]=-t;
1041     xg[npts-1-i]=t;
1042     wg[i]=2.0/((1.0-t*t)*pp*pp);
1043     wg[npts-1-i]=wg[i];
1044     }
1045
1046 //GaussLegendre::GaussRange(int N,int iop,double A,double B,
1047 //  double* xg1,double* wg1)
1048 //{
1049 //     transform gausspoints to other range than [-1;1]
1050 //     iop = 1  [A,B]       uniform
1051 //     iop = 2  [0,inf]     A is midpoint
1052 //     opt = 3  [-inf,inf]  scale is A
1053 //     opt = 4  [B,inf]     A+2B is midoint
1054 //     opt = 5  [0,B]     AB/(A+B)+ is midoint
1055
1056   int N=npts;
1057   double xp, wp;
1058   for(int i=0; i<N; i++) {
1059       if(iop == 1) {
1060           //...... A to B
1061           xp=(B+A)/2+(B-A)/2*xg[i];
1062           wp=(B-A)/2*wg[i];
1063       } else if(iop == -1) {
1064           //......   A to B
1065           xp=(B+A)/2+(B-A)/2*xg[i];
1066           if(i <= N/2)
1067             xp=(A+B)/2-(xp-A);
1068           else
1069             xp=(A+B)/2+(B-xp);
1070           wp=(B-A)/2*wg[i];
1071       } else if(iop == 2) {
1072           //...... zero to infinity
1073           xp=A*(1+xg[i])/(1-xg[i]);
1074           double tmp=(1-xg[i]);
1075           wp=2.*A/(tmp*tmp)*wg[i];
1076       } else if(iop ==  3) {
1077           //...... -inf to inf scale A
1078           xp=A*(xg[i])/(1-xg[i]*xg[i]);
1079           double tmp=1-xg[i]*xg[i];
1080           wp=A*(1+xg[i]*xg[i])/(tmp*tmp)*wg[i];
1081       } else if(iop == 4) {
1082           //......  B to inf,  A+2B is midoint
1083           xp=(A+2*B+A*xg[i])/(1-xg[i]);
1084           double tmp=1-xg[i];
1085           wp=2.*(B+A)/(tmp*tmp)*wg[i];
1086       } else if(iop == 5) {
1087           //...... -A to A , scale B
1088           //xp=A*pow(abs(xg[i]),B) *sign(1.0,xg(i));
1089           double tmp = xg[i] >= 0 ? 1.0 : -1.0;
1090           xp=A*pow(abs(xg[i]),B) * tmp;
1091           //xp=A*pow(abs(xg[i]),B) *sign(1.0,xg(i));
1092           wp=A*B*pow(abs(xg[i]),(B-1))*wg[i];
1093       } else if(iop ==  6) {
1094           //...... 0 to B , AB/(A+B) is midpoint
1095           xp=A*B*(1+xg[i])/(B+A-(B-A)*xg[i]);
1096           double tmp = B+A-(B-A)*xg[i];
1097           wp=2*A*B*B/(tmp*tmp)*wg[i];
1098       } else {
1099           cerr << " invalid option iop = " << iop << endl;
1100           exit(-1);
1101       }
1102       xg[i]=xp;
1103       wg[i]=wp;
1104   }
1105 }
1106
1107
1108
1109
1110
1111 //**********************************************************************
1112 void get_bin_average_and_count(istream& is, ostream& os, vector<double>* bins, long col_to_bin, void (*func)(vector<double>*), long wanted_data_columns, bool silence)
1113 // Group data into bins by set by the "bins". The data in the column
1114 // "col_to_bin" read from "is" are the ones used to determine the binning.
1115 // Once the binning is decided, the averages of all data are calculated.
1116 // The result, which has the same number of rows and the number of bins,
1117 // will be outputted to "os". The i-th line of the output data contains the
1118 // average of data from each column, for the i-th bin. The output will
1119 // have 2 more columns; the 1st being the number count N and the 2nd being
1120 // dN/dx where dx is the bin width.
1121 // The values given in "bins" define the boundaries of bins and is assumed
1122 // to be non-decreasing.
1123 // After each line is decided to go into which bin, the function specified
1124 // by "func" will be called to transform data. It is the transformed data
1125 // that will be averaged. The transformed data can have different number of
1126 // columns than the data passed in, in which case its number of columns
1127 // is specified by "wanted_data_columns". The counting info will still
1128 // be always written in the last two columns.
1129 // The function specified by "func" should accepts a vector of doubles
1130 // which is one line of data, and then modify it as returned result. The
1131 // data vector passed in already has the desired length so it can be modified
1132 // directly.
1133 // The argument "col_to_bin" starts with 1.
1134 // Refer also to getBinnedAverageAndCount MATLAB program.
1135 {
1136   // initialization
1137   char* buffer = new char[99999]; // each line should be shorter than this
1138   // get first line and continue initialization
1139   is.getline(buffer, 99999);
1140   vector<double> line_data = stringToDoubles(buffer);
1141   long number_of_cols = line_data.size();
1142   long number_of_bins = bins->size()-1;
1143
1144   if (number_of_cols==0)
1145   {
1146     cout << "get_bin_average_and_count error: the input data is empty!" << endl;
1147     exit(-1);
1148   }
1149
1150   // create the counting array
1151   if (wanted_data_columns>0) number_of_cols = wanted_data_columns;
1152   double bin_total_and_count[number_of_bins][number_of_cols+2];
1153   for (long i=0; i<number_of_bins; i++)
1154   for (long j=0; j<number_of_cols+2; j++)
1155   {
1156     bin_total_and_count[i][j] = 0;
1157   }
1158
1159   // add up all data
1160   long number_of_lines=1;
1161   while (is.eof()==false)
1162   {
1163     // determine which bin
1164     long bin_idx = binarySearch(bins, line_data[col_to_bin-1], true);
1165     if (bin_idx!=-1)
1166     {
1167       // transform data
1168       line_data.resize(number_of_cols, 0);
1169       if (func) (*func) (&line_data);
1170       // add to the counting matrix
1171       for (long j=0; j<number_of_cols; j++) bin_total_and_count[bin_idx][j] += line_data[j];
1172       // also the counting column
1173       bin_total_and_count[bin_idx][number_of_cols] ++;
1174     }
1175     // next iteration
1176     is.getline(buffer, 99999);
1177     line_data = stringToDoubles(buffer);
1178     if (number_of_lines % 100000 == 0 && !silence) cout << "Line " << number_of_lines << " reached." << endl;
1179     number_of_lines++;
1180   }
1181
1182   // find the averages
1183   for (long i=0; i<number_of_bins; i++)
1184   for (long j=0; j<number_of_cols; j++)
1185   {
1186     if (bin_total_and_count[i][number_of_cols]<1e-15) continue;
1187     bin_total_and_count[i][j] /= bin_total_and_count[i][number_of_cols];
1188   }
1189
1190
1191   // get dN/(d bin_width)
1192   for (long i=0; i<number_of_bins; i++)
1193   {
1194     if (bin_total_and_count[i][number_of_cols]<1e-15) continue;
1195     bin_total_and_count[i][number_of_cols+1] = bin_total_and_count[i][number_of_cols]/((*bins)[i+1]-(*bins)[i]);
1196   }
1197
1198   // output
1199   for (long i=0; i<number_of_bins; i++)
1200   {
1201     for (long j=0; j<number_of_cols+2; j++)
1202     {
1203       os << scientific << scientific << setprecision(OUTPUT_PRECISION) << bin_total_and_count[i][j] << "  ";
1204     }
1205     os << endl;
1206   }
1207
1208 }