/** \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;
	EnergyMatrix Hcorrection;
	ForceMatrix Force;
	ForceMatrix Shielding;
	ForceMatrix ShieldingPAS;
	EnergyMatrix Time;
	EnergyMatrix EnergyFragments;
	EnergyMatrix HcorrectionFragments;
	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;
	stringstream yrange;
	char *dir = NULL;
	bool Hcorrected = true;
	double norm;
	int counter;

	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;
	} else {
		dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
		strcpy(dir, "/");
		strcat(dir, argv[2]);
	}

	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.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
	Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
	if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
	if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
	if (periode != NULL) { // also look for PAS values
		if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
		if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
	}

	// ---------- Parse the TE Factors into an array -----------------
	if (!Energy.ParseIndices()) return 1;
	if (Hcorrected) Hcorrection.ParseIndices();

	// ---------- Parse the Force indices into an array ---------------
	if (!Force.ParseIndices(argv[1])) return 1;
	if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
	if (!ForceFragments.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;
		if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
		if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
		if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
		if(!ShieldingPASFragments.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.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
	if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
	if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
	if (periode != NULL) { // also look for PAS values
		if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
		if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 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;
	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();
	for(int k=0;k<Time.ColumnCounter;k++)
		Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];

	// +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++

	cout << "Analyzing ..." << endl;

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

	// +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
	// +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
	if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
	if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
	for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
		for(int k=Time.ColumnCounter;k--;) {
			Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
		}
	counter = 0;
	output << "#Order\tFrag.No.\t" << Time.Header << endl;
	output2 << "#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];
				}
		counter += KeySet.FragmentsPerOrder[BondOrder];
		output << BondOrder+1 << "\t" << counter;
		output2 << BondOrder+1 << "\t" << counter;
		for(int k=0;k<Time.ColumnCounter;k++) {
			output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
			if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
				output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
			else
				output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
		}
		output << endl;
		output2 << endl;
	}
	output.close();
	output2.close();

	// +++++++++++++++++++++++++++++++++++++++ Plotting shieldings

	if (periode != NULL) { // also look for PAS values
		if (!CreateDataDeltaForcesOrderPerAtom(ShieldingPAS, ShieldingPASFragments, KeySet, argv[3], "DeltaShieldingsPAS-Order", "Plot of error between approximated shieldings and full shieldings versus the Bond Order", datum)) return 1;
		if (!CreateDataForcesOrderPerAtom(ShieldingPASFragments, KeySet, argv[3], "ShieldingsPAS-Order", "Plot of approximated shieldings versus the Bond Order", datum)) return 1;
		if (!AppendOutputFile(output, argv[3], "ShieldingsPAS-Order.dat" )) return false;
		output << endl << "# Full" << endl;
		for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
			output << j << "\t";
			for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
				output << scientific <<	ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
			output << endl;
		}
	}
	output.close();


	// +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
	if (!CreateDataDeltaEnergyOrder(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 (!CreateDataEnergyOrder(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 (!CreateDataDeltaForcesOrderPerAtom(Force, ForceFragments, KeySet, argv[3], "DeltaForces-Order", "Plot of error between approximated forces and full forces versus the Bond Order", datum)) return 1;

	// min force
	if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMinForces-Order", "Plot of min error between approximated forces and full forces versus the Bond Order", datum, CreateMinimumForce)) return 1;

	// mean force
	if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMeanForces-Order", "Plot of mean error between approximated forces and full forces versus the Bond Order", datum, CreateMeanForce)) return 1;

	// max force
	if (!CreateDataDeltaForcesOrder(Force, ForceFragments, KeySet, argv[3], "DeltaMaxForces-Order", "Plot of max error between approximated forces and full forces versus the Bond Order", datum, CreateMaximumForce)) return 1;

	// ++++++++++++++++++++++++++++++++++++++Plotting Forces vs. Order
	if (!CreateDataForcesOrderPerAtom(ForceFragments, KeySet, argv[3], "Forces-Order", "Plot of approximated forces versus the Bond Order", datum)) return 1;
	if (!AppendOutputFile(output, argv[3], "Forces-Order.dat" )) return false;
	output << endl << "# Full" << endl;
	for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
		output << j << "\t";
		for(int k=0;k<Force.ColumnCounter;k++)
			output << scientific <<	Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
		output << endl;
	}
	output.close();
	// min force
	if (!CreateDataForcesOrder(ForceFragments, KeySet, argv[3], "MinForces-Order", "Plot of min approximated forces versus the Bond Order", datum, CreateMinimumForce)) return 1;

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

	// max force
	if (!CreateDataForcesOrder(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(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 << "]";
	yrange.str("[1e-8:1e+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(), yrange.str().c_str(), "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
	yrange.str("[1e-8:1e+0]");
	// 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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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
	yrange.str("");
	yrange << "[" << EnergyFragments.FindMinValue() << ":" << EnergyFragments.FindMaxValue() << "]";
	if (!CreatePlotOrder(EnergyFragments, KeySet, argv[3], "Energies-Fragment", 5, "below", "y", "", 1, 5, "fragment number", "Energies of each fragment [Ht]", Fragmentxrange.str().c_str(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;

	// +++++++++++++++++++++++++++++++=Ploting min/mean/max forces for each fragment
	yrange.str("");
	yrange << "[" << ForceFragments.FindMinValue() << ":" << ForceFragments.FindMaxValue()<< "]";
	// 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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "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(), yrange.str().c_str(), "1" , "with boxes fillcolor", BoxesForcePlotLine)) return 1;

	// +++++++++++++++++++++++++++++++=Ploting approximated and true shielding for each atom
	if (periode != NULL) { // also look for PAS values
		if(!OpenOutputFile(output, argv[3], "ShieldingsPAS-Order.pyx")) return 1;
		if(!OpenOutputFile(output2, argv[3], "DeltaShieldingsPAS-Order.pyx")) return 1;
		CreatePlotHeader(output, "ShieldingsPAS-Order", 1, "top right", NULL, NULL,	1, 5, "nuclei index", "iso chemical shielding value [ppm]");
		CreatePlotHeader(output2, "DeltaShieldingsPAS-Order", 1, "top right", NULL, NULL,	1, 5, "nuclei index", "iso chemical shielding value [ppm]");
		double step=0.8/KeySet.Order;
		output << "set boxwidth " << step << endl;
		output << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
		output2 << "plot [0:" << ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter]+10 << "]\\" << endl;
		for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
			output << "'ShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
			output2 << "'DeltaShieldingsPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
			if (BondOrder-1 != KeySet.Order)
				output2 << ", \\" << endl;
		}
		output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
		output.close();
		output2.close();
	}

	// 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);
	Free((void **)&dir, "main: *dir");
	cout << "done." << endl;
	return 0;
};

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