#include "utilities.h"

//---------------------
bool my_mkdir(string s)
	{
	bool result =false;
	string outstr;
	outstr = "mkdir -p " + s;
	//cout << outstr << endl;
	int res = system(outstr.c_str());
	if(res != 0)
		{
		cerr << "directory creation error! - directory: " << s << endl;
		result =false;
		}
	else
		{
		result = true;
		}
	return result;
	/*
	DIR *pdir;
	string dir_str = s;
	pdir = opendir(dir_str.c_str());
	if (!pdir)
		{
		mkdir(dir_str.c_str(), 0777);
		}
	else
		{
		closedir(pdir);
		}
	*/
	}
//---------------------
Array<double,3> calculate_simpson_13_weights(int Nx,int Ny,int Nz,double Dx,double Dy,double Dz)
	{
	if ( (Nx<=0) || (Ny <=0) || (Nz <=0) )
		{
		cerr << "input error: from method calculate_simpson_13_weights\n\tnumber of input points should be positive!\n\taborting!\n";
		exit(-1);
		}
	if ( (Nx%2 == 0) || (Ny%2 == 0) || (Nz%2 == 0))
		{
		cerr << "input error: from method calculate_simpson_13_weights\n\tnumber of input points should be odd!\n\taborting!\n";
		exit(-1);
		}
	Array<double,1> simpson_13_weights_x(Nx);
	Array<double,1> simpson_13_weights_y(Ny);
	Array<double,1> simpson_13_weights_z(Nz);

	simpson_13_weights_x(Range::all()) = 1.0*Dx/3.0;
	if(Nx >= 2)
		{
		simpson_13_weights_x(Range(1,Nx-2,2)) = 4.0*Dx/3.0;
		}
	if(Nx >= 3)
		{
		simpson_13_weights_x(Range(2,Nx-3,2)) = 2.0*Dx/3.0;
		}

	simpson_13_weights_y(Range::all()) = 1.0*Dy/3.0;
	if(Ny >= 2)
		{
		simpson_13_weights_y(Range(1,Ny-2,2)) = 4.0*Dy/3.0;
		}
	if(Ny >= 3)
		{
		simpson_13_weights_y(Range(2,Ny-3,2)) = 2.0*Dy/3.0;
		}

	simpson_13_weights_z(Range::all()) = 1.0*Dz/3.0;
	if(Nz >= 2)
		{
		simpson_13_weights_z(Range(1,Nz-2,2)) = 4.0*Dz/3.0;
		}
	if(Nz >= 3)
		{
		simpson_13_weights_z(Range(2,Nz-3,2)) = 2.0*Dz/3.0;
		}

	Array<double,3> simpson_13_weights_xyz(Nx,Ny,Nz);

	simpson_13_weights_xyz = simpson_13_weights_x(tensor::i) *
	                         simpson_13_weights_y(tensor::j) * simpson_13_weights_z(tensor::k);
	return simpson_13_weights_xyz;
	}
//---------------------
Array<double,2> calculate_simpson_13_weights(int Ntheta,int Nphi,double Dtheta,double Dphi)
	{
	if ( (Ntheta<=0) || (Nphi <=0) )
		{
		cerr << "input error: from method calculate_simpson_13_weights\n\tnumber of input points should be positive!\n\taborting!\n";
		exit(-1);
		}
	if ( (Ntheta%2 == 0) || (Nphi%2 == 0) )
		{
		cerr << "input error: from method calculate_simpson_13_weights\n\tnumber of input points should be odd!\n\taborting!\n";
		exit(-1);
		}
	Array<double,1> simpson_13_weights_t(Ntheta);
	Array<double,1> simpson_13_weights_p(Nphi);

	simpson_13_weights_t(Range::all()) 				= 1.0*Dtheta/3.0;
	if(Ntheta >= 2)
		{
		simpson_13_weights_t(Range(1,Ntheta-2,2))  = 4.0*Dtheta/3.0;
		}
	if(Ntheta >=3 )
		{
		simpson_13_weights_t(Range(2,Ntheta-3,2))  = 2.0*Dtheta/3.0;
		}

	simpson_13_weights_p(Range::all()) 				= 1.0*Dphi/3.0;
	if(Nphi >= 2)
		{
		simpson_13_weights_p(Range(1,Nphi-2,2)) 		= 4.0*Dphi/3.0;
		}
	if(Nphi >= 3)
		{
		simpson_13_weights_p(Range(2,Nphi-3,2)) 		= 2.0*Dphi/3.0;
		}

	Array<double,2> simpson_13_weights_tp(Ntheta,Nphi);

	simpson_13_weights_tp = simpson_13_weights_t(tensor::i) *
	                        simpson_13_weights_p(tensor::j);
	return simpson_13_weights_tp;
	}
//-----------------------
void file_read_error(string s)
	{
	cerr << "File error: Unable to read " << s.c_str() << "!\n";
	exit(-1);
	}
//-----------------------
Array<double,1> load_polynomial_coeffs(string filename, int poly_fit_order)
	{
	Array<double,1> polycoeffs(poly_fit_order+1);
	ifstream fin(filename.c_str());
	if(!fin)
		{
		file_read_error(filename);
		}
	for(int i = 0; i<poly_fit_order+1; i++)
		{
		fin >> polycoeffs(i);
		}
	fin.close();
	return polycoeffs;
	}
//-------------------------
Array<complex<double>,2> make_angular_weights( Array<double,1> ampl_coeffs,
    Array<double,1> phase_coeffs, focus3d_parameters param)
	{
	complex<double> I = complex<double>(0,1);

	int poly_fit_order = ampl_coeffs.extent(0)-1;
	Array<complex<double>,2> aweights(param.Ntheta,param.Nphi);
	Array<double,1> theta(param.Ntheta);
	double Dtheta = (param.theta_max-param.theta_min)/(param.Ntheta-1);
	theta = param.theta_min + tensor::i*Dtheta;

	complex<double> temp = 0;
	double temp_ampl;
	double temp_phase;

	for (int l =0; l<param.Ntheta;l++)
		{
		temp = 0;
		temp_ampl = 0;
		temp_phase = 0;
		for(int i =0; i<poly_fit_order+1;i++)
			{
			temp_ampl += ampl_coeffs(i)*pow(theta(l), poly_fit_order-i);
			temp_phase += phase_coeffs(i)*pow(theta(l), poly_fit_order-i);
			}
		temp = temp_ampl*exp(I*temp_phase);
		aweights(l,Range::all()) = temp;
		//cout << theta(l) << "\t" 	<< real(temp) << "\t" << imag(temp) << "\n";
		//cout << abs(temp) << endl;
		}
	return aweights;
	}
//----------------------------
int determine_poly_fit_order(string filename)
	{
	fstream fin(filename.c_str());

	if(!fin) {file_read_error(filename);}
	fin.seekg(0,ios::end);
	int file_length = fin.tellg(), current_pos = 0;
	fin.seekg(0,ios::beg);

	double temp;
	int index = 0;
	while( current_pos != file_length )
		{
		fin >> temp;
		current_pos = fin.tellg();
		if(current_pos == -1)
			{
			break;
			}
		++index;
		}
	fin.close();

	int poly_fit_order = index-1;
	cout << "File: " << filename << " -  polynomial fit order " << poly_fit_order << endl;
	return poly_fit_order;
	}
//----------------------------
Array<complex<double>,2> make_angular_weights( string ampl_coeff_file,
    string phase_coeff_file, focus3d_parameters param)
	{
	complex<double> I = complex<double>(0,1);

	int temp1 = determine_poly_fit_order(ampl_coeff_file);
	int temp2 = determine_poly_fit_order(phase_coeff_file);

	if(temp1 != temp2)
		{
		cerr << "Input error: Polynomial coefficient files are not of equal lengths!\n";
		exit(-1);
		}

	int poly_fit_order = temp1;

	Array<double,1> poly_coeffs_ampl(poly_fit_order+1), poly_coeffs_phase(poly_fit_order+1);

	poly_coeffs_ampl 	= load_polynomial_coeffs(ampl_coeff_file,poly_fit_order); // load the ampl polynomial coeffs file
	poly_coeffs_phase	= load_polynomial_coeffs(phase_coeff_file,poly_fit_order); // load the phase polynomial coeffs file
	cout << "Amplitude coeffs:\n" << poly_coeffs_ampl << endl
	<< "Phase coeffs:\n" << poly_coeffs_phase << endl;

	Array<complex<double>,2> aweights(param.Ntheta,param.Nphi);
	Array<double,1> theta(param.Ntheta);
	double Dtheta = (param.theta_max-param.theta_min)/(param.Ntheta-1);
	theta = param.theta_min + tensor::i*Dtheta;

	double temp_ampl;
	double temp_phase;

	for (int l =0; l<param.Ntheta;l++)
		{
		temp_ampl = 0;
		temp_phase = 0;
		for(int i =0; i< poly_fit_order+1;i++)
			{
			temp_ampl += poly_coeffs_ampl(i)*pow(theta(l), poly_fit_order-i);
			temp_phase += poly_coeffs_phase(i)*pow(theta(l), poly_fit_order-i);
			}
		aweights(l,Range::all()) = temp_ampl*exp(I*temp_phase);
		}
	return aweights;
	}
//----------------------------
int get_number_of_data_points(string filename)
	{
	ifstream fin(filename.c_str());

	if(!fin) {file_read_error(filename);}
	fin.seekg(0,ios::end);
	int file_length = fin.tellg(), current_pos = 0;
	fin.seekg(0,ios::beg);

	double temp;
	int index = 0;
	while( current_pos != file_length )
		{
		fin >> temp;
		current_pos = fin.tellg();
		if(current_pos == -1)
			{
			if(!fin.good() ) 	{	fin.clear(); }
			break;
			}
		++index;
		}
	fin.close();
	return index;
	}
//----------------------------
void interpolate_data(string filename, Array<double,1> theta, Array<complex<double>,1>& output)
	{
	ifstream fin(filename.c_str());

	if(!fin) {file_read_error(filename);}
	fin.seekg(0,ios::end);
	int file_length = fin.tellg(), current_pos = 0;
	fin.seekg(0,ios::beg);

	double temp;
	int index = 0;
	while( current_pos != file_length )
		{
		fin >> temp;
		current_pos = fin.tellg();
		if(current_pos == -1)
			{
			if(!fin.good() ) 	{	fin.clear(); }
			break;
			}
		++index;
		}
	if( (index%3) != 0  )
		{
		cerr << "File error! " << filename << "should contain three columns of data corresponding to theta, amplitude and phase of G(theta) respectively! Exiting!";
		exit(-1);
		}

	int numpoints = index/3;

	Array<double,1> theta_in(numpoints), ampl_in(numpoints), phase_in(numpoints);
	Array<double,1> ampl_out(theta.extent(0)), phase_out(theta.extent(0));

	fin.seekg(0,ios::beg);
	for(int i =0; i<numpoints;i++)
		{
		fin >> theta_in(i);
		fin >> ampl_in(i);
		fin >> phase_in(i);
		}
	fin.close();

	int low_index = INT_MIN, up_index = INT_MAX;
	for(int i =0; i<theta.extent(0); i++)
		{
		if( (first( theta_in > theta(i) ) != INT_MIN) && (last( theta_in < theta(i) ) != INT_MAX) )
			{
			//cout << theta(i) << "\t"<<"I am in low-up index\n" << "\n";
			low_index = first( theta_in > theta(i) );
			up_index = last( theta_in < theta(i) );
			//cout << low_index << "\t" << up_index << "\n";
			ampl_out(i) = ampl_in(low_index) + (theta(i)-theta_in(low_index)) * (ampl_in(up_index) - ampl_in(low_index)) / (theta_in(up_index) - theta_in(low_index));

			phase_out(i) = phase_in(low_index) + (theta(i)-theta_in(low_index)) * (phase_in(up_index) - phase_in(low_index)) / (theta_in(up_index) - theta_in(low_index));
			}
		else if ( (first( theta_in == theta(i) ) != INT_MIN) )
			{
			//cout << theta(i) << "\t" << "I am in low index\n" << "\n";
			//cout << low_index << "\t" << up_index << "\n";
			low_index = first( theta_in == theta(i) );
			ampl_out(i) = ampl_in(low_index);
			phase_out(i) = phase_in(low_index);
			}
		else if ( (last( theta_in == theta(i) ) != INT_MAX) )
			{
			//cout << theta(i) << "\t" << "I am in up index\n" << "\n";
			//cout << low_index << "\t" << up_index << "\n";
			up_index = last( theta_in == theta(i) );
			ampl_out(i) = ampl_in(up_index);
			phase_out(i) = phase_in(up_index);
			}
		else
			{
			cerr << "Data error! The range of input data in "<< filename <<" is not sufficient for interpolation!\n";
			cerr << theta(i) << " is not within the limits of the input " << theta_in;
			cerr << "Exiting!\n";
			exit(-1);
			}
		//cout << low_index << "\t" << up_index << "\t";
		//cout << theta(i) << "\t" << ampl_out(i) << "\t" << phase_out(i) << "\n";
		}
	complex<double> I = complex<double>(0.0,1.0);
	output = ampl_out*exp(I*phase_out);
	}




