/** \file analyzer.cpp
 *
 * Takes evaluated fragments (energy and forces) and does evaluation of how sensible the BOSSANOVA
 * approach was, e.g. in the decay of the many-body-contributions.
 *    
 */

//============================ INCLUDES ===========================

#include "datacreator.hpp"
#include "helpers.hpp" 
#include "parser.hpp"
#include "periodentafel.hpp" 

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif


//============================== MAIN =============================

int main(int argc, char **argv)
{
  periodentafel *periode = NULL; // and a period table of all elements
  EnergyMatrix Energy;
  ForceMatrix Force;
  ForceMatrix Shielding;
  ForceMatrix ShieldingPAS;
  MatrixContainer Time;
  EnergyMatrix EnergyFragments;
  ForceMatrix ForceFragments;
  ForceMatrix ShieldingFragments;
  ForceMatrix ShieldingPASFragments;
  KeySetsContainer KeySet;
  ofstream output;
  ofstream output2;
  ofstream output3;
  ofstream output4;
  ifstream input;
  stringstream filename;
  time_t t = time(NULL);
  struct tm *ts = localtime(&t);
  char *datum = asctime(ts);
  stringstream Orderxrange;
  stringstream Fragmentxrange;

  cout << "ANOVA Analyzer" << endl;
  cout << "==============" << endl;
  
  // Get the command line options
  if (argc < 4) {
    cout << "Usage: " << argv[0] << " <inputdir> <prefix> <outputdir> [elementsdb]" << endl;
    cout << "<inputdir>\ttherein the output of a molecuilder fragmentation is expected, each fragment with a subdir containing an energy.all and a forces.all file." << endl;
    cout << "<prefix>\tprefix of energy and forces file." << endl;
    cout << "<outputdir>\tcreated plotfiles and datafiles are placed into this directory " << endl;
    cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
    return 1;
  }
  if (argc > 4) {
    cout << "Loading periodentafel." << endl;
    periode = new periodentafel;
    periode->LoadPeriodentafel(argv[4]);
  }
  
  // Test the given directory
  if (!TestParams(argc, argv))
    return 1;

  // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++
  
  // ------------- Parse through all Fragment subdirs --------
  if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix,0,0)) return 1;
  if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix,0,0)) return 1;
  if (!Time.ParseMatrix(argv[1], argv[2], TimeSuffix, 10,1)) return 1;
  if (periode != NULL) { // also look for PAS values
    if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
    if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
  }

  // ---------- Parse the TE Factors into an array -----------------
  if (!Energy.ParseIndices()) return 1;
  
  // ---------- Parse the Force indices into an array ---------------
  if (!Force.ParseIndices(argv[1])) return 1;

  // ---------- Parse the shielding indices into an array ---------------
  if (periode != NULL) { // also look for PAS values
    if(!Shielding.ParseIndices(argv[1])) return 1;
    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
  }

  // ---------- Parse the KeySets into an array ---------------
  if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
  if (!KeySet.ParseManyBodyTerms()) return 1;
  
  // ---------- Parse fragment files created by 'joiner' into an array -------------
  if (!EnergyFragments.ParseMatrix(argv[1], argv[2], EnergyFragmentSuffix,0,0)) return 1;
  if (!ForceFragments.ParseMatrix(argv[1], argv[2], ForceFragmentSuffix,0,0)) return 1;
  if (periode != NULL) { // also look for PAS values
    if (!ShieldingFragments.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
    if (!ShieldingPASFragments.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
  }

  // +++++++++++++++ TESTING ++++++++++++++++++++++++++++++
  
  // print energy and forces to file
  filename.str("");
  filename << argv[3] << "/" << "energy-forces.all";
  output.open(filename.str().c_str(), ios::out);
  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
    for(int k=0;k<Energy.ColumnCounter;k++)
      output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
    output << endl;
  }
  output << endl;
  
  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    for(int k=0;k<Force.ColumnCounter;k++)
      output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    output << endl;
  }
  output << endl;

  if (periode != NULL) { // also look for PAS values
    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
      for(int k=0;k<Shielding.ColumnCounter;k++)
        output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;
  
    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
      for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
      output << endl;
    }
    output << endl;
  }
  
  output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
  Time.SetLastMatrix(0., 0);
  for (int BondOrder=KeySet.Order;BondOrder--;)
    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
        for(int k=Time.ColumnCounter;k--;) {
          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
        }
  for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
    for(int k=0;k<Time.ColumnCounter;k++)
      output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
    output << endl;
  }
  output << endl;
  output.close();

  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
  
  cout << "Analyzing ..." << endl;

  // ======================================= Creating the data files ==============================================================

  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
  if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
  Time.SetLastMatrix(0., 0);
  output << "#Order\tFrag.No.\t" << Time.Header << endl;
  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
        for(int k=Time.ColumnCounter;k--;) {
          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
        }
    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    for(int k=0;k<Time.ColumnCounter;k++)
      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
    output << endl;
  }
  output.close();

  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
  if (!EnergyFragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0)) return 1;
  if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "DeltaEnergies-Order", "Plot of error between approximated and full energies energies versus the Bond Order", datum)) return 1;
  
  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
  if (!EnergyFragments.SetLastMatrix(1.,0)) return 1;
  if (!CreateDataEnergyOrder(Energy, EnergyFragments, KeySet, argv[3], "Energies-Order", "Plot of approximated energies versus the Bond Order", datum)) return 1;
  
  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
  if (!ForceFragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0)) return 1;
  if (!OpenOutputFile(output, argv[3], "DeltaForces-Order.dat" )) return false;
  if (!OpenOutputFile(output2, argv[3], "DeltaMinForces-Order.dat" )) return false;
  if (!OpenOutputFile(output3, argv[3], "DeltaMeanForces-Order.dat" )) return false;
  if (!OpenOutputFile(output4, argv[3], "DeltaMaxForces-Order.dat" )) return false;
  cout << "Plot of per atom and min/mean/max error between approximated forces and full forces versus the Bond Order ... " << endl;
  output << "# Plot of error between approximated forces and full forces versus the Bond Order, created on " << datum;
  output << "# AtomNo" << Force.Header << endl;
  output2 << "# Plot of min error between approximated forces and full forces versus the Bond Order, created on " << datum;
  output2 << "# Order\tFrag.No.\t" << Force.Header << endl;
  output3 << "# Plot of mean error between approximated forces and full forces versus the Bond Order, created on " << datum;
  output3 << "# Order\tFrag.No.\t" << Force.Header << endl;
  output4 << "# Plot of max error between approximated forces and full forces versus the Bond Order, created on " << datum;
  output4 << "# Order\tFrag.No.\t" << Force.Header << endl;
  int FragmentCounter = 0;
  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
        if (j > Force.RowCounter[Force.MatrixCounter]) {
          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
          return 1;
        }
        if (j != -1)
          for(int k=0;k<Force.ColumnCounter;k++) {
            Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
          }
      }
      FragmentCounter++;
    }
    // errors per atom
    output << "#Order\t" << BondOrder+1 << endl;
    for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) { 
      output << Force.Indices[Force.MatrixCounter][j] << "\t";
      for (int l=0;l<Force.ColumnCounter;l++)
        if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) < MYEPSILON)
          output << scientific << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
        else
          output << scientific << (Force.Matrix[Force.MatrixCounter][ j ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ j ][l]) << "\t";
      output << endl;
    }
    output << endl;
    
    // min error
    output2 << BondOrder+1 << "\t" << FragmentCounter;
    CreateMinimumForce(Force, Force.MatrixCounter);
    for (int l=0;l<Force.ColumnCounter;l++)
      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
        output2 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
      else
        output2 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
    output2 << endl;
    
    // mean error
    output3 << BondOrder+1 << "\t" << FragmentCounter;
    CreateMeanForce(Force, Force.MatrixCounter);
    for (int l=0;l<Force.ColumnCounter;l++)
      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
        output3 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
      else
        output3 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
    output3 << endl;
    
    // max error
    output4 << BondOrder+1 << "\t" << FragmentCounter;
    CreateMaximumForce(Force, Force.MatrixCounter);
    for (int l=0;l<Force.ColumnCounter;l++)
      if (fabs(ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]) < MYEPSILON)
        output4 << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
      else
        output4 << scientific << "\t" << (Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l] / ForceFragments.Matrix[ForceFragments.MatrixCounter][ ForceFragments.RowCounter[ForceFragments.MatrixCounter] ][l]);
    output4 << endl;
  }
  output.close();
  output2.close();
  output3.close();
  output4.close();

  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
  if(!OpenOutputFile(output, argv[3], "Forces-Order.dat")) return 1;
  cout << "Plot of approximated forces versus the Bond Order ... " << endl;
  output << "# Plot of approximated forces versus the Bond Order, created on " << datum;
  output << "# AtomNo" << Force.Header << endl;
  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
      for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
        int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
        if (j > Force.RowCounter[Force.MatrixCounter]) {
          cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
          return 1;
        }
        if (j != -1)
          for(int k=0;k<Force.ColumnCounter;k++) {
            Force.Matrix[Force.MatrixCounter][j][k] -= ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
          }
      }
    }
    // errors per atom
    output << "#Order\t" << BondOrder+1 << endl;
    for(int j=0;j<Force.RowCounter[ Force.MatrixCounter ];j++) { 
      output << Force.Indices[Force.MatrixCounter][j] << "\t";
      for (int l=0;l<Force.ColumnCounter;l++)
         output << Force.Matrix[Force.MatrixCounter][ j ][l] << "\t";
      output << endl;
    }
    output << endl;
  }
  output.close();

  // min force
  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MeanForces-Order", "Plot of mean approximated forces versus the Bond Order", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "MaxForces-Order", "Plot of max approximated forces versus the Bond Order", datum, CreateMaximumForce)) return 1;

  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum (should be 0) vs. bond order
  if (!CreateDataForcesOrder(Force, ForceFragments, KeySet, argv[3], "VectorSum-Order", "Plot of vector sum of the approximated forces versus the Bond Order", datum, CreateVectorSumForce)) return 1;

  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-Fragment", "Plot of fragment energy versus the Fragment No", datum, CreateEnergy)) return 1;
  if (!CreateDataFragment(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", "Plot of fragment energy of each Fragment No vs. Bond Order", datum, CreateEnergy)) return 1;
  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", "Plot of maximum of fragment energy vs. Bond Order", datum, CreateMaxFragmentOrder)) return 1;
  if (!CreateDataFragmentOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", "Plot of minimum of fragment energy vs. Bond Order", datum, CreateMinFragmentOrder)) return 1;

  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment
  // min force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-Fragment", "Plot of min approximated forces versus the Fragment No", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", "Plot of mean approximated forces versus the Fragment No", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", "Plot of max approximated forces versus the Fragment No", datum, CreateMaximumForce)) return 1;

  // +++++++++++++++++++++++++++++++Ploting min/mean/max forces for each fragment per order
  // min force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", "Plot of min approximated forces of each Fragment No vs. Bond Order", datum, CreateMinimumForce)) return 1;

  // mean force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", "Plot of mean approximated forces of each Fragment No vs. Bond Order", datum, CreateMeanForce)) return 1;

  // max force
  if (!CreateDataFragment(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", "Plot of max approximated forces of each Fragment No vs. Bond Order", datum, CreateMaximumForce)) return 1;

  // ======================================= Creating the plot files ==============================================================
 
  Orderxrange << "[1:" << KeySet.Order << "]";
  Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
  
  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
  if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
  
  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
  if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), "[1e-8:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
  
  // +++++++++++++++++++++++++++++++++++Plotting Energies vs. Order
  if (!CreatePlotOrder(Energy, KeySet, argv[3], "Energies-Order", 1, "outside", "", "",  1, 1, "bond order k", "approximate energy [Ht]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;

  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in forces to full QM
  // min force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMinForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in min force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // mean force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMeanForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in mean force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;

  // max force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "DeltaMaxForces-Order", 2, "top right", "y", "",  1, 1, "bond order k", "absolute error in max force [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-8:1e+0]", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // min/mean/max comparison for total force
  if(!OpenOutputFile(output, argv[3], "DeltaMinMeanMaxTotalForce-Order.pyx")) return 1;
  CreatePlotHeader(output, "DeltaMinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL,  1, 1, "bond order k", "absolute error in total forces [Ht/a.u.]");
  output << "plot " << Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;  
  output << "'DeltaMinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
  output << "'DeltaMeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
  output << "'DeltaMaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
  output.close();  

  // ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
  // min force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "MinForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated min force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // mean force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "MeanForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated mean force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", AbsFirstForceValuePlotLine)) return 1;

  // max force
  if (!CreatePlotOrder(Force, KeySet, argv[3], "MaxForces-Order", 2, "bottom right", "y", "",  1, 1, "bond order k", "absolute approximated max force [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // min/mean/max comparison for total force
  if(!OpenOutputFile(output, argv[3],"MinMeanMaxTotalForce-Order.pyx")) return 1;
  CreatePlotHeader(output, "MinMeanMaxTotalForce-Order", 1, "bottom left", "y", NULL, 1, 1, "bond order k", "absolute total force [Ht/a.u.]");
  output << "plot "<< Orderxrange.str().c_str() << " [1e-8:1e+0] \\" << endl;  
  output << "'MinForces-Order.dat' title 'minimum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints, \\" << endl;
  output << "'MeanForces-Order.dat' title 'mean' using 1:(abs($" << 8 << ")) with linespoints, \\" << endl;
  output << "'MaxForces-Order.dat' title 'maximum' using 1:(sqrt($" << 8 << "*$" << 8 << "+$" << 8+1 << "*$" << 8+1 << "+$" << 8+2 << "*$" << 8+2 << ")) with linespoints" << endl;
  output.close();  

  // ++++++++++++++++++++++++++++++++++++++Plotting vector sum vs. Order

  if (!CreatePlotOrder(Force, KeySet, argv[3], "VectorSum-Order", 2, "bottom right", "y" ,"", 1, 1, "bond order k", "vector sum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "", "1" , "with linespoints", ForceMagnitudePlotLine)) return 1;

  // +++++++++++++++++++++++++++++++Plotting energyfragments vs. order
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with points", AbsEnergyPlotLine)) return 1;
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Energies of each fragment [Ht]", Orderxrange.str().c_str(), "[1e-10:1e+1]", "1" , "with points", AbsEnergyPlotLine)) return 1;
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MaxEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Maximum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
  if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "MinEnergies-FragmentOrder", 5, "below", "y", "", 1, 1, "bond order", "Minimum fragment energy [Ht]", Orderxrange.str().c_str(), "[1e-6:1e+1]", "1" , "with linespoints", AbsEnergyPlotLine)) return 1;

  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
  // min
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "minimum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
    
  // mean
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "mean of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
    
  // max
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-Fragment", 5, "below", "y", "set boxwidth 0.2", 1, 5, "fragment number", "maximum of approximated forces [Ht/a.u.]", Fragmentxrange.str().c_str(), "[1e-16:1e+1]", "2" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;

  // +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment per bond order
  // min
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MinForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "minimum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
    
  // mean
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MeanForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "mean of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesFirstForceValuePlotLine)) return 1;
    
  // max
  if (!CreatePlotOrder(ForceFragments, KeySet, argv[3], "MaxForces-FragmentOrder", 5, "below", "y", "set boxwidth 0.2", 1, 1, "bond order", "maximum of approximated forces [Ht/a.u.]", Orderxrange.str().c_str(), "[1e-16:1e+1]", "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;
 
  // create Makefile
  if(!OpenOutputFile(output, argv[3], "Makefile")) return 1;
  output << "PYX = $(shell ls *.pyx)" << endl << endl;
  output << "EPS = $(PYX:.pyx=.eps)" << endl << endl;
  output << "%.eps: %.pyx" << endl;
  output << "\t~/build/pyxplot/pyxplot $<" << endl << endl;
  output << "all: $(EPS)" << endl << endl;
  output << ".PHONY: clean" << endl;
  output << "clean:" << endl;
  output << "\trm -rf $(EPS)" << endl;
  output.close();
 
  // ++++++++++++++++ exit ++++++++++++++++++++++++++++++++++
  delete(periode);
  cout << "done." << endl;
  return 0;
};

//============================ END ===========================