//////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <string>
#include <cmath>
#include <vector>
#include<algorithm>
#include <ctime>

using namespace std;

/////////////////////////////////////////////////////

typedef double DP;
typedef unsigned long int ULI;

template <class T>
class NRVec {
private:
	int nn;
	T *v;
public:
	NRVec();
	explicit NRVec(int n);
	NRVec(const T &a, int n);
	NRVec(const T *a, int n);
	NRVec(const NRVec &rhs);
	NRVec & operator=(const NRVec &rhs);
	NRVec & operator=(const T &a);
	inline T & operator[](const int i);
	inline const T & operator[](const int i) const;
	inline int size() const;
	~NRVec();
};

template <class T>
NRVec<T>::NRVec(int n) : nn(n), v(new T[n]) {}

template <class T>
inline T & NRVec<T>::operator[](const int i)	//subscripting
{
	return v[i];
}

template <class T>
inline const T & NRVec<T>::operator[](const int i) const	//subscripting
{
	return v[i];
}

template <class T>
inline int NRVec<T>::size() const
{
	return nn;
}

template <class T>
NRVec<T>::~NRVec()
{
	if (v != 0)
		delete[] (v);
}

typedef const NRVec<DP> Vec_I_DP;
typedef NRVec<DP> Vec_DP, Vec_O_DP, Vec_IO_DP;

inline void nrerror(const string error_text)
	{
		cerr << "Run-time error..." << endl;
		cerr << error_text << endl;
		cerr << "...now exiting to system..." << endl;
		exit(1);
	}

void polint(Vec_I_DP &xa, Vec_I_DP &ya, const DP x, DP &y, DP &dy)
{
	int i,m,ns=0;
	DP den,dif,dift,ho,hp,w;

	int n=xa.size();
	Vec_DP c(n),d(n);
	dif=fabs(x-xa[0]);
	for (i=0;i<n;i++) {
		if ((dift=fabs(x-xa[i])) < dif) {
			ns=i;
			dif=dift;
		}
		c[i]=ya[i];
		d[i]=ya[i];
	}
	y=ya[ns--];
	for (m=1;m<n;m++) {
		for (i=0;i<n-m;i++) {
			ho=xa[i]-x;
			hp=xa[i+m]-x;
			w=c[i+1]-d[i];
			if ((den=ho-hp) == 0.0) nrerror("Error in routine polint");
			den=w/den;
			d[i]=hp*den;
			c[i]=ho*den;
		}
		y += (dy=(2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
	}
}

///////////////////////////////////////////////////// INITIALISATION AND DATA

class State
{
public:

	int					Index;
	double				Energy;
	double				Gamma_a;
	double				Qn;
	string				Config;
	string				Term;
	float				TwoJ;
	bool				Parity;
	bool				Exp;
	char				Status;
	double				Gamma_r;
	double				Gamma_rT2;
	bool				useFlag;

	vector <double>		Nj;
	vector <double>		bu;
	double decRoutSahaCap;

private:
};

/////////////////////////////////////////////////////

class Container
{
public:
	double Value(int Column, int Row) {return Table[Column-1][Row-1];}
	double condVal(float tj1, float tj2, bool p1, bool p2, int col, int row)
	{
		if(tj1==TwoJ1 && tj2==TwoJ2 && p1==P1 && p2==P2)
			return Table[col-1][row-1];
		else if(tj1==TwoJ2 && tj2==TwoJ1 && p1==P2 && p2==P1)
			return Table[row-1][col-1];
		else
			return 0;
	}

	int							numOfColumns;
	float						TwoJ1;
	bool						P1;
	int							numOfRows;
	float						TwoJ2;
	bool						P2;

	vector < vector <double> >	Table;

private:
};

/////////////////////////////////////////////////////

class Symmetry
{
public:
	State *findState(float tj, bool p, int idx)	{if(tj==TwoJ && p==P) return States[idx-1];	}

	float				TwoJ;
	bool				P;
	int					numOfStates;

	vector <State*>		States;

private:
};

/////////////////////////////////////////////////////

class AstroData
{
public:

	string dataSource;
	double lam_obs;
	float F_lam;
	float I_lam;
	string Ion;
	double lam_lab;
	string Mult;
	string LowerTerm;
	string UpperTerm;
	float g1;
	float g2;
	int index;

private:
};

/////////////////////////////////////////////////////

//physical constants
const double Boltzmann			= 1.3806504E-23;
const double Planck				= 6.62606896E-34;
const double h_bar				= 1.05457266E-34;
const double speedOfLight		= 299792458;
const double electMass			= 9.10938215E-31;
const double sinSoLimitcm		= 196664.7;
const double alpha				= 0.0072973525376;
const double tau_o				= 2.418884326505E-17;

//conversion factors
const double invcmtoinvm		= 100;
const double mtoAng				= 1.00E+10;
const double Rydtoinvcm			= 109737.31568525;
const double invcmtoRyd			= 1/Rydtoinvcm;
const double RydtoJ				= 2.1798741E-18;
const double PI					= acos(-1.0);
double alphaValuesCovFac		= 1.0E-6;
double NormWLCovFac				= 1.0E-9;

//keywords
string keywords [] = {"Temperature(K)", "Ni(m^-3)", "Ne(m^-3)", "ResidualCharge", "AstronomicalData",
						"TestsFlag", "RecCoefficient", "NormalisationChoice", "TheoAstroLS", "DecayRoutes"};
bool   keywordsFlags [] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1};

//parameters (mainly associated with keywords)
vector <double> TEMP;
double    Ni = 1.0E10;
double    Ne = 1.0E10;
double	  residualCharge = 2;
double	  zScale;
double    gIon = 1.0;
int		  NormalisationChoice = 0;
int		  NormalisationEntry = 1;
vector	  <double> NormalisationValues;
vector	  <double> alphaValues;
double	  NormalisationWL;
bool	  RecCoef = 0;
bool      astroDataFlag = 0;
string    astroDataFileName;
AstroData *astroNormLine;

bool TheoAstroLSFlag = 0;
bool DecayRoutesFlag = 0;
string TheoAstroLSFileName;
string DecayRoutesConfig;
string DecayRoutesTerm;
int    DecayRoutes2J;
ULI maxRoutes = 10000;
bool TheoAstroLSEntFlag = 0;
int polationPointsNum;

//containers parameters
int       resTwoJ = -999;
bool      resPar  = 1;
int       testFlag = 0;

/////////////////////////////////////////////////////

char buffer[1500];

int numOfOddRes(0), numOfEvenRes(0), numOfAllRes(0);
int numOfOddFSym(0), numOfEvenFSym(0), numOfAllFSym(0);
vector <State*> oddResonances, evenResonances, Resonances;
vector <Container*> ResfValues;
vector <Container*> ResPhotonEnergy;
vector <Container*> FGamma_rul;
Container *FFfValues;
Container *FFGamma_rul = new Container;

int numOfOddBStates(0), numOfEvenBStates(0), numOfAllBStates(0);
int numOfOddBSym(0), numOfEvenBSym(0), numOfAllBSym(0);
int numOfBSym(0);
vector < Symmetry* > allBSym;
vector <Container*> BFVALUES;
vector <Container*> BGamma_rul;
vector < pair <double, State*> >					allFStates;
vector < pair <double, State*> >					allBStates;

int numOfFFTrans = 0;
int numOfFBTrans = 0;
int numOfBBTrans = 0;
vector < pair <double, pair<State*, State*> > >		Transitions;
vector < vector <double> >							Emissivities;
vector < pair < double, AstroData*> >				astroDataVector;

vector <double> dys;

vector <pair <int, int> > theoAstroLSIdxMap;
vector < pair <double, vector <double> > > theoAstroLSValMap;
vector <double> theoAstroLSs;
vector <double> alphaChiSq;
vector <double> ChiSq;
bool Chi2ConfFlag = 0;
int degreesOfFreedom;
double DeltaChi2 = 1.0;
vector < pair <double, vector <double> > > ChiSqComps;
vector < vector <double> > ChiSqCompsSome;

int decRoutCount = 0;

ofstream MapLines;

///////////////////////////////////////////////////// GENERAL-PURPOSE FUNCTIONS

double getContVal(bool Par1, bool Par2, int index1, int index2,
				  vector <State*> oddstates, vector <State*> evenstates, Container* gamma_rul)
{
	int c(-99), r(-99);

	if(Par1 == 1 && Par2 == 0)
	{
		for(int u=0; u<oddstates.size(); u++)
			if(oddstates[u]->Index == index1)
			{c = u; break;}

		for(int v=0; v<evenstates.size(); v++)
			if(evenstates[v]->Index == index2)
			{r = v; break;}
	}

	else if(Par1 == 0 && Par2 == 1)
	{
		for(int u=0; u<oddstates.size(); u++)
			if(oddstates[u]->Index == index2)
			{c = u; break;}

		for(int v=0; v<evenstates.size(); v++)
			if(evenstates[v]->Index == index1)
			{r = v; break;}
	}

	else
	{
		cout <<"Error" << endl;
		exit(-1);
	}

	return gamma_rul->Table[c][r];
}

/////////////////////////////////////////////////////

double getArrOfContVal(vector <Container*> arrOfCont, float tj1, float tj2, bool P1, bool P2, int row, int col)
{
	for(int a=0; a<arrOfCont.size(); a++)
	{
		if(arrOfCont[a]->TwoJ1 == tj1 && arrOfCont[a]->P1 == P1 && arrOfCont[a]->TwoJ2 == tj2 && arrOfCont[a]->P2 == P2 )
			return arrOfCont[a]->Value(row, col);

		if(arrOfCont[a]->TwoJ1 == tj2 && arrOfCont[a]->P1 == P2 && arrOfCont[a]->TwoJ2 == tj1 && arrOfCont[a]->P2 == P1 )
			return arrOfCont[a]->Value(col, row);
	}

	cout << "Unrecognised Value" << endl;
	exit(-1);
}

/////////////////////////////////////////////////////

double Gamma_rjl(double f, double nuRyd, double gl, double gu)
{
	return (alpha * alpha * alpha * nuRyd * nuRyd * gl * f) / (2 * gu * tau_o);
}

////////////////////////////////////////////////////

double Saha(double Ne, double Ni, double gu, double gi, double delE, double T)
{
	return  Ne * Ni * (gu/(2*gi)) * exp(-(delE/(Boltzmann*T))) * pow((Planck*Planck)/(2*PI*electMass*Boltzmann*T), 1.5);
}

/////////////////////////////////////////////////////

bool E1Possible(float tj1, float tj2, bool p1, bool p2)
{
	float d = fabs(tj1 - tj2);

	if((d == 2 || d == 0) && (p1 + p2) == 1)
		return true;
	else
		return false;
}

/////////////////////////////////////////////////////

double WnToWlVacAng(double waveNum) //converts from wavenumber (1/cm) in vacuum to wavelength in vacuum in Angstrom
{
	return 1.0E8 / waveNum;
}

/////////////////////////////////////////////////////

double WnToWlAirAng(double waveNum) //converts from wavenumber (1/cm) in vacuum to wavelength in air in Angstrom
{
	return (1.0E8)/(waveNum*(1.0 + 6432.8E-8 + 2949810.0/(146.0E8 - waveNum*waveNum)  + 25540./(41.0E8 - waveNum*waveNum)));
}

/////////////////////////////////////////////////////

double recCombToEmissivity(double NE, double NI, double ALPHA, double LAMBDA)		//converts effective recombination coefficients to emissivities
{
	return (NE * NI * ALPHA * Planck * speedOfLight) / LAMBDA;
}

/////////////////////////////////////////////////////

double emissivityToRecCoef(double NE, double NI, double EPSILON, double LAMBDA)		//converts emissivities to effective recombination coefficients
{
	return (EPSILON * LAMBDA) / (NE * NI * Planck * speedOfLight);
}


///////////////////////////////////////////////////// READING INPUT DATA FILES

void readInputFile()
{
	ifstream Input;
	Input.open("Input");
	if(Input.fail())
	{
		cout << "No Input file" << endl;
		exit(-1);
	}

	string keyword;

	TEMP.resize(1);
	TEMP[0] = 10000;

	int numOfT;

	while(!Input.eof())
	{
		Input >> keyword;

		if(keyword == keywords[0] && keywordsFlags[0])
		{
			bool tempFlag;

			Input >> tempFlag >> numOfT;
			TEMP.resize(numOfT);
			if(!tempFlag)
			{
				for(int z=0; z<numOfT; z++)
					Input >> TEMP[z];
			}
			else
			{
				double initTemp, tempInterval;
				Input >> initTemp >> tempInterval;
				for(int y=0; y<numOfT; ++y)
					TEMP[y] = initTemp + (y * tempInterval);
			}

			keywordsFlags[0] = 0;
		}

		else if(keyword == keywords[1] && keywordsFlags[1])
		{
			Input >> Ni;
			keywordsFlags[1] = 0;
		}

		else if(keyword == keywords[2] && keywordsFlags[2])
		{
			Input >> Ne;
			keywordsFlags[2] = 0;
		}

		else if(keyword == keywords[3] && keywordsFlags[3])
		{
			Input >> residualCharge;
			keywordsFlags[3] = 0;
		}

		else if(keyword == keywords[4] && keywordsFlags[4])
		{
			Input >> astroDataFlag;
			if(astroDataFlag)
				Input >> astroDataFileName;

			keywordsFlags[4] = 0;
		}

		else if(keyword == keywords[5] && keywordsFlags[5])
		{
			Input >> testFlag;
			keywordsFlags[5] = 0;
		}

		else if(keyword == keywords[6] && keywordsFlags[6])
		{
			Input >> RecCoef;

			if(RecCoef != 1)
				RecCoef = 0;
			keywordsFlags[6] = 0;
		}

		else if(keyword == keywords[7] && keywordsFlags[7])
		{
			Input >> NormalisationChoice;
			keywordsFlags[7] = 0;

			if(NormalisationChoice == 0)
				Input.getline(buffer, 1500, '\n');

			else if(NormalisationChoice == 1)
			{
				Input >> NormalisationEntry;
				Input.getline(buffer, 1500, '\n');
			}

			else if(NormalisationChoice == 2)
			{
				Input.getline(buffer, 1500, '\n');

				double value;
				for(int a=0; a<TEMP.size(); a++)
				{
					Input >> value;
					NormalisationValues.push_back(value);
				}
				Input.getline(buffer, 1500, '\n');
			}

			else if(NormalisationChoice == 3)
			{
				Input.getline(buffer, 1500, '\n');

				double alphaValue;

				for(int a=0; a<TEMP.size(); a++)
				{
					Input >> alphaValue;
					alphaValue *= alphaValuesCovFac;
					alphaValues.push_back(alphaValue);
				}

				Input >> NormalisationWL;
				NormalisationWL *= NormWLCovFac;

				Input.getline(buffer, 1500, '\n');

				for(int b=0; b<TEMP.size(); b++)
				{
					alphaValue = alphaValues[b];
					double emiss = recCombToEmissivity(Ne, Ni, alphaValue, NormalisationWL);
					NormalisationValues.push_back(emiss);
				}
			}

			else if(NormalisationChoice == 4)
			{
				Input >> polationPointsNum;
				Input.getline(buffer, 1500, '\n');

				int numOfPoints;
				Input >> numOfPoints;

				if(numOfPoints < polationPointsNum)
				{
					cout << "Number of points is not suitable for the required interpolation scheme" << endl;
					polationPointsNum = numOfPoints < 3 ? numOfPoints : 3;
				}

				DP dy, temperature, alphaValue;
				Vec_DP xa(numOfPoints), xxa(polationPointsNum);
				Vec_DP ya(numOfPoints), yya(polationPointsNum);

				for(int a=0; a<numOfPoints; a++)
					Input >> xa[a] >> ya[a];

				Input >> NormalisationWL;
				NormalisationWL *= NormWLCovFac;

				/////////////////////////////////////////

				for(int b=0; b<TEMP.size(); b++)
				{
					temperature = TEMP[b];
					int index;

					if(temperature < xa[0])
						index = 0;
					else if(temperature > xa[xa.size()-1])
						index = xa.size() - polationPointsNum;
					else
					{
						for(int c=0; c<xa.size(); c++)
						{
							if(xa[c] == temperature)
							{
								alphaValue = ya[c];
								dy = 0;
								index = -999;
								break;
							}

							else if(xa[c] > temperature)
							{
								index = (int)(c - ceil(polationPointsNum/2.0));

								while(index < 0)
									index++;
								while(index > (xa.size() - polationPointsNum))
									index--;

								break;
							}
						}
					}

					if(index > -1)
					{
						for(int d=0; d<polationPointsNum; d++)
						{
							xxa[d] = xa[index+d];
							yya[d] = ya[index+d];
						}
						polint(xxa,yya,temperature,alphaValue,dy);
					}

					dys.push_back(dy);

					alphaValue *= alphaValuesCovFac;
					alphaValues.push_back(alphaValue);

					double emiss = recCombToEmissivity(Ne, Ni, alphaValue, NormalisationWL);
					NormalisationValues.push_back(emiss);
				}
			}

			else
				NormalisationChoice = 0;
		}

		else if(keyword == keywords[8] && keywordsFlags[8])
		{
			Input >> TheoAstroLSFlag;
			if(TheoAstroLSFlag)
			{
				Input >> TheoAstroLSFileName >> TheoAstroLSEntFlag >> Chi2ConfFlag;
				if(Chi2ConfFlag)
					Input >> degreesOfFreedom >> DeltaChi2;
			}

			keywordsFlags[8] = 0;
		}

		else if(keyword == keywords[9] && keywordsFlags[9])
		{
			Input >> DecayRoutesFlag;
			if(DecayRoutesFlag)
				Input >> DecayRoutesConfig >> DecayRoutesTerm >> DecayRoutes2J >> maxRoutes;

			keywordsFlags[9] = 0;
		}

		else
			Input.getline(buffer, 1500, '\n');

		/////////////////////////////

		if(Input.eof())
		{
			for(int z=0; z<sizeof(keywordsFlags); z++)
				if(keywordsFlags[z])
					cout << keywords[z] << " is missing. Default value will be used. See output file." << endl << endl;
		}

	}

	zScale = residualCharge * residualCharge;

	Input.close();
}

/////////////////////////////////////////////////////

void readResEmisFile()
{
	ifstream ResEmis;
	ResEmis.open("ResEmis");
	if(ResEmis.fail())
	{
		cout << "No ResEmis file" << endl;
		exit(-1);
	}

	ResEmis.getline(buffer, 1500, '\n');
	ResEmis >> numOfAllRes;

	Resonances.resize(numOfAllRes);
	double factor = zScale * RydtoJ/h_bar;

	for(int i=0; i<numOfAllRes; ++i)
	{
		Resonances[i] = new State;

		ResEmis >> Resonances[i]->Index
				>> Resonances[i]->Energy
				>> Resonances[i]->Gamma_a		//This is the width. It will be converted to Gamma_a
				>> Resonances[i]->Config
				>> Resonances[i]->Term
				>> Resonances[i]->TwoJ
				>> Resonances[i]->Parity
				>> Resonances[i]->Exp
				>> Resonances[i]->useFlag;		//remove

		Resonances[i]->Gamma_a = factor * Resonances[i]->Gamma_a;
		Resonances[i]->Status  = 'F';


		if(Resonances[i]->Parity == 1)
		{
			oddResonances.push_back(Resonances[i]);
			numOfOddRes++;
		}
		else
		{
			evenResonances.push_back(Resonances[i]);
			numOfEvenRes++;
		}
	}

// Read f-value data

	string comment = "abcd";

	while(comment != "EndOfComment")
		ResEmis >> comment;

	float twoJB = 0.0;
	bool parityB = 0;
	int numOfB = 0;

	ResEmis >> numOfBSym;
	ResEmis.getline(buffer, 1500, '\n');

	ResfValues.resize(numOfBSym);

	for(int k=0; k<numOfBSym; k++)
	{
		ResEmis >> twoJB >> parityB >> numOfB;
		ResEmis.getline(buffer, 1500, '\n');

		ResfValues[k] = new Container;

		ResfValues[k]->numOfColumns		= numOfAllRes;
		ResfValues[k]->numOfRows		= numOfB;


		ResfValues[k]->TwoJ1	= resTwoJ;
		ResfValues[k]->P1		= resPar;

		ResfValues[k]->TwoJ2	= twoJB;
		ResfValues[k]->P2		= parityB;

		ResfValues[k]->Table.resize(numOfAllRes);

		for(int z=0; z<ResfValues[k]->Table.size(); z++)
			ResfValues[k]->Table[z].resize(numOfB);

		for(int f=0; f<ResfValues[k]->numOfRows; f++)
			for(int g=0; g<ResfValues[k]->numOfColumns; g++)
				ResEmis >> ResfValues[k]->Table[g][f];

	} //end of: for(int k=0; k<numOfBSym; k++)


// Read photon energies data

	comment = "abcd";

	while(comment != "EndOfComment")
		ResEmis >> comment;

	twoJB = 0.0;
	parityB = 0;
	numOfB = 0;

	ResPhotonEnergy.resize(numOfBSym);

	for(int l=0; l<numOfBSym; l++)
	{
		ResEmis >> twoJB >> parityB >> numOfB;
		ResEmis.getline(buffer, 1500, '\n');

		ResPhotonEnergy[l] = new Container;

		ResPhotonEnergy[l]->numOfColumns	= numOfAllRes;
		ResPhotonEnergy[l]->numOfRows		= numOfB;

		ResPhotonEnergy[l]->TwoJ1	= resTwoJ;
		ResPhotonEnergy[l]->P1	= resPar;

		ResPhotonEnergy[l]->TwoJ2	= twoJB;
		ResPhotonEnergy[l]->P2		= parityB;

		ResPhotonEnergy[l]->Table.resize(numOfAllRes);

		for(int z=0; z<ResPhotonEnergy[l]->Table.size(); z++)
			ResPhotonEnergy[l]->Table[z].resize(numOfB);

		for(int f=0; f<ResfValues[l]->numOfRows; f++)
			for(int g=0; g<ResPhotonEnergy[l]->numOfColumns; g++)
				ResEmis >> ResPhotonEnergy[l]->Table[g][f];
	}  //end of: for(int l=0; l<numOfBSym; l++)

// Read FF f-values

	if(numOfOddRes > 0 && numOfEvenRes > 0)
	{
		comment = "abcd";

		while(comment != "EndOfComment")
			ResEmis >> comment;

		FFfValues = new Container;

		FFfValues->numOfColumns = numOfOddRes;
		FFfValues->numOfRows    = numOfEvenRes;

		FFfValues->Table.resize(numOfOddRes);

		for(int h=0; h<FFfValues->Table.size(); h++)
			FFfValues->Table[h].resize(numOfEvenRes);

		for(int i=0; i<FFfValues->numOfRows; i++)
			for(int j=0; j<FFfValues->numOfColumns; j++)
				ResEmis >> FFfValues->Table[j][i];
	}

	ResEmis.close();
}

/////////////////////////////////////////////////////

void readELEVEmisFile()
{

	ifstream ELEVEmis;
	ELEVEmis.open("ELEVEmis");
	if(ELEVEmis.fail())
	{
		cout << "No ELEVEmis file" << endl;
		exit(-1);
	}

	string comment = "abcd";

	while(comment != "EndOfComment")
		ELEVEmis >> comment;

	ELEVEmis >> numOfAllBSym;

	allBSym.resize(numOfAllBSym);

	for(int j=0; j<numOfAllBSym; ++j)
	{
		allBSym[j] = new Symmetry;

		ELEVEmis >> allBSym[j]->TwoJ
				 >> allBSym[j]->P
				 >> allBSym[j]->numOfStates;

		allBSym[j]->States.resize(allBSym[j]->numOfStates);

		for(int f=0; f<allBSym[j]->numOfStates; f++)
		{
			allBSym[j]->States[f] = new State;

			ELEVEmis >> allBSym[j]->States[f]->Index
					 >> allBSym[j]->States[f]->Energy
					 >> allBSym[j]->States[f]->Qn
					 >> allBSym[j]->States[f]->Config
					 >> allBSym[j]->States[f]->Term
					 >> allBSym[j]->States[f]->Exp;

			allBSym[j]->States[f]->TwoJ   = allBSym[j]->TwoJ;
			allBSym[j]->States[f]->Parity = allBSym[j]->P;
			allBSym[j]->States[f]->Status = 'B';
		}
	}

	ELEVEmis.close();
}

////////////////////////////////////////

void readFVALUEFile()
{
	ifstream FVALUE;
	FVALUE.open("FVALUE");
	if(FVALUE.fail())
	{
		cout << "No FVALUE file" << endl;
		exit(-1);
	}

	FVALUE.getline(buffer, 1500, '\n');

	int twoJ1, twoJ2, numOfTrans, temp, numOfCols, numOfRows, Idx1, Idx2;
	bool P1, P2;
	double fL, fV;

	while(!FVALUE.eof())
	{
		FVALUE >> temp >> twoJ1 >> P1 >> temp >> twoJ2 >> P2 >> numOfTrans;

		if(twoJ1 == 0 && twoJ2 == 0 && numOfTrans == 0)
			break;

		for(int m=0; m<allBSym.size(); m++)
		{
			if(allBSym[m]->TwoJ == twoJ1 && allBSym[m]->P == P1)
				numOfCols = allBSym[m]->numOfStates;

			if(allBSym[m]->TwoJ == twoJ2 && allBSym[m]->P == P2)
				numOfRows = allBSym[m]->numOfStates;
		}

		Container *container = new Container;

		container->TwoJ1 = twoJ1;
		container->TwoJ2 = twoJ2;
		container->P1 = P1;
		container->P2 = P2;

		container->Table.resize(numOfCols);

		for(int l=0; l<numOfCols; l++)
			container->Table[l].resize(numOfRows);

		for(int k=0; k<numOfTrans; k++)
		{
			FVALUE >> Idx1 >> Idx2 >> fL >> fV;
			container->Table[Idx1-1][Idx2-1] = fL;
		}

		BFVALUES.push_back(container);
	}

	FVALUE.close();
}

///////////////////////////////////////////////////// FREE-FREE SECTION

vector < pair <double, State*> > arrangeAllFStates(vector <State*> fstates, vector < pair <double, State*> > states)
{
	for(int f=0; f<fstates.size(); f++)
	{
		pair <double, State*> st;
		st.first  = fabs(fstates[f]->Energy);
		st.second = fstates[f];

		states.push_back(st);
	}

	for(int s=0; s<states.size(); s++)
		states[s].first *= -1.0;

	stable_sort(states.begin(), states.end());

	for(int m=0; m<states.size(); m++)
		states[m].first *= -1.0;

	return states;
}

//////////////////////////////////////////////////////

Container* findFFGamma_rul(Container* fffValues, vector <State*> fstates, Container* ffgamma_rul)
{
	ffgamma_rul->Table.resize(fffValues->Table.size());

	for(int l=0; l<ffgamma_rul->Table.size(); l++)
		ffgamma_rul->Table[l].resize(fffValues->Table[l].size());


	for(int k=0; k<ffgamma_rul->Table.size(); k++)
	{
		for(int s=0; s<ffgamma_rul->Table[k].size(); s++)
		{
			double f = fabs(fffValues->Table[k][s]);

			State *state1 = new State;
			State *state2 = new State;

			state1 = oddResonances[k];
			state2 = evenResonances[s];

			double nuRyd = zScale * fabs(fabs(state1->Energy) - fabs(state2->Energy));

			double gl = 1.0;
			double gu = state1->Energy > state2->Energy ? state1->TwoJ + 1 : state2->TwoJ + 1;

			ffgamma_rul->Table[k][s] = Gamma_rjl(f, nuRyd, gl, gu);
		}
	}

	return ffgamma_rul;
}

/////////////////////////////////////////////////////

void findFFGamma_r(vector < pair <double, State*> > states, Container* ffgamma_rul)
{
	for(int s=0; s<(states.size()-1); s++)
	{
		states[s].second->Gamma_r = 0.0;

		for(int t=(s+1); t<states.size(); t++)
		{
			float tj1  = states[s].second->TwoJ;
			float tj2  = states[t].second->TwoJ;
			bool  p1   = states[s].second->Parity;
			bool  p2   = states[t].second->Parity;
			int   idx1 = states[s].second->Index;
			int   idx2 = states[t].second->Index;

			bool pos = E1Possible(tj1, tj2, p1, p2);

			if(pos == true)
			{
				double gamrul = getContVal(p1, p2, idx1, idx2, oddResonances, evenResonances, ffgamma_rul);
				states[s].second->Gamma_r += gamrul;
			}
		}
	}
}

///////////////////////////////////////////////////// FREE-BOUND SECTION

vector <Container*> findFGamma_rul(int numOfSym, vector <State*> states, vector <Container*> fValues,
				   vector <Container*> photonEnergy, vector <Container*> gamma_rul)
{
	gamma_rul.resize(numOfSym);

	for(int k=0; k<numOfSym; k++)
	{
		gamma_rul[k] = new Container;

		gamma_rul[k]->numOfColumns	= fValues[k]->numOfColumns;
		gamma_rul[k]->numOfRows		= fValues[k]->numOfRows;

		gamma_rul[k]->TwoJ1		= fValues[k]->TwoJ1;
		gamma_rul[k]->P1		= fValues[k]->P1;

		gamma_rul[k]->TwoJ2		= fValues[k]->TwoJ2;
		gamma_rul[k]->P2		= fValues[k]->P2;

		gamma_rul[k]->Table.resize(gamma_rul[k]->numOfColumns);

		for(int z=0; z<gamma_rul[k]->Table.size(); z++)
			gamma_rul[k]->Table[z].resize(gamma_rul[k]->numOfRows);

		for(int f=0; f<gamma_rul[k]->numOfRows; f++)
			for(int g=0; g<gamma_rul[k]->numOfColumns; g++)
				gamma_rul[k]->Table[g][f] =
				Gamma_rjl(fValues[k]->Table[g][f],
						  photonEnergy[k]->Table[g][f],
						  (fValues[k]->TwoJ2 + 1),
						  (states[g]->TwoJ + 1));

	} // end of: for(int k=0; k<numOfSym; k++)

	return gamma_rul;
}

//////////////////////////////////////////////////////

void findFGamma_r(int numOfStates, int numOfSym, vector <State*> states, vector <Container*> gamma_rul)
{
	if((numOfOddRes * numOfEvenRes) == 0)
	{
		for(int n=0; n<numOfStates; n++)
			states[n]->Gamma_r = 0.0;
	}

	for(int n=0; n<numOfStates; n++)
		states[n]->Gamma_rT2 = 0.0;

	for(int k=0; k<numOfSym; k++)
		for(int f=0; f<gamma_rul[k]->numOfRows; f++)
			for(int g=0; g<gamma_rul[k]->numOfColumns; g++)
			{
				states[g]->Gamma_r   += gamma_rul[k]->Table[g][f];
				states[g]->Gamma_rT2 += gamma_rul[k]->Table[g][f];
			}
}

/////////////////////////////////////////////////////

void findbu(int numOfStates, vector <State*> states)
{
	for(int m=0; m<numOfStates; m++)
	{
		states[m]->bu.resize(TEMP.size());

		for(int h=0; h<TEMP.size(); h++)
			if(states[m]->Gamma_r == 0.0)
				states[m]->bu[h] = 0.0;
			else
				states[m]->bu[h] = states[m]->Gamma_a / (states[m]->Gamma_a + states[m]->Gamma_r);
	}
}

/////////////////////////////////////////////

void findFFNj(vector < pair <double, State*> > ffstates, vector <double> temperature, Container* ffgamma_rul)
{
	for(int s=0; s<ffstates.size(); s++)
	{
		double gu   = ffstates[s].second->TwoJ + 1;
		double gi   = gIon;
		double delE = RydtoJ * zScale * ffstates[s].second->Energy ;

		ffstates[s].second->Nj.resize(temperature.size());

		for(int k=0; k<temperature.size(); k++)
		{
			ffstates[s].second->Nj[k] = ffstates[s].second->bu[k] * Saha(Ne, Ni, gu, gi, delE, temperature[k]);

			if(k == 0)
				ffstates[s].second->decRoutSahaCap = ffstates[s].second->Nj[k];

			for(int l=0; l<s; l++)
			{
				float tj1  = ffstates[l].second->TwoJ;
				float tj2  = ffstates[s].second->TwoJ;
				bool  p1   = ffstates[l].second->Parity;
				bool  p2   = ffstates[s].second->Parity;
				int   idx1 = ffstates[l].second->Index;
				int   idx2 = ffstates[s].second->Index;

				bool pos = E1Possible(tj1, tj2, p1, p2);

				if(pos == true)
				{
					double Nk    = ffstates[l].second->Nj[k];
					double Gamkj = getContVal(p1, p2, idx1, idx2, oddResonances, evenResonances, ffgamma_rul);
					double Gamjr  = ffstates[s].second->Gamma_r;
					double Gamja = ffstates[s].second->Gamma_a;

					double radPop = (Nk * Gamkj) / (Gamjr + Gamja);

					ffstates[s].second->Nj[k] += radPop;
				}
			}
		}
	}
}

///////////////////////////////////////////////////// BOUND-BOUND SECTION

vector <Container*> findBGamma_rul(vector <Container*> fValues, vector < Symmetry* > symmetries, vector <Container*> gamma_rul)
{
	for(int m=0; m<fValues.size(); m++)
	{
		Container *container = new Container;

		container->TwoJ1	= fValues[m]->TwoJ1;
		container->TwoJ2	= fValues[m]->TwoJ2;
		container->P1		= fValues[m]->P1;
		container->P2		= fValues[m]->P2;

		container->Table.resize(fValues[m]->Table.size());

		for(int l=0; l<container->Table.size(); l++)
			container->Table[l].resize(fValues[m]->Table[l].size());

		for(int k=0; k<container->Table.size(); k++)
		{
			for(int s=0; s<container->Table[k].size(); s++)
			{
				double f = fabs(fValues[m]->Table[k][s]);

				State *state1 = new State;
				State *state2 = new State;

				for(int t=0; t<symmetries.size(); t++)
				{
					if(symmetries[t]->TwoJ == container->TwoJ1 && symmetries[t]->P == container->P1)
						state1 = symmetries[t]->States[k];

					if(symmetries[t]->TwoJ == container->TwoJ2 && symmetries[t]->P == container->P2)
						state2 = symmetries[t]->States[s];
				}

				double nuRyd = zScale * fabs(fabs(state1->Energy) - fabs(state2->Energy));

				double gl = 1.0;
				double gu = state1->Energy > state2->Energy ? state1->TwoJ + 1 : state2->TwoJ + 1;

				container->Table[k][s] = Gamma_rjl(f, nuRyd, gl, gu);
			}
		}

		gamma_rul.push_back(container);

	} //end of: for(int m=0; m<fValues.size(); m++)

	return gamma_rul;
}

/////////////////////////////////////////////////////

vector < pair <double, State*> > arrangeAllBStates(vector < Symmetry* > symmetries, vector < pair <double, State*> > states)
{
	for(int f=0; f<symmetries.size(); f++)
	{
		for(int g=0; g<symmetries[f]->States.size(); g++)
		{
			pair <double, State*> st;
			st.first  = fabs(symmetries[f]->States[g]->Energy);
			st.second = symmetries[f]->States[g];

			states.push_back(st);
		}
	}

	stable_sort(states.begin(), states.end());

	for(int s=0; s<states.size(); s++)
		states[s].first *= -1.0;

	return states;
}

/////////////////////////////////////////////////////

void findBGamma_r(vector < pair <double, State*> > states, vector <Container*> gamma_rul)
{
	for(int s=0; s<(states.size()-1); s++)
	{
		states[s].second->Gamma_r = 0.0;

		for(int t=(s+1); t<states.size(); t++)
		{
			float tj1  = states[s].second->TwoJ;
			float tj2  = states[t].second->TwoJ;
			bool  p1   = states[s].second->Parity;
			bool  p2   = states[t].second->Parity;
			int   idx1 = states[s].second->Index;
			int   idx2 = states[t].second->Index;

			bool pos = E1Possible(tj1, tj2, p1, p2);

			if(pos == true)
				states[s].second->Gamma_r += getArrOfContVal(gamma_rul, tj1, tj2, p1, p2, idx1, idx2);
		}
	}
}

/////////////////////////////////////////////

void findBNj(vector <State*> Fstates, vector < pair <double, State*> > Bstates, vector <double> temperature,
			 vector <Container*> gamma_rul)
{
	for(int s=0; s<Bstates.size(); s++)
	{
		Bstates[s].second->Nj.resize(temperature.size());

		for(int k=0; k<temperature.size(); k++)
		{
			Bstates[s].second->Nj[k] = 0.0;

			for(int n=0; n<Fstates.size(); n++)
			{
				float tj1  = Fstates[n]->TwoJ;
				float tj2  = Bstates[s].second->TwoJ;
				bool  p1   = Fstates[n]->Parity;
				bool  p2   = Bstates[s].second->Parity;
				int   idx1 = Fstates[n]->Index;
				int   idx2 = Bstates[s].second->Index;

				bool pos = E1Possible(tj1, tj2, p1, p2);

				if(pos == true)
				{
					double Nk    = Fstates[n]->Nj[k];
					double Gamkj = getArrOfContVal(FGamma_rul, resTwoJ, tj2, resPar, p2, idx1, idx2);
					double Gamjr = Bstates[s].second->Gamma_r;

					Bstates[s].second->Nj[k] += (Nk * Gamkj) / Gamjr;
				}
			}

			for(int l=0; l<s; l++)
			{
				float tj1  = Bstates[l].second->TwoJ;
				float tj2  = Bstates[s].second->TwoJ;
				bool  p1   = Bstates[l].second->Parity;
				bool  p2   = Bstates[s].second->Parity;
				int   idx1 = Bstates[l].second->Index;
				int   idx2 = Bstates[s].second->Index;

				bool pos = E1Possible(tj1, tj2, p1, p2);

				if(pos == true)
				{
					double Nk    = Bstates[l].second->Nj[k];
					double Gamkj = getArrOfContVal(gamma_rul, tj1, tj2, p1, p2, idx1, idx2);
					double Gamjr = Bstates[s].second->Gamma_r;

					Bstates[s].second->Nj[k] += (Nk * Gamkj) / Gamjr;
				}
			}
		}
	}
}
///////////////////////////////////////////////////// MIXED SECTION

vector < pair <double, pair<State*, State*> > > findTransitions(vector <State*> fstates, vector < pair <double, State*> > allbstates,
					 vector < pair <double, pair<State*, State*> > > transits, vector < pair <double, State*> > allfstates)
{
	for(int u=0; u<(allfstates.size()-1); u++)
	{
		for(int v=(u+1); v<allfstates.size(); v++)
		{
			float tj1  = allfstates[u].second->TwoJ;
			float tj2  = allfstates[v].second->TwoJ;
			bool  p1   = allfstates[u].second->Parity;
			bool  p2   = allfstates[v].second->Parity;

			bool pos = E1Possible(tj1, tj2, p1, p2);

			if(pos == true)
			{
				pair <double, pair<State*, State*> > trans;

				trans.first         = zScale * fabs(fabs(allfstates[u].second->Energy) - fabs(allfstates[v].second->Energy));
				trans.second.first  = allfstates[u].second;
				trans.second.second = allfstates[v].second;

				transits.push_back(trans);
				numOfFFTrans++;
			}
		}
	}

	for(int b=0; b<fstates.size(); b++)
	{
		for(int c=0; c<allbstates.size(); c++)
		{
			float tj1  = fstates[b]->TwoJ;
			float tj2  = allbstates[c].second->TwoJ;
			bool  p1   = fstates[b]->Parity;
			bool  p2   = allbstates[c].second->Parity;
			int   idx1 = fstates[b]->Index;
			int   idx2 = allbstates[c].second->Index;

			bool pos = E1Possible(tj1, tj2, p1, p2);

			if(pos == true)
			{
				pair <double, pair<State*, State*> > trans;

				trans.first         = zScale * (fstates[b]->Energy + fabs(allbstates[c].second->Energy));
				trans.second.first  = fstates[b];
				trans.second.second = allbstates[c].second;

				transits.push_back(trans);
				numOfFBTrans++;
			}
		}
	}

	for(int s=0; s<(allbstates.size()-1); s++)
	{
		for(int t=(s+1); t<allbstates.size(); t++)
		{
			float tj1  = allbstates[s].second->TwoJ;
			float tj2  = allbstates[t].second->TwoJ;
			bool  p1   = allbstates[s].second->Parity;
			bool  p2   = allbstates[t].second->Parity;

			bool pos = E1Possible(tj1, tj2, p1, p2);

			if(pos == true)
			{
				pair <double, pair<State*, State*> > trans;

				trans.first         = zScale * fabs(fabs(allbstates[s].second->Energy) - fabs(allbstates[t].second->Energy));
				trans.second.first  = allbstates[s].second;
				trans.second.second = allbstates[t].second;

				transits.push_back(trans);
				numOfBBTrans++;
			}
		}
	}

	stable_sort(transits.begin(), transits.end());

	return transits;
}

/////////////////////////////////////////////////////

vector < vector <double> > findEmissivities(vector < pair <double, pair<State*, State*> > > transitions,
											Container* ffgamma_rul, vector <Container*> fgamma_rul,
											vector <Container*> bgamma_rul, vector < vector <double> > emis)
{
	for(int j=0; j<transitions.size(); j++)
	{
		vector <double> EmissivityVec;

		for(int m=0; m<TEMP.size(); m++)		//here
		{
			float tj1  = transitions[j].second.first->TwoJ;
			float tj2  = transitions[j].second.second->TwoJ;
			bool  p1   = transitions[j].second.first->Parity;
			bool  p2   = transitions[j].second.second->Parity;
			int   idx1 = transitions[j].second.first->Index;
			int   idx2 = transitions[j].second.second->Index;

			double Nk    = transitions[j].second.first->Nj[m];
			double E     = RydtoJ * zScale * (transitions[j].second.first->Energy - transitions[j].second.second->Energy);
			double Gamkj;

			if(transitions[j].second.first->Status == 'F' && transitions[j].second.second->Status == 'F')
				Gamkj = getContVal(p1, p2, idx1, idx2, oddResonances, evenResonances, ffgamma_rul);
			else if(transitions[j].second.first->Status == 'F' && transitions[j].second.second->Status == 'B')
				Gamkj = getArrOfContVal(fgamma_rul, resTwoJ, tj2, resPar, p2, idx1, idx2);
			else if(transitions[j].second.first->Status == 'B' && transitions[j].second.second->Status == 'B')
				Gamkj = getArrOfContVal(bgamma_rul, tj1, tj2, p1, p2, idx1, idx2);
			else
			{
				cout <<"Error!" << endl;
				exit(-1);
			}

			double Emissivity = Nk * Gamkj * E;

			EmissivityVec.push_back(Emissivity);
		}

		emis.push_back(EmissivityVec);
	}

	return emis;
}

/////////////////////////////////////////////////////

void readAstroData()
{
	ifstream astroData(astroDataFileName.c_str());;
	if(astroData.fail())
	{
		cout << "No " << astroDataFileName << " file" << endl;
		exit(-1);
	}

	string datSource, comment;
	astroData >> datSource;

	while(comment != "EndOfComment")
		astroData >> comment;

	astroNormLine = new AstroData;

	astroNormLine->dataSource = datSource;

	astroData >> astroNormLine->lam_obs
			  >> astroNormLine->F_lam
			  >> astroNormLine->I_lam
			  >> astroNormLine->Ion
			  >> astroNormLine->lam_lab
			  >> astroNormLine->Mult
			  >> astroNormLine->LowerTerm
			  >> astroNormLine->UpperTerm
			  >> astroNormLine->g1
			  >> astroNormLine->g2;

	int astroDataCount = 0;

	while(!astroData.eof())
	{
		AstroData *astd = new AstroData;

		astroDataCount++;
		astd->index = astroDataCount;

		astd->dataSource = datSource;

		astroData >> astd->lam_obs
			      >> astd->F_lam
				  >> astd->I_lam
				  >> astd->Ion;

		if(astroData.eof())
		{
			astd = NULL;
			break;
		}

		astroData >> astd->lam_lab
				  >> astd->Mult
				  >> astd->LowerTerm
				  >> astd->UpperTerm
				  >> astd->g1
				  >> astd->g2;

	    pair <double, AstroData*> astDatP;

		astDatP.first  = astd->lam_obs;   //for changing priority
		astDatP.second = astd;
		astroDataVector.push_back(astDatP);
	}

	stable_sort(astroDataVector.begin(), astroDataVector.end());

	astroData.close();
}

/////////////////////////////////////////////////////

void TheoAstroLS()
{
	if(NormalisationChoice <= 0 || !astroDataFlag)
		cout << "Unable to compute least squares as there is no normalisation or observational data." << endl;

	else
	{
		ifstream MapFile(TheoAstroLSFileName.c_str());;
		if(MapFile.fail())
		{
			cout << "No '" << TheoAstroLSFileName << "' file" << endl;
			exit(-1);
		}

		for(;;)
		{
			int Idx1, Idx2;

			MapFile >> Idx1;
			if(Idx1 == -999 || MapFile.eof())
				break;
			else
				MapFile >> Idx2;

			if(Transitions[Transitions.size()-Idx2].second.first->Status == 'B' ||
			   Transitions[Transitions.size()-Idx2].second.first->useFlag )
			{
				pair <int, int> map;
				map.first  = Idx1;
				map.second = Idx2;

				theoAstroLSIdxMap.push_back(map);
			}
		}

		MapFile.close();

		theoAstroLSValMap.resize(theoAstroLSIdxMap.size());
		ChiSqComps.resize(theoAstroLSIdxMap.size());

		for(int z=0; z<theoAstroLSValMap.size(); z++)
		{
			theoAstroLSValMap[z].second.resize(TEMP.size());
			ChiSqComps[z].second.resize(TEMP.size());
		}

		for(int a=0; a<theoAstroLSIdxMap.size(); a++)
		{
			int astroIdx = theoAstroLSIdxMap[a].first;
			int theoIdx  = theoAstroLSIdxMap[a].second;

			double astroNorm = astroDataVector[astroIdx-1].second->I_lam / astroNormLine->I_lam;
			theoAstroLSValMap[a].first = astroNorm;

			if(NormalisationChoice == 1)
			{
				for(int p=0; p<TEMP.size(); p++)
				{
					double nEmiss = Emissivities[Transitions.size() - NormalisationEntry][p];
					double theoNorm = Emissivities[Transitions.size() - theoIdx][p] / nEmiss;

					theoAstroLSValMap[a].second[p] = theoNorm;
				}
			}
			else
			{
				for(int q=0; q<NormalisationValues.size(); q++)
				{
					double nEmiss = NormalisationValues[q];
					double theoNorm = Emissivities[Transitions.size() - theoIdx][q] / nEmiss;

					theoAstroLSValMap[a].second[q] = theoNorm;
				}
			}
		}

		for(int d=0; d<theoAstroLSIdxMap.size(); d++)
			for(int o=0; o<TEMP.size(); o++)
				ChiSqComps[d].second[o] = pow((theoAstroLSValMap[d].first - theoAstroLSValMap[d].second[o]), 2.0) / theoAstroLSValMap[d].first;

		for(int j=0; j<theoAstroLSIdxMap.size()-1; j++)
			if(theoAstroLSIdxMap[j].first == theoAstroLSIdxMap[j+1].first)
			{
				for(int k=0; k<TEMP.size(); k++)
					theoAstroLSValMap[j+1].second[k] += theoAstroLSValMap[j].second[k];

				theoAstroLSValMap[j].second.resize(0);
			}

		ChiSqCompsSome.resize(TEMP.size());

		for(int b=0; b<TEMP.size(); b++)
		{
			double LS = 0.0;
			double Ch = 0.0;
			double LSComp = 0.0;
			double ChComp = 0.0;

			for(int c=0; c<theoAstroLSValMap.size(); c++)
			{
				if(theoAstroLSValMap[c].second.size() == 0)
				{}
				else
				{
					double diff = theoAstroLSValMap[c].first - theoAstroLSValMap[c].second[b];
					LSComp = (diff * diff);
					ChComp = LSComp / theoAstroLSValMap[c].first;

					LS += LSComp;
					Ch += ChComp;

					ChiSqCompsSome[b].push_back(ChComp);
				}
			}

			theoAstroLSs.push_back(LS);
			alphaChiSq.push_back(Ch);
		}
	}
}

/////////////////////////////////////////////////////

void writeDataFile()
{
	int dataCounter = 0;

	ofstream out;
	out.open("Transitions");

	out.flags(ios::left);

	out.width(35);
	out << "Temperatures(K) ";
	for(int k=0; k<TEMP.size(); k++)
	{
		out.width(17);
		out << TEMP[k];
	}
	out << endl;

	out.width(35);
	out << "Ni(m^-3) "       << Ni << endl;

	out.width(35);
	out << "Ne(m^-3) "       << Ne << endl;

	out.width(35);
	out << "Residual charge " << residualCharge << endl;

	out.width(35);
	out << "Astronomical Data " << astroDataFlag;
	if(astroDataFlag) out << "\t" << astroDataFileName;
	out << endl;

	out.width(35);
	out << "Number of FF transitions " << numOfFFTrans << endl;

	out.width(35);
	out << "Number of FB transitions " << numOfFBTrans << endl;

	out.width(35);
	out << "Number of BB transitions " << numOfBBTrans << endl;

	out.width(35);
	out << "Number of all transitions " << (numOfFFTrans + numOfFBTrans + numOfBBTrans) << endl;

	out.width(35);
	out << "Normalisation Choice " << NormalisationChoice;

	if(NormalisationChoice == 1)
	{
		out << " (inside normalisation)" << endl;

		if((NormalisationEntry < 1) || (NormalisationEntry > (numOfFFTrans + numOfFBTrans + numOfBBTrans)))
		{
			cout << "Invalid value for Normalisation Entry " << endl;
			cout << "First Entry is chosen for normalisation " << endl;
			NormalisationEntry = 1;
		}

		out.width(35);
		out << "Normalisation Entry " << NormalisationEntry;
	}

	else if(NormalisationChoice == 2)
	{
		out << " (outside normalisation with emissivity data)" << endl;

		out.width(35);
		out << "Normalisation Values(W/m^3)";

		for(int q=0; q<NormalisationValues.size(); q++)
		{
			out.setf(ios::scientific);
			out  << NormalisationValues[q] << "    ";
			out.unsetf(ios::scientific);
		}
	}

	else if(NormalisationChoice == 3 || NormalisationChoice == 4)
	{
		if(NormalisationChoice == 3)
			out << " (outside normalisation with recombination coefficient data)" << endl;

		if(NormalisationChoice == 4)
			out << " (outside normalisation with interpolated recombination coefficient data using "
				<< polationPointsNum << " points interpolation scheme)" << endl;

		out.width(35);
		out << "Recombination Coefficients(cm^3/s)";

		for(int n=0; n<NormalisationValues.size(); n++)
		{
			out.setf(ios::scientific);
			out.width(17);
			out  << (alphaValues[n] / alphaValuesCovFac);
			out.unsetf(ios::scientific);
		}

		out << endl;

		////////////////////////////////////////

		out.width(35);
		out << "Interpolation Error Estimate ";

		for(int x=0; x<dys.size(); x++)
		{
			out.setf(ios::scientific);
			out.width(17);
			out  << dys[x];
			out.unsetf(ios::scientific);
		}

		out << endl;

		////////////////////////////////////////

		out.width(35);
		out << "Normalisation Wavelength(nm)";
		out << (NormalisationWL / NormWLCovFac) << endl;

		out.width(35);
		out << "Normalisation Values(W/m^3)";

		for(int q=0; q<NormalisationValues.size(); q++)
		{
			out.setf(ios::scientific);
			out  << NormalisationValues[q] << "    ";
			out.unsetf(ios::scientific);
		}
	}

	else
	{}

	out << endl;

	out.width(35);
	out << "Write Recombination Coefficients? " << RecCoef << endl;

	out.width(35);
	out << "Computing Least Squares? " << TheoAstroLSFlag;
	if(TheoAstroLSFlag) out << "\t" << TheoAstroLSFileName;
	out << endl;

	if(TheoAstroLSFlag)
	{
		if(Chi2ConfFlag)
		{
			ofstream OU;
			OU.open("TNV");

			double minAlphaChiSq = alphaChiSq.front();

			for(int s=1; s<alphaChiSq.size(); ++s)
				minAlphaChiSq = (alphaChiSq[s] <= minAlphaChiSq) ? alphaChiSq[s] : minAlphaChiSq;

			double alphaSq = minAlphaChiSq / degreesOfFreedom;

			ChiSq.resize(alphaChiSq.size());

			for(int j=0; j<alphaChiSq.size(); ++j)
				ChiSq[j] = alphaChiSq[j] / alphaSq;

			double minChiSq = ChiSq.front();
			int indexOfMinChi2 = 0;

			double Chi2Limit = degreesOfFreedom + DeltaChi2;

			for(int f=1; f<ChiSq.size(); ++f)
				if(ChiSq[f] <= minChiSq)
				{
					minChiSq = ChiSq[f];
					indexOfMinChi2 = f;
				}

			OU << "Degrees of freedom: " << degreesOfFreedom << endl;
			OU << "DeltaChi^2: "		 << DeltaChi2 << endl;
			OU << "alpha^2: "			 << alphaSq << endl;

			OU << "Temperature at minimum Chi^2: " << TEMP[indexOfMinChi2] << endl << endl;

			OU << "Temperature limits for confidence interval: " << endl;

			for(int r=1; r<ChiSq.size(); ++r)
			{
				if(	(ChiSq[r-1] >= Chi2Limit) && (ChiSq[r] <= Chi2Limit) ||
					(ChiSq[r-1] <= Chi2Limit) && (ChiSq[r] >= Chi2Limit)  )
				{
					double tempLimit = TEMP[r] + (TEMP[r-1] - TEMP[r]) * ((Chi2Limit - ChiSq[r]) / (ChiSq[r-1] - ChiSq[r]));
					OU << tempLimit << "\t\t";
				}
			}

			OU << endl << endl;

			OU.flags(ios::left);
			OU.width(10);	OU << "Temp(K)";
			OU.width(15);	OU << "Least squares";
			OU.width(15);	OU << "(alphaChi)^2";
			OU.width(15);	OU << "Chi^2";

			OU << endl << endl;

			for(int x=0; x<TEMP.size(); x++)
			{
				OU.width(10);
				OU << TEMP[x];
				OU.width(15);
				OU << theoAstroLSs[x];
				OU.width(15);
				OU << alphaChiSq[x];
				OU.width(15);
				OU << ChiSq[x];
				OU << endl;
			}

			OU.close();

			////////////////////////////// Test block
			/*
			ofstream str;
			str.open("STR");

			str << "(Idx_obs, Idx_theo, Gamma_rul)" << endl;
			str.width(12);
			str << "";
			str << "(Upper: Config, Term, 2J, Status, Gamma_r, Gamma_a)" << endl;
			str.width(12);
			str << "";
			str << "(Lower: Config, Term, 2J, Status, Gamma_r, Gamma_a)" << endl << endl;

			//start some

			int COUNT = -1;

			for(int h=0; h<theoAstroLSValMap.size(); ++h)
			{
				if(theoAstroLSValMap[h].second.size() == 0)
					continue;

				COUNT++;

				int astroIdx = theoAstroLSIdxMap[h].first;
				int theoIdx  = theoAstroLSIdxMap[h].second;

				int Idx = Transitions.size() - theoIdx;

				float tj1  = Transitions[Idx].second.first->TwoJ;
				float tj2  = Transitions[Idx].second.second->TwoJ;
				bool  p1   = Transitions[Idx].second.first->Parity;
				bool  p2   = Transitions[Idx].second.second->Parity;
				int   idx1 = Transitions[Idx].second.first->Index;
				int   idx2 = Transitions[Idx].second.second->Index;

				double Gamkj;

				if(Transitions[Idx].second.first->Status == 'F' && Transitions[Idx].second.second->Status == 'F')
					Gamkj = getContVal(p1, p2, idx1, idx2, oddResonances, evenResonances, FFGamma_rul);
				else if(Transitions[Idx].second.first->Status == 'F' && Transitions[Idx].second.second->Status == 'B')
					Gamkj = getArrOfContVal(FGamma_rul, resTwoJ, tj2, resPar, p2, idx1, idx2);
				else if(Transitions[Idx].second.first->Status == 'B' && Transitions[Idx].second.second->Status == 'B')
					Gamkj = getArrOfContVal(BGamma_rul, tj1, tj2, p1, p2, idx1, idx2);
				else
				{
					cout <<"Error!" << endl;
					exit(-1);
				}

				str << endl;
				str.flags(ios::left);
				str.width(5);  str << astroIdx;
				str.width(7);  str << theoIdx;
				str.setf(ios::scientific);
				str.precision(2);
				str.width(14); str << Gamkj;
				str.unsetf(ios::scientific);

				str << endl;

				str.width(12); str << "";
				str.width(14); str << Transitions[Idx].second.first->Config;
				str.width(5);  str << Transitions[Idx].second.first->Term;
				str.width(4);  str << Transitions[Idx].second.first->TwoJ;
				str.width(4);  str << Transitions[Idx].second.first->Status;
				str.setf(ios::scientific);
				str.precision(2);
				str.width(12);  str << Transitions[Idx].second.first->Gamma_r;
				if(Transitions[Idx].second.first->Status == 'F')
				{str.width(12);  str << Transitions[Idx].second.first->Gamma_a;}
				str.unsetf(ios::scientific);

				str << endl;

				str.width(12); str << "";
				str.width(14); str << Transitions[Idx].second.second->Config;
				str.width(5);  str << Transitions[Idx].second.second->Term;
				str.width(4);  str << Transitions[Idx].second.second->TwoJ;
				str.width(4);  str << Transitions[Idx].second.second->Status;
				str.setf(ios::scientific);
				str.width(12);  str << Transitions[Idx].second.second->Gamma_r;
				if(Transitions[Idx].second.second->Status == 'F')
				{str.width(12);  str << Transitions[Idx].second.second->Gamma_a;}
				str.unsetf(ios::scientific);
				str << endl;

				str.precision(6);
				//for(int y=0; y<TEMP.size(); ++y)
				//	str << ChiSqCompsSome[y][COUNT] / alphaSq << endl;

				str << endl;
			}

			str << endl << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" << endl << endl;
			str << "These are Chi^2 as function of temperature." << endl << endl;

			str.setf(ios::fixed);

			for(int k=0; k<ChiSqCompsSome.size(); ++k)
			{
				double sum = 0.0;
				for(int j=0; j<ChiSqCompsSome[k].size(); ++j)
					sum += ChiSqCompsSome[k][j] / alphaSq;

				str.precision(0);
				str.width(10);
				str << TEMP[k];

				str.precision(4);
				str << sum << endl;
			}

			str << endl << endl;
			str << "These are the components of Chi^2 for each contributing line." << endl << endl;

			str.precision(4);
			str.setf(ios::fixed);

			if(ChiSqCompsSome.size() > 0)
				for(int k=0; k<ChiSqCompsSome[0].size(); ++k)
				{
					double sum = 0.0;
					for(int j=0; j<ChiSqCompsSome.size(); ++j)
						sum += ChiSqCompsSome[j][k] / alphaSq;

					str.precision(6);
					str << sum << endl;
				}

			str << endl << endl;
			str << "Table of weighted LS components for contributing lines vs. temperature." << endl;	//for weighted LS components
			//str << "Table of Chi^2 components for contributing lines vs. temperature." << endl;	//for Chi^2 components
			str << "Columns stand for lines while rows for temperature." << endl << endl;

			str.precision(4);
			str.setf(ios::fixed);

			for(int d=0; d<ChiSqCompsSome.size(); ++d)
			{
				str.width(10);
				str.precision(0);
				str << TEMP[d];

				for(int j=0; j<ChiSqCompsSome[d].size(); ++j)
				{
					str.width(10);
					str.precision(4);
					str << ChiSqCompsSome[d][j];					//for weighted LS components
					//str << ChiSqCompsSome[d][j] / alphaSq;		//for Chi^2 components
				}

				str << endl;
			}

			//end some

			//start all

			for(int h=0; h<theoAstroLSIdxMap.size(); ++h)
			{
				int astroIdx = theoAstroLSIdxMap[h].first;
				int theoIdx  = theoAstroLSIdxMap[h].second;

				int Idx = Transitions.size() - theoIdx;

				float tj1  = Transitions[Idx].second.first->TwoJ;
				float tj2  = Transitions[Idx].second.second->TwoJ;
				bool  p1   = Transitions[Idx].second.first->Parity;
				bool  p2   = Transitions[Idx].second.second->Parity;
				int   idx1 = Transitions[Idx].second.first->Index;
				int   idx2 = Transitions[Idx].second.second->Index;

				double Gamkj;

				if(Transitions[Idx].second.first->Status == 'F' && Transitions[Idx].second.second->Status == 'F')
					Gamkj = getContVal(p1, p2, idx1, idx2, oddResonances, evenResonances, FFGamma_rul);
				else if(Transitions[Idx].second.first->Status == 'F' && Transitions[Idx].second.second->Status == 'B')
					Gamkj = getArrOfContVal(FGamma_rul, resTwoJ, tj2, resPar, p2, idx1, idx2);
				else if(Transitions[Idx].second.first->Status == 'B' && Transitions[Idx].second.second->Status == 'B')
					Gamkj = getArrOfContVal(BGamma_rul, tj1, tj2, p1, p2, idx1, idx2);
				else
				{
					cout <<"Error!" << endl;
					exit(-1);
				}

				str << endl;
				str.flags(ios::left);
				str.width(5);  str << astroIdx;
				str.width(7);  str << theoIdx;
				str.setf(ios::scientific);
				str.precision(2);
				str.width(14); str << Gamkj;
				str.unsetf(ios::scientific);

				str << endl;

				str.width(12); str << "";
				str.width(14); str << Transitions[Idx].second.first->Config;
				str.width(5);  str << Transitions[Idx].second.first->Term;
				str.width(4);  str << Transitions[Idx].second.first->TwoJ;
				str.width(4);  str << Transitions[Idx].second.first->Status;
				str.setf(ios::scientific);
				str.precision(2);
				str.width(12);  str << Transitions[Idx].second.first->Gamma_r;
				if(Transitions[Idx].second.first->Status == 'F')
				{str.width(12);  str << Transitions[Idx].second.first->Gamma_a;}
				str.unsetf(ios::scientific);

				str << endl;

				str.width(12); str << "";
				str.width(14); str << Transitions[Idx].second.second->Config;
				str.width(5);  str << Transitions[Idx].second.second->Term;
				str.width(4);  str << Transitions[Idx].second.second->TwoJ;
				str.width(4);  str << Transitions[Idx].second.second->Status;
				str.setf(ios::scientific);
				str.width(12);  str << Transitions[Idx].second.second->Gamma_r;
				if(Transitions[Idx].second.second->Status == 'F')
				{str.width(12);  str << Transitions[Idx].second.second->Gamma_a;}
				str.unsetf(ios::scientific);
				str << endl;

				for(int y=0; y<TEMP.size(); ++y)
					str << ChiSqComps[h].second[y] / alphaSq << endl;

				str << endl;
			}

			//end all

			str.close();
			*/
			//////////////////////////////
		}

		out.width(35);
		out << "Least Squares ";


		for(int x=0; x<TEMP.size(); x++)
		{
			out.width(17);
			out.setf(ios::fixed);
			out.precision(8);
			out << theoAstroLSs[x];
			out.unsetf(ios::fixed);
		}
		out << endl;

		if(TheoAstroLSEntFlag)
		{
			out << "Least Square Terms:" << endl;
			out << "(Observational index, Theoretical index, Terms)" << endl;
			for(int d=0; d<theoAstroLSIdxMap.size(); d++)
			{
				out.width(5);
				out << theoAstroLSIdxMap[d].first;

				out.width(10);
				out << theoAstroLSIdxMap[d].second;

				out.width(20);
				out << "";

				for(int g=0; g<theoAstroLSValMap[d].second.size(); g++)
				{
					out.setf(ios::fixed);
					out.precision(8);
					out.width(17);
					out << pow((theoAstroLSValMap[d].first - theoAstroLSValMap[d].second[g]), 2.0);
					out.unsetf(ios::fixed);
				}

				out << endl;
			}

			out << endl;

			out << "(Observational index, Theoretical index, Percent dfference {(EmN_obs - EmN_theo) / EmN_obs})" << endl;
			for(int h=0; h<theoAstroLSIdxMap.size(); h++)
			{
				out.width(5);
				out << theoAstroLSIdxMap[h].first;

				out.width(10);
				out << theoAstroLSIdxMap[h].second;

				out.width(20);
				out << "";

				for(int i=0; i<theoAstroLSValMap[h].second.size(); i++)
				{
					out.setf(ios::fixed);
					out.precision(2);
					out.width(17);
					out << 100 * (theoAstroLSValMap[h].first - theoAstroLSValMap[h].second[i]) / theoAstroLSValMap[h].first;
					out.unsetf(ios::fixed);
				}

				out << endl;
			}
		}
	}

	out << endl << endl;

	out << "Stat = Status (Bound or Free);   Exp = Experimental;   Config = Configuration;   P = Parity" << endl;
	out << "WlVac(Ang) = Wavelength in vacuum (Angstrom);   WlAir(Ang) = Wavelength in air (Angstrom)" << endl;
	out << "Em_5000 = Emissivity at temperature 5000K (W/m^3)" << endl;

	if(NormalisationChoice != 0)
		out << "EmNor_5000 = Normalised emissivity at temperature 5000K" << endl;

	if(RecCoef)
		out << "RecCoef_5000 = Recombination coefficient at temperature 5000K for vacuum wavelength (m^3/s)" << endl;

	out << endl;

	out.width(7);
	out<< "Index";

	out.width(6);
	out<< "Stat";

	out.width(6);
	out<< "Exp";

	out.width(15);
	out<< "Config";

	out.width(6);
	out<< "Term";

	out.width(3);
	out<< "2J";

	out.width(3);
	out<< "P";

	out.width(10);
	out << "   -> ";

	out.width(15);
	out<< "Config";

	out.width(6);
	out<< "Term";

	out.width(3);
	out<< "2J";

	out.width(5);
	out<< "P";

	out.width(15);
	out << "WlVac(Ang)";

	out.width(15);
	out << "WlAir(Ang)";

	for(int d=0; d<TEMP.size(); d++)
	{
		out << "Em_";
		out.width(12);
		out << TEMP[d];
	}

	if(NormalisationChoice > 0)
	{
		out << "     ";

		for(int e=0; e<TEMP.size(); e++)
		{
			out << "EmNor_";
			out.width(9);
			out << TEMP[e];
		}
	}

	if(RecCoef)
	{
		out << "     ";

		for(int e=0; e<TEMP.size(); e++)
		{
			out << "RecCoef_";
			out.width(7);
			out << TEMP[e];
		}
	}

	out << endl << endl;

	MapLines.open("MapLines");

	for(int j=Transitions.size(); j>0; j--)
	{
		bool mapline = 0;
		for(int P=0; P<theoAstroLSIdxMap.size(); ++P)
			if(theoAstroLSIdxMap[P].second == (Transitions.size() - j + 1))
				mapline = 1;

		if(mapline)
		{
			MapLines << endl << endl;

			MapLines.width(7);
			MapLines << (Transitions.size() - j + 1) ;

			MapLines << Transitions[j-1].second.first->Status;
			MapLines.width(5);
			MapLines << Transitions[j-1].second.second->Status;

			MapLines << Transitions[j-1].second.first->Exp;
			MapLines.width(5);
			MapLines << Transitions[j-1].second.second->Exp;

			//////////////////////

			MapLines.width(15);
			MapLines  << Transitions[j-1].second.first->Config;

			MapLines.width(6);
			MapLines  << Transitions[j-1].second.first->Term;

			MapLines.width(3);
			MapLines  << Transitions[j-1].second.first->TwoJ;

			MapLines.width(3);
			MapLines  << Transitions[j-1].second.first->Parity;

			MapLines.width(10);
			MapLines << "   -> ";

			MapLines.width(15);
			MapLines  << Transitions[j-1].second.second->Config;

			MapLines.width(6);
			MapLines  << Transitions[j-1].second.second->Term;

			MapLines.width(3);
			MapLines  << Transitions[j-1].second.second->TwoJ;

			MapLines.width(5);
			MapLines  << Transitions[j-1].second.second->Parity;

			MapLines << endl << endl;
		}

		out.width(7);
		out << (Transitions.size() - j + 1) ;

		out << Transitions[j-1].second.first->Status;
		out.width(5);
		out << Transitions[j-1].second.second->Status;

		out << Transitions[j-1].second.first->Exp;
		out.width(5);
		out << Transitions[j-1].second.second->Exp;

		//////////////////////

		out.width(15);
		out  << Transitions[j-1].second.first->Config;

		out.width(6);
		out  << Transitions[j-1].second.first->Term;

		out.width(3);
		out  << Transitions[j-1].second.first->TwoJ;

		out.width(3);
		out  << Transitions[j-1].second.first->Parity;

		out.width(10);
		out << "   -> ";

		out.width(15);
		out  << Transitions[j-1].second.second->Config;

		out.width(6);
		out  << Transitions[j-1].second.second->Term;

		out.width(3);
		out  << Transitions[j-1].second.second->TwoJ;

		out.width(5);
		out  << Transitions[j-1].second.second->Parity;

		double wn = Rydtoinvcm * Transitions[j-1].first;

		out.width(15);
		out.setf(ios::fixed);
		out.precision(2);
		out << WnToWlVacAng(wn);
		out.precision(6);
		out.unsetf(ios::fixed);

		out.width(15);
		if(WnToWlAirAng(wn)<2000)
			out << "----";
		else
		{
			out.setf(ios::fixed);
			out.precision(2);
			out << WnToWlAirAng(wn);
			out.precision(6);
			out.unsetf(ios::fixed);
		}

		for(int r=0; r<TEMP.size(); r++)
		{
			out.width(15);
			out.setf(ios::scientific);
			out  << Emissivities[j-1][r];
			out.unsetf(ios::scientific);

			if(mapline)
			{
				MapLines.width(20);
				MapLines.setf(ios::scientific);
				MapLines  << Emissivities[j-1][r];
				if(NormalisationChoice > 0)
				{
					MapLines.width(20);
					if(NormalisationChoice == 1)
						MapLines  << Emissivities[j-1][r] / (Emissivities[Transitions.size() - NormalisationEntry][r]);
					else
						MapLines  << Emissivities[j-1][r] / NormalisationValues[r];
				}

				MapLines.unsetf(ios::scientific);
				MapLines << endl;
			}
		}

		if(NormalisationChoice > 0)
		{
			out << "  |  ";

			if(NormalisationChoice == 1)
			{
				for(int p=0; p<TEMP.size(); p++)
				{
					double nEmiss = Emissivities[Transitions.size() - NormalisationEntry][p];
					out.width(15);
					out.setf(ios::scientific);
					out  << Emissivities[j-1][p] / nEmiss;
					out.unsetf(ios::scientific);
				}
			}
			else
			{
				for(int q=0; q<NormalisationValues.size(); q++)
				{
					double nEmiss = NormalisationValues[q];
					out.width(15);
					out.setf(ios::scientific);
					out  << Emissivities[j-1][q] / nEmiss;
					out.unsetf(ios::scientific);
				}
			}
		}

		if(RecCoef)
		{
			out << "  |  ";

			double wlvac = WnToWlVacAng(wn) / mtoAng;
			double reccoef, emis;

			for(int x=0; x<TEMP.size(); x++)
			{
				emis	= Emissivities[j-1][x];
				reccoef = emissivityToRecCoef(Ne, Ni, emis, wlvac);

				out.width(15);
				out.setf(ios::scientific);
				out  << reccoef;
				out.unsetf(ios::scientific);
			}
		}

		out << endl;

		double curr = WnToWlAirAng(wn);
		int nnn = (j-2) < 0 ? 0 : (j-2);
		double next = WnToWlAirAng(Rydtoinvcm * Transitions[nnn].first);

		if( (astroDataFlag == 1) &&
			(
			((curr <= astroDataVector[dataCounter].first) && (next >= astroDataVector[dataCounter].first)) ||
			((j==Transitions.size()) && (curr > astroDataVector[dataCounter].first)) ||
			((j==1) && (next < astroDataVector[dataCounter].first))
			)
		  )
		{
			out << endl;

			do{
                if(dataCounter >= astroDataVector.size())
                   break;

				out.width(5);
				out << astroDataVector[dataCounter].second->index;

				out.width(4);
				out << astroDataVector[dataCounter].second->dataSource;

				out.width(10);
				out << astroDataVector[dataCounter].second->Ion;

				out << "(U)";
				out.width(12);
				out  << astroDataVector[dataCounter].second->UpperTerm;

				out.width(4);
				out  << " ";

				out  << "(g2)";
				out.width(4);
				out  << astroDataVector[dataCounter].second->g2;

				out.width(10);
				out << "   -> ";

				out << "(L)";
				out.width(12);
				out  << astroDataVector[dataCounter].second->LowerTerm;

				out.width(4);
				out  << " ";

				out  << "(g1)";
				out.width(6);
				out  << astroDataVector[dataCounter].second->g1;

				out << "(o)";
				out.width(12);
				out << astroDataVector[dataCounter].second->lam_obs;

				out << "(l)";
				out.width(12);
				out << astroDataVector[dataCounter].second->lam_lab;

				out  << "(F)";
				out.width(12);
				out  << astroDataVector[dataCounter].second->F_lam;

				out  << "(I)";
				out.width(12);
				out  << astroDataVector[dataCounter].second->I_lam;

				out  << "(Mult)";
				out.width(9);
				out  << astroDataVector[dataCounter].second->Mult;

				if(NormalisationChoice > 0)
				{
					out << "  |  ";
					double norm = astroDataVector[dataCounter].second->I_lam / astroNormLine->I_lam;
					out.width(15);
					out << norm;
					out << "(relative to lam_obs = " << astroNormLine->lam_obs << " line)";
				}

				out << endl;

				dataCounter ++;

			} while((curr <= astroDataVector[dataCounter].first) && (next >= astroDataVector[dataCounter].first));

			out << endl;
		}
	}

	out.close();
}

///////////////////////////////////////////////////// TESTS SECTION

void firstTest()
{
	cout << "First Test:" << endl << endl;

	for(int x=0; x<TEMP.size(); x++)
	{
		bool flag = 0;

		for(int f=0; f<allBStates.size()-2; f++)
		{
			double NjAji  = 0.0;
			double   Aji  = 0.0;
			double NiSAik = 0.0;

			State *cur = allBStates[f].second;

			for(int g=0; g<Transitions.size(); g++)
			{
				if(Transitions[g].second.second == cur)
				{
					if(Transitions[g].second.first->Status == 'F')
					{
						Aji = getArrOfContVal(FGamma_rul, resTwoJ, Transitions[g].second.second->TwoJ,
															resPar,  Transitions[g].second.second->Parity,
															Transitions[g].second.first->Index, Transitions[g].second.second->Index);
					}

					else if(Transitions[g].second.first->Status == 'B')
					{
						Aji = getArrOfContVal(BGamma_rul, Transitions[g].second.first->TwoJ, Transitions[g].second.second->TwoJ,
														  Transitions[g].second.first->Parity,  Transitions[g].second.second->Parity,
														  Transitions[g].second.first->Index, Transitions[g].second.second->Index);
					}
					else
						Aji = -1.0E-99;

					NjAji += (Transitions[g].second.first->Nj[x] * Aji);
				}
			}

			double   Aik = 0.0;

			for(int w=0; w<Transitions.size(); w++)
			{
				if(Transitions[w].second.first == cur)
				{
					Aik = getArrOfContVal(BGamma_rul, Transitions[w].second.first->TwoJ, Transitions[w].second.second->TwoJ,
													  Transitions[w].second.first->Parity,  Transitions[w].second.second->Parity,
													  Transitions[w].second.first->Index, Transitions[w].second.second->Index);
					NiSAik += cur->Nj[x] * Aik;
				}
			}

			if((NjAji - NiSAik)/NjAji > 1.0E-14)
				flag = 1;
		}

		if(flag)
			cout << "T = " << TEMP[x] << "\tFail" <<endl;
		else
			cout << "T = " << TEMP[x] << "\tPass" <<endl;
	}

	cout << endl;
}

/////////////////////////////////////////////////////

void secondTest()
{
	cout << "Second Test:" << endl << endl;

	for(int s=0; s<TEMP.size(); s++)
	{
		double resNjAj = 0.0;
		double metaNkAki = 0.0;

		for(int f=0; f<Resonances.size(); f++)
			resNjAj += Resonances[f]->Nj[s] * Resonances[f]->Gamma_rT2;

		for(int h=0; h<allBStates.size(); h++)
		{
			if(allBStates[h].second->Gamma_r == 0.0)
			{
				State *curr = allBStates[h].second;

				for(int i=0; i<Transitions.size(); i++)
				{
					if(Transitions[i].second.second == curr)
					{
						double Nk = Transitions[i].second.first->Nj[s];

						double Aki;

						if(Transitions[i].second.first->Status == 'F')
						Aki = getArrOfContVal(FGamma_rul, resTwoJ, Transitions[i].second.second->TwoJ,
															resPar,  Transitions[i].second.second->Parity,
															Transitions[i].second.first->Index, Transitions[i].second.second->Index);

						else if(Transitions[i].second.first->Status == 'B')
						Aki = getArrOfContVal(BGamma_rul, Transitions[i].second.first->TwoJ, Transitions[i].second.second->TwoJ,
													  Transitions[i].second.first->Parity,  Transitions[i].second.second->Parity,
													  Transitions[i].second.first->Index, Transitions[i].second.second->Index);

						else
							Aki = -1.0E99;

						 metaNkAki += Nk * Aki;
					}
				}
			}
		}

		if((resNjAj - metaNkAki)/metaNkAki > 1.0E-14)
			cout << "T = " << TEMP[s] << "\tFail" <<endl;
		else
			cout << "T = " << TEMP[s] << "\tPass" <<endl;
	}

	cout << endl;
}

/////////////////////////////////////////////////////

int decayRoutes(Container* ffgamma_rul)
{
	string decRouFileName;
	stringstream decRouTemp;
	decRouTemp << "Z_DecayRoutes" << decRoutCount;
	decRouTemp >> decRouFileName;
	ofstream decayRoutes;
	decayRoutes.open(decRouFileName.c_str());

	int index;
	vector < pair <double, State*> > allStates;
	vector < pair <double, State*> > allUpperStates;

	for(int a=0; a<allFStates.size(); ++a)
		allStates.push_back(allFStates[a]);

	for(int b=0; b<allBStates.size(); ++b)
		allStates.push_back(allBStates[b]);

	bool flag = 1;

	for(int c=0; c<allStates.size(); ++c)
	{
		if(allStates[c].second->Config == DecayRoutesConfig &&
		   allStates[c].second->Term   == DecayRoutesTerm   &&
		   allStates[c].second->TwoJ   == DecayRoutes2J)
		{
			index = c;
			flag = 0;
			break;
		}
	}

	if(flag)
	{
		decayRoutes << "The state is not recognised!";
		return 0;
	}

	for(int d=index; d>=0; --d)
		allUpperStates.push_back(allStates[d]);

	//////////////////////////////////////////////////

	vector < vector < vector <int> > > ppVec;			//past-present vector
	vector < vector <int> > shortRoutes;

	ppVec.resize(2);
	ppVec[0].resize(1);
	ppVec[0][0].push_back(0);

	ULI count = 0;

	do
	{
		flag = 0;

		for(ULI e=0; e<ppVec[0].size(); ++e)
		{
			int idx = ppVec[0][e].back();
			if(idx == allUpperStates.size()-1 && ppVec[0][e].size() > 1)
			{
				shortRoutes.push_back(ppVec[0][e]);
				count++;
				continue;
			}

			else if(idx < allUpperStates.size()-1)
			{
				bool FF =0;

				for(ULI f=idx+1; f<allUpperStates.size(); ++f)
				{
					float tj1 = allUpperStates[idx].second->TwoJ;
					float tj2 = allUpperStates[f].second->TwoJ;
					bool  p1  = allUpperStates[idx].second->Parity;
					bool  p2  = allUpperStates[f].second->Parity;

					if(E1Possible(tj1, tj2, p1, p2))
					{
						flag = 1;
						FF = 1;
						vector <int> route;
						route = ppVec[0][e];
						route.push_back(f);
						ppVec[1].push_back(route);
						count++;
					}

					if( (f == allUpperStates.size()-1) && (FF == 0 ) && ppVec[0][e].size() > 1)
					{
						shortRoutes.push_back(ppVec[0][e]);
						count++;
					}
				}
			}

			if(count >= maxRoutes)
				break;
		}

		if(count >= maxRoutes)
		{
			decayRoutes << "The number of decay routes > " << maxRoutes << ". This is the first " << maxRoutes << " routes " << endl;
			ppVec.erase(ppVec.begin(), ppVec.begin()+1);
			ppVec[0].resize(maxRoutes);
			break;
		}

		else if(flag)
		{
			ppVec.erase(ppVec.begin(), ppVec.begin()+1);
			ppVec.resize(2);
			count = shortRoutes.size();
		}

		if(count >= maxRoutes)
				break;

	}while(flag);

	if(count == 0)
	{
		decayRoutes << "No decay routes are detected to the state: "
					<< allUpperStates[ppVec[0][0][0]].second->Config << " "
					<< allUpperStates[ppVec[0][0][0]].second->Term   << " "
					<< allUpperStates[ppVec[0][0][0]].second->TwoJ;
	}

	else
	{
		ULI completeR  = shortRoutes.size() >= maxRoutes ? maxRoutes : shortRoutes.size();
		ULI uCompleteR = count >= maxRoutes ? (maxRoutes - completeR) : 0;

		decayRoutes << "Number of complete decay routes in this file: " << completeR << endl;
		decayRoutes << "Number of uncomplete decay routes in this file: " << uCompleteR << endl << endl;

		decayRoutes << "Data for each state: configuration, term, 2J." << endl;
		decayRoutes << "For resonances, Saha capture and radiative population are added." << endl;
		decayRoutes << "The temperature dependent data are for T = " << TEMP[0] << "K." << endl;

		if(completeR > 0)
		{
			decayRoutes << endl << "Complete routes:" << endl << endl;

			for(ULI u=0; u<completeR; ++u)
			{
				for(ULI v=0; v<shortRoutes[u].size(); ++v)
				{
					if(v > 0)
						decayRoutes	<< " <- ";

					decayRoutes << allUpperStates[shortRoutes[u][v]].second->Config << " "
								<< allUpperStates[shortRoutes[u][v]].second->Term   << " "
								<< allUpperStates[shortRoutes[u][v]].second->TwoJ;

					if(allUpperStates[shortRoutes[u][v]].second->Status == 'F')
					{
						decayRoutes << " (";
						decayRoutes.setf(ios::scientific);
						decayRoutes.precision(2);
						decayRoutes.width(10);
						decayRoutes << allUpperStates[shortRoutes[u][v]].second->decRoutSahaCap;

						decayRoutes << ",";
						decayRoutes.precision(2);
						decayRoutes.width(10);

						if(v < shortRoutes[u].size()-1)
						{
							float tj1  = allUpperStates[shortRoutes[u][v+1]].second->TwoJ;		//upper
							float tj2  = allUpperStates[shortRoutes[u][v]].second->TwoJ;		//lower
							bool  p1   = allUpperStates[shortRoutes[u][v+1]].second->Parity;
							bool  p2   = allUpperStates[shortRoutes[u][v]].second->Parity;
							int   idx1 = allUpperStates[shortRoutes[u][v+1]].second->Index;
							int   idx2 = allUpperStates[shortRoutes[u][v]].second->Index;

							bool pos = E1Possible(tj1, tj2, p1, p2);
							double radPop;

							if(pos == true)
							{
								double Nk    = allUpperStates[shortRoutes[u][v+1]].second->Nj[0];
								double Gamkj = getContVal(p1, p2, idx1, idx2, oddResonances, evenResonances, ffgamma_rul);
								double Gamjr  = allUpperStates[shortRoutes[u][v]].second->Gamma_r;
								double Gamja = allUpperStates[shortRoutes[u][v]].second->Gamma_a;

								radPop = (Nk * Gamkj) / (Gamjr + Gamja);

								decayRoutes << radPop;
							}

							else
								decayRoutes << "NP";			// E1 not possible
						}

						else
							decayRoutes << "NA";				// Not applicable

						decayRoutes.unsetf(ios::scientific);
						decayRoutes << ")";
					}
				}

				decayRoutes << endl;
			}
		}

		if(uCompleteR > 0)
		{
			decayRoutes << endl << "Uncomplete routes:" << endl << endl;

			for(ULI s=0; s<uCompleteR; ++s)
			{
				for(ULI t=0; t<ppVec[0][s].size(); ++t)
				{
					if(t > 0)
						decayRoutes	<< " <- ";

					decayRoutes << allUpperStates[ppVec[0][s][t]].second->Config << " "
								<< allUpperStates[ppVec[0][s][t]].second->Term   << " "
								<< allUpperStates[ppVec[0][s][t]].second->TwoJ;
				}

				decayRoutes << endl;
			}
		}
	}

	/////////////////////////////////////////////////

	// Repetition test block
	/*
	for(ULI w=0; w<ppVec[0].size()-1; ++w)
		for(ULI x=w+1; x<ppVec[0].size(); ++x)
		{
			bool F = 0;

			if(ppVec[0][w].size() == ppVec[0][x].size())
			{
				for(int q=0; q<ppVec[0][w].size(); ++q)
				{
					if(ppVec[0][w][q] != ppVec[0][x][q])
					{
						F = 1;
						break;
					}
				}
			}

			if(!F)
            {
				cout << "Decay repetition test failed if all decay routes are complete!" << endl;
				return 0;
            }
		}
		cout << "Decay repetition test passed successfully!" << endl << endl;
	*/

	allStates.clear();
	allUpperStates.clear();
	ppVec.clear();

	return 0;
}

/////////////////////////////////////////////////////

int main()
{
	clock_t startProgram(clock());

///////////////////////////////////////////////////// READING INPUT DATA FILES

	readInputFile();										// Reading input file
	readResEmisFile();										// Reading and storing ResEmis data
	readELEVEmisFile();										// Reading and storing ELEVEmis data
	readFVALUEFile();										// Reading and storing FVALUE data

///////////////////////////////////////////////////// FREE-FREE SECTION

	allFStates = arrangeAllFStates(Resonances, allFStates); // Arrange all free states in desending order

	if(numOfOddRes > 0 && numOfEvenRes > 0)
	{
		FFGamma_rul = findFFGamma_rul(FFfValues, Resonances, FFGamma_rul);			// Find Gamma_rul for free states
		findFFGamma_r(allFStates, FFGamma_rul);
	}

///////////////////////////////////////////////////// FREE-BOUND SECTION

	FGamma_rul = findFGamma_rul(numOfBSym, Resonances, ResfValues, ResPhotonEnergy, FGamma_rul); // Find Gamma_rul for FB
	findFGamma_r(numOfAllRes, numOfBSym, Resonances, FGamma_rul);								 // Find Gamma_r for resonances
	findbu(numOfAllRes, Resonances);															 // Find bu for resonances
	findFFNj(allFStates, TEMP, FFGamma_rul);

///////////////////////////////////////////////////// BOUND-BOUND SECTION

	BGamma_rul = findBGamma_rul(BFVALUES, allBSym, BGamma_rul);			// Find Gamma_rul for bound states
	allBStates = arrangeAllBStates(allBSym, allBStates);				// Arrange all bound states in desending order
	findBGamma_r(allBStates, BGamma_rul);								// Find Gamma_r for bound states
	findBNj(Resonances, allBStates, TEMP, BGamma_rul);						// Find Nj for bound states

///////////////////////////////////////////////////// MIXED SECTION

	Transitions = findTransitions(Resonances, allBStates, Transitions, allFStates);		// Find possible transitions
	Emissivities = findEmissivities(Transitions, FFGamma_rul, FGamma_rul, BGamma_rul, Emissivities); // Find emissivities
	if(astroDataFlag == 1)																// Read asronomical data from AstroData file
		readAstroData();
	if(TheoAstroLSFlag == 1)															// Perform LS
		TheoAstroLS();
	writeDataFile();																	// Write data to output file

///////////////////////////////////////////////////// TESTS SECTION

	if(testFlag == 1 || testFlag == 3)
		firstTest();
	if(testFlag == 2 || testFlag == 3)
		secondTest();

///////////////////////////////////////////////////// DECAY ROUTES

	if(DecayRoutesFlag)
		decayRoutes(FFGamma_rul);

///////////////////////////////////////////////////// END OF PROGRAM

	double elapsedProgram = (double)(clock() - startProgram) / CLOCKS_PER_SEC;
	int min = (int)(elapsedProgram/60.0);
	int sec = (int)(elapsedProgram - min * 60.0);
	double dmsec = (elapsedProgram - (min*60.0) - (sec*1.0)) * 1000.0;
	int msec = (int) dmsec;

	cout   <<"Program Duration (min:sec:msec):      " << min << ":" << sec << ":" << msec << endl << endl << endl;
	cout << "Press any key to exit" << endl;
	getchar();

	return 0;
}