source: util/average.cpp@ 774ae8

using namespace std;
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>


int main(int argc, char **argv) {

        double avg, dev;
        double tmp;
        int i, j, zaehler;
  int firstcol, lastcol;
        double ecut;
        int flag;
        string zeile;
        
        if (argc < 4) {
                cout << "Usage: " << argv[0] << " <firstCol> <lastCol> <X> <File>\n";
                return(1);
        } 
        firstcol = atoi(argv[1]);
  lastcol = atoi(argv[2]);
  ecut = atof(argv[3]);

        // get average
        avg = 0.;
        zaehler=0;
        ifstream test(argv[4]);
        if (test == NULL) { cout << "Can't open File " << argv[4] << "\n"; return(255); }
        flag=1;
        //cout << "Looking for " << ecut << " \n";
        while (getline(test, zeile, '\n')) {
                //cout << zeile;
                istringstream input(zeile);
                input >> ws >> tmp;                     
                if (tmp==(double)ecut) {        // found correct line!
                        for (j=2;j<=lastcol;j++)
                                if (!input.eof()) {
                                        input >> ws >> tmp;
          if (j >= firstcol) {
            avg += tmp;
            zaehler++;
          }
        } 
                } 
        }
  test.clear();
        test.seekg(ios::beg);
        if (zaehler != 0) avg /= zaehler;       

        // get deviation
        dev = 0.;
        zaehler=0;
        flag=1;
        while (getline(test, zeile, '\n')) {
                istringstream input(zeile);
                input >> ws >> tmp;                     
                if (tmp==(double)ecut) {        // found correct line!
                        for (j=2;j<=lastcol;j++)
                                if (!input.eof()) {
                                        input >> ws >> tmp;
          if (j >= firstcol) {
            dev += (tmp - avg)*(tmp - avg);
            zaehler++;
          }
        }
                }
        }
        test.close();
        //dev = 1/(zaehler)*dev;
        if (dev != 0) dev /= zaehler;
                        
        cout << setprecision(8) << avg << "\t" << setprecision(8) << sqrt(dev) << "\n";
        return(0);
}