#include "focus3dws.h"

//---------------------
bool focus3dws::test_parameters(focus3d_parameters initial, focus3d_parameters scrambled)
	{
	return (initial.test() && scrambled.test());
	}
//---------------------
void focus3dws::load_parameters(focus3d_parameters initial, focus3d_parameters scrambled)
	{
	if ( (initial.lambda == scrambled.lambda) && (initial.f == scrambled.f) && (initial.Nphi == scrambled.Nphi) && (initial.phi_min == scrambled.phi_min) && (initial.phi_max == scrambled.phi_max) )
		{
		param_ini.Nx 				= initial.Nx;
		param_ini.Ny 				= initial.Ny;
		param_ini.Nz 				= initial.Nz;
		param_ini.Ntheta 		= initial.Ntheta;
		param_ini.Nphi 			= initial.Nphi;
		param_ini.xmin 			= initial.xmin;
		param_ini.xmax 			= initial.xmax;
		param_ini.ymin 			= initial.ymin;
		param_ini.ymax 			= initial.ymax;
		param_ini.zmin 			= initial.zmin;
		param_ini.zmax 			= initial.zmax;
		param_ini.theta_min	= initial.theta_min;
		param_ini.theta_max	= initial.theta_max;
		param_ini.phi_min 	= initial.phi_min;
		param_ini.phi_max 	= initial.phi_max;
		param_ini.f 				= initial.f;
		param_ini.lambda 		= initial.lambda;

		param_scr.Nx 				= scrambled.Nx;
		param_scr.Ny 				= scrambled.Ny;
		param_scr.Nz 				= scrambled.Nz;
		param_scr.Ntheta 		= scrambled.Ntheta;
		param_scr.Nphi 			= scrambled.Nphi;
		param_scr.xmin 			= scrambled.xmin;
		param_scr.xmax 			= scrambled.xmax;
		param_scr.ymin 			= scrambled.ymin;
		param_scr.ymax 			= scrambled.ymax;
		param_scr.zmin 			= scrambled.zmin;
		param_scr.zmax 			= scrambled.zmax;
		param_scr.theta_min	= scrambled.theta_min;
		param_scr.theta_max	= scrambled.theta_max;
		param_scr.phi_min 	= scrambled.phi_min;
		param_scr.phi_max 	= scrambled.phi_max;
		param_scr.f 				= scrambled.f;
		param_scr.lambda 		= scrambled.lambda;

		if(param_scr.Nx != 1)
			{
			Dx 	= (param_scr.xmax-param_scr.xmin)/double(param_scr.Nx-1);
			}
		else
			{
			Dx = 0;
			}
		if(param_scr.Ny != 1)
			{Dy = (param_scr.ymax-param_scr.ymin)/double(param_scr.Ny-1);
			}
		else
			{
			Dy = 0;
			}
		if(param_scr.Nz != 1)
			{
			Dz	= (param_scr.zmax-param_scr.zmin)/double(param_scr.Nz-1);
			}
		else
			{
			Dz = 0;
			}
		if(param_scr.Ntheta != 1)
			{
			Dtheta_scr 	= (param_scr.theta_max-param_scr.theta_min)/double(param_scr.Ntheta-1);
			}
		else
			{
			Dtheta_scr = 0;
			}
		if(param_scr.Nphi != 1)
			{
			Dphi_scr 	= (param_scr.phi_max-param_scr.phi_min)/double(param_scr.Nphi-1);
			}
		else
			{
			Dphi_scr = 0;
			}

		if(param_ini.Ntheta != 1)
			{
			Dtheta_ini 	= (param_ini.theta_max-param_ini.theta_min)/double(param_ini.Ntheta-1);
			}
		else
			{
			Dtheta_ini = 0;
			}
		if(param_ini.Nphi != 1)
			{
			Dphi_ini 	= (param_ini.phi_max-param_ini.phi_min)/double(param_ini.Nphi-1);
			}
		else
			{
			Dphi_ini = 0;
			}
		outputstream.str(string());
		outputstream << "Input parameters:\n******************\nlambda       = " << param_scr.lambda
		<< " µm\nfocal length = " << param_scr.f
		<< " µm\nNx     = " << param_scr.Nx 			<< "\txmin      = " << param_scr.xmin 			<< " µm\txmax      = " << param_scr.xmax      << " µm\tDx     = " << Dx
		<< " µm\nNy     = " << param_scr.Ny 			<< "\tymin      = " << param_scr.ymin 			<< " µm\txmax      = " << param_scr.ymax      << " µm\tDy     = " << Dy
		<< " µm\nNz     = " << param_scr.Nz 			<< "\tzmin      = " << param_scr.zmin 			<< " µm\tzmax      = " << param_scr.zmax      << " µm\tDz     = " << Dz
		<< " µm\nNtheta = " << param_scr.Ntheta 	<< "\ttheta_min = " << param_scr.theta_min << " rad\ttheta_max = " << param_scr.theta_max << " rad\tDtheta = " << Dtheta_scr
		<< " rad\nNphi   = " << param_scr.Nphi 		<< "\tphi_min   = " << param_scr.phi_min 		<< " rad\tphi_max   = " <<  param_scr.phi_max 	<< " rad\tDphi   = " << Dphi_scr
		<< " rad\nNtheta' = " << param_ini.Ntheta 	<< "\ttheta_min' = " << param_ini.theta_min << " rad\ttheta_max' = " << param_ini.theta_max << " rad\tDtheta' = " << Dtheta_ini
		<< " rad\nNphi'   = " << param_ini.Nphi 		<< "\tphi_min'   = " << param_ini.phi_min 		<< " rad\tphi_max'   = " <<  param_ini.phi_max 	<< " rad\tDphi'   = " << Dphi_ini << " rad\n\n";

		show_output(outputstream.str());
		}
	else
		{
		outputstream.str(string());
		outputstream << "Input error! phi_min, phi_max, Nphi, lambda and f should be same for both the focus3d_parameters!\n" << "Current inputs are \n\t initial\t scrambled\nNphi   = \t" << initial.Nphi << "\t" << scrambled.Nphi << "\nlambda = \t" << initial.lambda << "\t" << scrambled.lambda << "\nf      = \t" << initial.f << "\t" << scrambled.f << "\nphi_max= \t" << initial.phi_max << "\t" << scrambled.phi_max << "\nphi_min= \t" << initial.phi_min << "\t" << scrambled.phi_min << "\n";
		show_error(outputstream.str());
		}
	}
//---------------------
focus3dws::focus3dws(focus3d_parameters initial, focus3d_parameters scrambled)
	{
	I = complex<double>(0.0,1.0);
	reset_parameters(initial,scrambled);
	}
//---------------------
void focus3dws::reset_parameters(focus3d_parameters initial, focus3d_parameters scrambled)
	{
	if (test_parameters(initial,scrambled) == true)
		{
		load_parameters(initial,scrambled);
		initialize();
		}
	else
		{
		errorstream.str(string());
		errorstream << "Input error: Check your parameter input!\n";
		show_error(errorstream.str());
		}
	}
//----------------------
focus3dws::~focus3dws()
	{
	x.free();
	y.free();
	z.free();
	theta_ini.free();
	theta_scr.free();
	phi_ini.free();
	phi_scr.free();
	simpson_13_weights.free();
	theta2_ini.free();
	theta2_scr.free();
	phi2_ini.free();
	phi2_scr.free();
	E_inc_x.free();
	E_inc_y.free();
	E_inc_z.free();
	ini_integrand_amplitude_x.free();
	ini_integrand_amplitude_y.free();
	ini_integrand_amplitude_z.free();
	scr_integrand_amplitude_x.free();
	scr_integrand_amplitude_y.free();
	scr_integrand_amplitude_z.free();
	exp_arg1.free();
	exp_arg2.free();
	exp_arg3.free();
	integrand_exp_term.free();
	E_focus_x.free();
	E_focus_y.free();
	E_focus_z.free();
	}
//----------------------
void focus3dws::populate_x()
	{
	x.resize(param_scr.Nx);
	x  = param_scr.xmin + tensor::i*Dx;
	}
//----------------------
void focus3dws::populate_y()
	{
	y.resize(param_scr.Ny);
	y  = param_scr.ymin + tensor::i*Dy;
	}
//----------------------
void focus3dws::populate_z()
	{
	z.resize(param_scr.Nz);
	z  = param_scr.zmin + tensor::i*Dz;
	}
//----------------------
void focus3dws::populate_theta()
	{
	theta_scr.resize(param_scr.Ntheta);
	theta_scr  = param_scr.theta_min + tensor::i*Dtheta_scr;
	theta2_scr.resize(param_scr.Ntheta,param_scr.Nphi);
	for (int l =0; l<param_scr.Ntheta;l++)
		{
		theta2_scr(l,Range::all()) = theta_scr(l);
		}

	theta_ini.resize(param_ini.Ntheta);
	theta_ini  = param_ini.theta_min + tensor::i*Dtheta_ini;
	theta2_ini.resize(param_ini.Ntheta,param_ini.Nphi);
	for (int l =0; l<param_ini.Ntheta;l++)
		{
		theta2_ini(l,Range::all()) = theta_ini(l);
		}
	}
//----------------------
void focus3dws::populate_phi()
	{
	phi_scr.resize(param_scr.Nphi);
	phi_scr  = param_scr.phi_min + tensor::i*Dphi_scr;
	phi2_scr.resize(param_scr.Ntheta,param_scr.Nphi);
	for (int m =0; m<param_scr.Nphi;m++)
		{
		phi2_scr(Range::all(),m) = phi_scr(m);
		}

	phi_ini.resize(param_ini.Nphi);
	phi_ini  = param_ini.phi_min + tensor::i*Dphi_ini;
	phi2_ini.resize(param_ini.Ntheta,param_ini.Nphi);
	for (int m =0; m<param_ini.Nphi;m++)
		{
		phi2_ini(Range::all(),m) = phi_ini(m);
		}
	}
//----------------------
void focus3dws::calculate_integrand_amplitudes()
	{
	// !!!!!!! call this function only after populating theta2, phi2, and E_inc_x,y,z
	ini_integrand_amplitude_x.resize(param_ini.Ntheta,param_ini.Nphi);
	ini_integrand_amplitude_y.resize(param_ini.Ntheta,param_ini.Nphi);
	ini_integrand_amplitude_z.resize(param_ini.Ntheta,param_ini.Nphi);

	ini_integrand_amplitude_x =
	  sqrt(cos(theta2_ini))
	  * ( + E_inc_x * ( sqr(cos(phi2_ini))*cos(theta2_ini) + sqr(sin(phi2_ini)) )
	      + E_inc_y * ( -(1-cos(theta2_ini))*sin(phi2_ini)*cos(phi2_ini) )
	      + E_inc_z * ( -cos(phi2_ini)*sin(theta2_ini) )
	    );
	ini_integrand_amplitude_y =
	  sqrt(cos(theta2_ini))
	  * ( + E_inc_x * ( -(1-cos(theta2_ini))*sin(phi2_ini)*cos(phi2_ini) )
	      + E_inc_y * ( sqr(sin(phi2_ini))*cos(theta2_ini) + sqr(cos(phi2_ini)) )
	      + E_inc_z * ( -sin(phi2_ini)*sin(theta2_ini) )
	    );
	ini_integrand_amplitude_z =
	  sqrt(cos(theta2_ini))
	  * ( + E_inc_x * ( sin(theta2_ini)*cos(phi2_ini) )
	      + E_inc_y * ( sin(theta2_ini)*sin(phi2_ini) )
	      + E_inc_z * ( cos(theta2_ini) )
	    );

	scr_integrand_amplitude_x.resize(param_scr.Ntheta,param_scr.Nphi);
	scr_integrand_amplitude_y.resize(param_scr.Ntheta,param_scr.Nphi);
	scr_integrand_amplitude_z.resize(param_scr.Ntheta,param_scr.Nphi);
	scr_integrand_amplitude_x = 0.0;
	scr_integrand_amplitude_y = 0.0;
	scr_integrand_amplitude_z = 0.0;

	for(int l = 0; l < param_scr.Ntheta; l++)
		{
		for(int lp = 0; lp < param_ini.Ntheta; lp++)
			{
			scr_integrand_amplitude_x(l,Range::all()) += G_ttp(l,lp)*ini_integrand_amplitude_x(lp,Range::all());
			scr_integrand_amplitude_y(l,Range::all()) += G_ttp(l,lp)*ini_integrand_amplitude_y(lp,Range::all());
			scr_integrand_amplitude_z(l,Range::all()) += G_ttp(l,lp)*ini_integrand_amplitude_z(lp,Range::all());
			}
		}

	scr_integrand_amplitude_x = sin(theta2_scr)*scr_integrand_amplitude_x;
	scr_integrand_amplitude_y = sin(theta2_scr)*scr_integrand_amplitude_y;
	scr_integrand_amplitude_z = sin(theta2_scr)*scr_integrand_amplitude_z;
	}

//----------------------
void focus3dws::calculate_integrand_exp_args()
	{
	// !!!!!!! call this function only after populating theta2 and phi2

	exp_arg1.resize(param_scr.Ntheta,param_scr.Nphi);
	exp_arg2.resize(param_scr.Ntheta,param_scr.Nphi);
	exp_arg3.resize(param_scr.Ntheta,param_scr.Nphi);

	exp_arg1 = I * (2.0*M_PI/param_scr.lambda) * sin(theta2_scr) * cos(phi2_scr); // i.k.sin(t).cos(f)
	exp_arg2 = I * (2.0*M_PI/param_scr.lambda) * sin(theta2_scr) * sin(phi2_scr);	// i.k.sin(t).sin(f)
	exp_arg3 = I * (2.0*M_PI/param_scr.lambda) * cos(theta2_scr);									// i.k.cos(t)

	prefactor = 0;
	prefactor = I * (1.0/param_scr.lambda) * param_scr.f * exp(I * param_scr.f * (2.0*M_PI/param_scr.lambda));
	}
//----------------------
void focus3dws::calculate_simpson_13_weights()
	{
	simpson_13_weights.resize(param_scr.Ntheta,param_scr.Nphi);
	double a,b;
	(Dtheta_scr == 0) ? (a = 3.0) : (a = Dtheta_scr); // Dtheta =0 only when Ntheta = 1
	(Dphi_scr == 0) ? (b = 3.0) : (b = Dphi_scr); // Dphi =0 only when Nphi = 1
	//cout << "a = " << a << "\tb = " << b <<  "\n";
	simpson_13_weights = ::calculate_simpson_13_weights(param_scr.Ntheta,param_scr.Nphi,a,b);
	}
//----------------------
void focus3dws::initialize_E_incs()
	{
	E_inc_x.resize(param_ini.Ntheta,param_ini.Nphi); E_inc_x = 1.0; // x-component - default input polarization
	E_inc_y.resize(param_ini.Ntheta,param_ini.Nphi); E_inc_y = 0.0; // y-component
	E_inc_z.resize(param_ini.Ntheta,param_ini.Nphi); E_inc_y = 0.0; // z-component
	}
//----------------------
void focus3dws::set_incident_field(Array<complex<double>,2> E_x,
                                   Array<complex<double>,2> E_y,
                                   Array<complex<double>,2> E_z)
	{
	if( (E_x.extent(0) == param_ini.Ntheta) && (E_x.extent(1) == param_ini.Nphi)
	        && (E_y.extent(0) == param_ini.Ntheta) && (E_y.extent(1) == param_ini.Nphi)
	        && (E_z.extent(0) == param_ini.Ntheta) && (E_z.extent(1) == param_ini.Nphi) )
		{
		E_inc_x = E_x;
		E_inc_y = E_y;
		E_inc_z = E_z;
		}
	else
		{
		errorstream.str(string());
		errorstream << "Input error: Input matrices should be of size " << param_ini.Ntheta << " x " << param_ini.Nphi << "\n";
		show_error(errorstream.str());
		}
	}
//----------------------
void focus3dws::set_incident_field(complex<double> E_x,
                                   complex<double> E_y,
                                   complex<double> E_z)
	{
	E_inc_x = E_x;
	E_inc_y = E_y;
	E_inc_z = E_z;
	}
//----------------------
void focus3dws::set_scrambler(string filename_real_scrambler, string filename_imag_scrambler)
	{
	//-------------- load real part of G_ttp
	fin.open(filename_real_scrambler.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_real_scrambler << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> real_G_ttp;
	fin >> real_G_ttp;
	fin.close();
	// check the sanity of the data
	if( (real_G_ttp.extent(0) != param_scr.Ntheta) && (real_G_ttp.extent(1) != param_ini.Ntheta)  )
		{
		errorstream.str(string());
		errorstream << "Dimensions of data in file " << filename_real_scrambler << " do not match with " << param_scr.Ntheta << " x " << param_ini.Ntheta << "!" << endl;
		show_error(errorstream.str());
		}

	//-------------- load imaginary part of G_ttp
	fin.open(filename_imag_scrambler.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_imag_scrambler << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> imag_G_ttp;
	fin >> imag_G_ttp;
	fin.close();
	// check the sanity of the data
	if( (imag_G_ttp.extent(0) != param_scr.Ntheta) && (imag_G_ttp.extent(1) != param_ini.Ntheta)  )
		{
		errorstream.str(string());
		errorstream << "Dimensions of data in file " << filename_imag_scrambler << " do not match with " << param_scr.Ntheta << " x " << param_ini.Ntheta << "!" << endl;
		show_error(errorstream.str());
		}

	//---------------- setting G(theta_out,theta_in)
	G_ttp = real_G_ttp + I*imag_G_ttp;
	}
//----------------------
void focus3dws::set_scrambler_from_ascii(string filename_real_scrambler, string filename_imag_scrambler)
	{
	//-------------- load real part of G_ttp
	fin.open(filename_real_scrambler.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_real_scrambler << endl;
		show_error(errorstream.str());
		}
	int num_elements= get_number_of_data_points(filename_real_scrambler);
	if(num_elements != (param_scr.Ntheta * param_ini.Ntheta))
		{
		fin.close();
		errorstream.str(string());
		errorstream << "Number of elements in file " << filename_real_scrambler << " do not match " << param_scr.Ntheta << " x " << param_ini.Ntheta << " = " << (param_scr.Ntheta * param_ini.Ntheta) << "!" << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> real_G_ttp(param_scr.Ntheta,param_ini.Ntheta);
	for(int l = 0; l<param_scr.Ntheta; l++)
		{
		for(int lp =0; lp<param_ini.Ntheta; lp++)
			{
			fin >> real_G_ttp(l,lp);
			}
		}
	fin.close();


	//-------------- load imaginary part of G_ttp
	fin.open(filename_imag_scrambler.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_imag_scrambler << endl;
		show_error(errorstream.str());
		}
	num_elements= get_number_of_data_points(filename_imag_scrambler);
	if(num_elements != (param_scr.Ntheta * param_ini.Ntheta))
		{
		fin.close();
		errorstream.str(string());
		errorstream << "Number of elements in file " << filename_imag_scrambler << " do not match " << param_scr.Ntheta << " x " << param_ini.Ntheta << " = " << (param_scr.Ntheta * param_ini.Ntheta) << "!" << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> imag_G_ttp(param_scr.Ntheta,param_ini.Ntheta);
	for(int l = 0; l<param_scr.Ntheta; l++)
		{
		for(int lp =0; lp<param_ini.Ntheta; lp++)
			{
			fin >> imag_G_ttp(l,lp);
			}
		}
	fin.close();

	//---------------- setting G(theta_out,theta_in)
	G_ttp = real_G_ttp + I*imag_G_ttp;
	}
//----------------------
void focus3dws::set_incident_field(string filename_xcomp, string filename_ycomp, string filename_zcomp)
	{
	// -------------- load xcomp from data file
	fin.open(filename_xcomp.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_xcomp << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> xcomp;
	fin >> xcomp;
	fin.close();
	// check the sanity of the data
	if( (xcomp.extent(0) != param_ini.Ntheta) && (xcomp.extent(1) != param_ini.Nphi)  )
		{
		errorstream.str(string());
		errorstream << "Dimensions of data in file " << filename_xcomp << " do not match!" << endl;
		show_error(errorstream.str());
		}

	// -------------- load ycomp from data file
	fin.open(filename_ycomp.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_ycomp << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> ycomp;
	fin >> ycomp;
	fin.close();
	// check the sanity of the data
	if( (ycomp.extent(0) != param_ini.Ntheta) && (ycomp.extent(1) != param_ini.Nphi)  )
		{
		errorstream.str(string());
		errorstream << "Dimensions of data in file " << filename_ycomp << " do not match!" << endl;
		show_error(errorstream.str());
		}

	// ---------------- load zcomp from data file
	fin.open(filename_zcomp.c_str());
	if(!fin)
		{
		errorstream.str(string());
		errorstream << "Unable to read file " << filename_zcomp << endl;
		show_error(errorstream.str());
		}
	Array<complex<double>,2> zcomp;
	fin >> zcomp;
	fin.close();
	// check the sanity of the data
	if( (zcomp.extent(0) != param_ini.Ntheta) && (zcomp.extent(1) != param_ini.Nphi)  )
		{
		errorstream.str(string());
		errorstream << "Dimensions of data in file " << filename_zcomp << " do not match!" << endl;
		show_error(errorstream.str());
		}

	E_inc_x = xcomp;
	E_inc_y = ycomp;
	E_inc_z = zcomp;
	}

//----------------------
void focus3dws::initialize_scrambler()
	{
	G_ttp.resize(param_scr.Ntheta,param_ini.Ntheta);
	G_ttp = 0.0;
	int diagonal = (param_scr.Ntheta < param_ini.Ntheta) ? param_scr.Ntheta : param_ini.Ntheta;
	//cout << "Number of Diagonal elements = " << diagonal << "\n";
	for(int l =0; l < diagonal;l++)
		{
		G_ttp(l,l) = 1.0;
		}
	}
//----------------------
void focus3dws::set_scrambler(Array<complex<double>,2> G)
	{
	if( (G.extent(0) == param_scr.Ntheta) && (G.extent(1) == param_ini.Ntheta) )
		{
		G_ttp = G;
		}
	else
		{
		errorstream.str(string());
		errorstream << "Input error: Scrambler matrix should be of size " << param_scr.Ntheta << " x " << param_ini.Ntheta << "\n";
		show_error(errorstream.str());
		}
	}
//----------------------
void focus3dws::initialize()
	{
	populate_theta();
	populate_phi();
	populate_x();
	populate_y();
	populate_z();
	initialize_E_incs(); 				//default input polarization is along x-axis
	initialize_scrambler();
	calculate_simpson_13_weights();
	calculate_integrand_exp_args();
	}
//----------------------
void focus3dws::calculate_focus()
	{
	calculate_integrand_amplitudes();

	E_focus_x.resize(param_scr.Nx,param_scr.Ny,param_scr.Nz);
	E_focus_y.resize(param_scr.Nx,param_scr.Ny,param_scr.Nz);
	E_focus_z.resize(param_scr.Nx,param_scr.Ny,param_scr.Nz);
	integrand_exp_term.resize(param_scr.Ntheta,param_scr.Nphi);

	E_focus_x = 0;
	E_focus_y = 0;
	E_focus_z = 0;

	int loop_var = 1;
	for(int i=0; i<param_scr.Nx; i++)
		{
		cout << round(double(loop_var)*100.00/double(param_scr.Nx*param_scr.Ny*param_scr.Nz)) << "%\t";
		cout.flush();
		for(int j=0; j<param_scr.Ny; j++)
			{
			for(int k=0; k<param_scr.Nz; k++)
				{
				integrand_exp_term = exp( + exp_arg1*x(i)
				                          + exp_arg2*y(j)
				                          + exp_arg3*z(k) );
				E_focus_x(i,j,k) = sum( scr_integrand_amplitude_x
				                        * integrand_exp_term
				                        * simpson_13_weights );
				E_focus_y(i,j,k) = sum( scr_integrand_amplitude_y
				                        * integrand_exp_term
				                        * simpson_13_weights );
				E_focus_z(i,j,k) = sum( scr_integrand_amplitude_z
				                        * integrand_exp_term
				                        * simpson_13_weights );
				E_focus_x(i,j,k) = prefactor*E_focus_x(i,j,k);
				E_focus_y(i,j,k) = prefactor*E_focus_y(i,j,k);
				E_focus_z(i,j,k) = prefactor*E_focus_z(i,j,k);
				//cout << E_focus_x(i,j,k) << "\t" << E_focus_y(i,j,k) << "\t" << E_focus_z(i,j,k) << "\n";
				loop_var++;
				}
			}
		}
	cout << "Done\n";
	cout.flush();
	loop_var = 0;
	}
//----------------------
void focus3dws::calculate_focus(Array<complex<double>,3>& x_comp,
                                Array<complex<double>,3>& y_comp,
                                Array<complex<double>,3>& z_comp)
	{
	calculate_focus();
	x_comp = E_focus_x;
	y_comp = E_focus_y;
	z_comp = E_focus_z;
	}
//----------------------
void focus3dws::save_results(string idir)
	{
	namestream.str(string());
	namestream << idir << "_focus_lambda_" << param_ini.lambda << "/";
	string outputdir = namestream.str();
	my_mkdir(outputdir);

	outputstream.str(string());
	outputstream << "\nSaving the results to " << outputdir << "\n";
	show_output(outputstream.str());

	outfilename.str(string());	outfilename << outputdir << "simulation_parameters_initial.dat";
	param_ini.save(outfilename.str());

	outfilename.str(string());	outfilename << outputdir << "simulation_parameters_scrambled.dat";
	param_scr.save(outfilename.str());

	outfilename.str(string());	outfilename << outputdir << "x.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << x;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "y.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << y;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "z.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << z;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "theta_ini.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << theta_ini;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "theta_scr.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << theta_scr;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "phi_ini.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << phi_ini;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "phi_scr.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << phi_scr;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "E_focus_x.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << E_focus_x;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "E_focus_y.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << E_focus_y;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "E_focus_z.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << E_focus_z;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "E_incident_x.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << E_inc_x;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "E_incident_y.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << E_inc_y;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "E_incident_z.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << E_inc_z;
	fout.close();

	outfilename.str(string());	outfilename << outputdir << "scrambler.bar";
	fout.open(outfilename.str().c_str(), ios::out);
	if(!fout)
		{
		show_outfile_error(outfilename.str());
		}
	fout << G_ttp;
	fout.close();

	cout << "Finished!\n";
	}
//------------------
void focus3dws::show_outfile_error(string filename)
	{
	errorstream.str(string());
	errorstream << "Unable to write file " << filename << "\n";
	show_error(errorstream.str());
	}
//-------------------
void focus3dws::show_error(string err_str)
	{
	cerr << "!!!! " << err_str;
	exit(-1);
	}
//-------------------
void focus3dws::show_output(string out_str)
	{
	cout << out_str;
	}

