#include "successive_focus_radiate_focus.h"

//! For calculating the first focus, the far-field radiation, refraction and the second focus.
int main(int argc, char** argv)
	{
	if(argc < 3)
		{
		cerr << "Usage: " << argv[0] << " <parameters file> <output directory>\n";
		exit(-1);
		}
	string param_file = argv[1];
	string outputdir = argv[2];

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

	Array<complex<double>,2> E_inc_x, E_inc_y, E_inc_z;
	Array<complex<double>,3> E_focus1_x, E_focus1_y, E_focus1_z;
	Array<complex<double>,2> E_rad_x, E_rad_y, E_rad_z;
	Array<complex<double>,3> E_focus2_x, E_focus2_y, E_focus2_z;

	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 = 1.0;
	E_inc_y = 0.0;
	E_inc_z = 0.0;

	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_rad_x.resize(p.Ntheta,p.Nphi); E_rad_x = 0;
	E_rad_y.resize(p.Ntheta,p.Nphi); E_rad_y = 0;
	E_rad_z.resize(p.Ntheta,p.Nphi); E_rad_z = 0;

	E_focus2_x.resize(p.Nx,p.Ny,p.Nz); E_focus2_x = 0;
	E_focus2_y.resize(p.Nx,p.Ny,p.Nz); E_focus2_y = 0;
	E_focus2_z.resize(p.Nx,p.Ny,p.Nz); E_focus2_z = 0;

	cout << "---------------------------------------------\n";
	cout << "---------------- first focus ----------------\n";
	cout << "---------------------------------------------\n";

	focus3d focus1(p);
	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(string(outputdir + "/coarse_focus1"));

	cout << "---------------------------------------------\n";
	cout << "------------------ radiation ----------------\n";
	cout << "---------------------------------------------\n";

	dipole_radiation drad(p);
	drad.set_dipole_field(E_focus1_x,E_focus1_y,E_focus1_z);
	drad.calculate_refracted_radiation(E_rad_x,E_rad_y,E_rad_z);
	drad.save_results(string(outputdir + "/coarse_radiation"));

	cout << "---------------------------------------------\n";
	cout << "---------------- second focus ---------------\n";
	cout << "---------------------------------------------\n";

	focus3d focus2(p);
	focus2.set_incident_field(E_rad_x,E_rad_y,E_rad_z);
	focus2.calculate_focus(E_focus2_x,E_focus2_y,E_focus2_z);
	focus2.save_results(string(outputdir+"/coarse_focus2"));

	return(0);
	}
