#include "focal_volume_HG01.h"

int main(int argc, char** argv)
	{
	if(argc < 5)
		{
		cerr << "Usage: " << argv[0] << " <parameters file> <amplitude polynomial coefficients file> <phase polynomial coefficients file> <output directory>\n";
		exit(-1);
		}
	string param_file = argv[1];
	string amp_coeff_file = argv[2];
	string phase_coeff_file = argv[3];
	string outputdir = argv[4];

	focus3d_parameters p;
	if(!p.load(param_file)) {exit (-1);} // load the parameters file
	p.print();

	Array<complex<double>,3> E_focus1_x, E_focus1_y, E_focus1_z;
	Array<complex<double>,2> E_inc_x, E_inc_y, E_inc_z;

	E_focus1_x.resize(p.Nx,p.Ny,p.Nz); E_focus1_x = 0;
	E_focus1_y.resize(p.Nx,p.Ny,p.Nz); E_focus1_y = 0;
	E_focus1_z.resize(p.Nx,p.Ny,p.Nz); E_focus1_z = 0;

	E_inc_x.resize(p.Ntheta,p.Nphi);
	E_inc_y.resize(p.Ntheta,p.Nphi);
	E_inc_z.resize(p.Ntheta,p.Nphi);
	E_inc_x = 0.0;
	E_inc_y = 0.0;
	E_inc_z = 0.0;

	double Dphi;
	if(p.Nphi == 1) {		Dphi = 0; }
	else if(p.Nphi > 1) {		Dphi = (p.phi_max-p.phi_min)/double(p.Nphi-1); }
	else	{	cerr << "Nphi = " << p.Nphi << " seems to set false! Exiting!\n";	exit(-1);	}

	Array<double,1> phi(p.Nphi);
	phi = p.phi_min + tensor::i*Dphi;

	for(int l =0; l<p.Ntheta; l++)
		{
		for(int m=0; m<p.Nphi; m++)
			{
			if( (phi(m) >=0 ) && (phi(m) < M_PI) )
				{
				E_inc_x(l,m) = +1.0;
				}
			else
				{
				E_inc_x(l,m) = -1.0;
				}
			}
		}

	focus3d focus1(p);
	//focus1.set_angular_weights(amp_coeff_file,phase_coeff_file);
	focus1.set_incident_field(E_inc_x,E_inc_y,E_inc_z);
	focus1.calculate_focus(E_focus1_x,E_focus1_y,E_focus1_z);
	focus1.save_results(outputdir);

	return(0);
	}




