//a routine for the plotting the radial distribution function g(r) from the output of the Verlet.exe executable. This code can be run directly inside of ROOT, which is a plotting package distributed by CERN, at http://root.cern.ch/
{

  gROOT->Reset();

  const int NUMBER_OF_STEPS = 65536;
  const int GOFR_FREQUENCY = 8;
  const int NUMBER_OF_BINS = 256;
  const double DELTA_R = 1.03328056486624947e-2;

  ifstream radialFile( "RadialDistribution.dat", ios_base::in ); 
  
  double dummyDouble;

  for ( int i=0; i<NUMBER_OF_STEPS/GOFR_FREQUENCY-1; i++ ) {
    for ( int j=0; j<NUMBER_OF_BINS; j++ ) {
      radialFile >> dummyDouble;
      //cout << dummyDouble << endl;
    }
  }

  double gOfR[NUMBER_OF_BINS];
  double radius[NUMBER_OF_BINS];

  for ( int i=0; i<NUMBER_OF_BINS; i++ ) {
     radialFile >> gOfR[i];
     radius[i] = (0.5+((double)(i)))*DELTA_R;
     cout << gOfR[i] << "  ";
  }

  cout << endl;

  radialFile.close();

  TCanvas * c1 = new TCanvas( "c1", "Radial Distribution Function" );
  TGraph * radialGraph = new TGraph( NUMBER_OF_BINS, radius, gOfR );

  radialGraph->SetTitle("Radial Distribution Function");
  (radialGraph->GetXaxis())->SetTitle("r");
  (radialGraph->GetYaxis())->SetTitle("g(r)");
  radialGraph->Draw("ALP");

}