/*! \file dipole_radiation.h
		\author Vishnu Krishnamachari
		\date June 2006 - May 2008
		\brief Defines dipole_radiation class.
	
		Defines a class dipole_radiation. It is used for calculating the vectorial far-field distribution due to the radiation from dipole distribution near the focus of an objective. The field distribution at a far-field point R = (X,Y,Z) from a dipole 'p' located at r = (x,y,z) is -[1/4.pi.R^3].exp(i.k.f).( (R-r) x (R-r) x p). Note that R and r are coordinate vectors from the origin joining the corresponding points. If dipoles are distributed through out a region of [xmin:xmax],[ymin:ymax] and [zmin:zmax] then the far-field amplitude at R would be a three-dimensional integration of amplitudes from all the individual points (x,y,z). The 3d integration is implemented ucing on Simpson 1/3 rule. If the surface of the collecting lens is described in (f,\f$\theta\f$,\f$\phi\f$) coordinates, where f is the focal length of the lens, the above integration would provide the amplitude on this surface. This class also provides a method to calculate the radiation pattern after refraction at the spherical surface of the lens. 
*/

#ifndef __DIPOLE_RADIATION_H__
#define __DIPOLE_RADIATION_H__

#include "focus3d.h"

	/*! \def dipole_parameters
			dipole_parameters is an alias to focus3d_parameters class
	*/
#define dipole_parameters focus3d_parameters


	//! A class.
	/*!
			Create an instance of this class, pass the simulation parameters and the source dipole strengths to calculate the far field radiation on the sphereical surface of the collecting lens. A method is also provided to calculate the radiation after the refraction at the collecting lens' spherical surface. This class uses the arrays from the blitz++ library. Hence to use this class and the focus3d class it would be a good idea to read up the documentation of blitz++ arrays.
		
			Using this class in your code is pretty simple. The following piece of code generates the vectorial far field due to a dipole distribution with all the dipoles being oriented along x-direction:
			- \#include "dipole_radiation.h"
			- int main() 
			-	{
					-	dipole_parameters param("parameters.dat"); // parameters.dat is a valid parameters file (see focus3d_parameters) 
					- // define two-dimensional matrices to receive the calculated radiation
					- Array<complex<double>,2> E_rad_x(param.Ntheta,param.Nphi);
					- Array<complex<double>,2> E_rad_y(param.Ntheta,param.Nphi);
					- Array<complex<double>,2> E_rad_z(param.Ntheta,param.Nphi);					
					- dipole_radiation drad(input_params);
					- drad.calculate_refracted_radiation(E_rad_x, E_rad_y, E_rad_z);
					- drad.save_results("directory_prefix");
			-	}
		
	*/
class dipole_radiation
{
	public:
		//------------------------------------
		//! A Constructor.
			/*! 
			Takes the input parameters through an object of dipole_parameters class which is an alias to focus3d_parameters class and initializes all the internal matrices. It calls internally the dipole_parameters::test() member function to test the sanity of the parameters. If insane, the program quits. Internally, this constructor calls the reset_parameters() method.
			
			\param p - A focus3d_parameters object
			*/
		dipole_radiation(dipole_parameters p);
		//------------------------------------
		
		//------------------------------------
		//! A Destructor.
			/*! 
			Releases the the memory associated with all the allocated internal matrices.
			*/
		~dipole_radiation();
		//------------------------------------
		
		//------------------------------------
		//! A method for resetting parameters.
			/*! 
			 Reloads the new parameters and redefines the sizes of the internal matrices.
			 \param p an instance of focus3d_parameters class. Note the effect of using this function is that it resets the effects of previous calls to set_dipole_field.
			 \return void.
			*/
		void reset_parameters(dipole_parameters p);
		//------------------------------------
		
		//------------------------------------
		//! A method for loading the input dipole field from arrays.
			/*! 
			 	The dipole field is the field that is at the focus of the focusing objective. It has three components corresponding to x, y and z polarizations. If we denote them as E_x, E_y and E_z, then the inputs should be three matrices of the size of (Nx, Ny, Nz), i.e. E_x = E_x(x,y,z), E_z = E_y(x,y,z) and E_z = E_z(x,y,z). The default case is E_x = 1; E_y =0; E_z = 0. 
				
				\param[in] E_x - E_x(x,y,z)
				\param[in] E_y - E_y(x,y,z)
				\param[in] E_z - E_z(x,y,z)
				\return void.		
			*/
		void set_dipole_field(Array<complex<double>,3> E_x,
								Array<complex<double>,3> E_y,
								Array<complex<double>,3> E_z);
		//------------------------------------
		
		//------------------------------------
		//! A method for loading the input dipole field from files.
			/*! 
			 	The dipole field is the field that is at the focus of the focusing objective. It has three components corresponding to x, y and z polarizations. If we denote them as E_x, E_y and E_z, then the inputs should be three matrices of the size of (Nx, Ny, Nz), i.e. E_x = E_x(x,y,z), E_z = E_y(x,y,z) and E_z = E_z(x,y,z). (The default case is E_x = 1; E_y =0; E_z = 0.) If these three components are stored in three different files in blitz array format in column major ordering then, this function can be used to load them. 
				
				\param[in] file_xcomp - name of the file containing x-component of the dipole field
				\param[in] file_ycomp - name of the file containing y-component of the dipole field
				\param[in] file_zcomp - name of the file containing z-component of the dipole field
				\return void.		
			*/
		void set_dipole_field(string file_xcomp, string file_ycomp, string file_zcomp);
		//------------------------------------
		
		//------------------------------------
		//! A method for calculating far-field radiation on the spherical surface of the collecting lens.
			/*! 
			 Calculates the far-field distributions for all the three components and saves the results in the variables called by reference. xcomp, ycomp and z_comp are matrices of size (Ntheta,Nphi) allocated before passing to this function by reference. 
			 
			 Note: if the far-field is desired at just one point (\f$\theta\f$,\f$\phi\f$) then set the parameter variables Ntheta = Nphi = 1 and theta_min = \f$\theta\f$, and phi_min = \f$\phi\f$.
			 
			 \param[out] E_rad_x - x-component of the far-field
			 \param[out] E_rad_y - y-component of the far-field
			 \param[out] E_rad_z - z-component of the far-field
			 \return void.
			*/
		void calculate_radiation(Array<complex<double>, 2>& E_rad_x,
								Array<complex<double>, 2>& E_rad_y,
								Array<complex<double>, 2>& E_rad_z);
		//------------------------------------
		
		//------------------------------------
		//! A method for calculating far-field radiation after refraction at the spherical surface of the collecting lens.
			/*! 
			 Calculates the far-field distributions after refraction for all the three components and saves the results in the variables called by reference. E_refr_x, E_refr_y and E_refr_z are matrices of size (Ntheta,Nphi) allocated before passing to this function by reference. 
			 
			 Note: 
			 -# if the far-field is desired at just one point (\f$\theta\f$,\f$\phi\f$) then set the parameter variables Ntheta = Nphi = 1 and theta_min = \f$\theta\f$, and phi_min = \f$\phi\f$.
			 -# call this method only after invoking calculate_radiation(). Internally, this function just performs a coordinate transformation operations on the matrices corresponding to the radiation on the spherical surface of the collecting lens. Hence if this function is called before invoking the calculate_radiation() method, the results would be spurious.   
			 
			 \param[out] E_refr_x - x-component of the refracted far-field
			 \param[out] E_refr_y - y-component of the refracted far-field
			 \param[out] E_refr_z - z-component of the refracted far-field
			 \return void.
			*/
		void calculate_refraction(Array<complex<double>, 2>& E_refr_x,
								Array<complex<double>, 2>& E_refr_y,
								Array<complex<double>, 2>& E_refr_z);
		//------------------------------------
		
		//------------------------------------
		//! A method for calculating far-field radiation after refraction at the spherical surface of the collecting lens.
			/*! 
			 Calculates the far-field distributions after refraction for all the three components and saves the results in the variables called by reference. E_refr_x, E_refr_y and E_refr_z are matrices of size (Ntheta,Nphi) allocated before passing to this function by reference. 
			 
			 Note: 
			 -# if the far-field is desired at just one point (\f$\theta\f$,\f$\phi\f$) then set the parameter variables Ntheta = Nphi = 1 and theta_min = \f$\theta\f$, and phi_min = \f$\phi\f$.
			 -# Internally, this function invokes calculate_radiation() and calculate_refraction() in series in that order.    
			 
			 \param[out] E_rad_x - x-component of the refracted far-field
			 \param[out] E_rad_y - y-component of the refracted far-field
			 \param[out] E_rad_z - z-component of the refracted far-field
			 \return void.
			*/
		void calculate_refracted_radiation(Array<complex<double>, 2>& E_rad_x,
											Array<complex<double>, 2>& E_rad_y,
											Array<complex<double>, 2>& E_rad_z);
		//------------------------------------
		
		//------------------------------------
		//! A method for saving the calculated results in a file.
			/*! 
			 saves the calculated results to different files. The generated output files are:
			 -# x.bar 								: blitz array corresponding to x coordinates
			 -# y.bar 								: blitz array corresponding to y coordinates
			 -# z.bar 								: blitz array corresponding to z coordinates
			 -# theta.bar 						: blitz array corresponding to \f$\theta\f$ coordinates
			 -# phi.bar								: blitz array corresponding to \f$\phi\f$ coordinates
			 -# E_dipole_x.bar 				: blitz array cooresponding to the x-component of the dipole field
			 -# E_dipole_y.bar 				: blitz array cooresponding to the y-component of the dipole field
			 -# E_dipole_z.bar 				: blitz array cooresponding to the z-component of the dipole field
			 -# E_far_x.bar 					: blitz array cooresponding to the x-component of the far-field
			 -# E_far_y.bar 					: blitz array cooresponding to the y-component of the far-field
			 -# E_far_z.bar 					: blitz array cooresponding to the z-component of the far-field
			 -# E_refracted_x.bar 		: blitz array cooresponding to the x-component of the refracted far-field
			 -# E_refracted_y.bar 		: blitz array cooresponding to the y-component of the refracted far-field
			 -# E_refracted_z.bar 		: blitz array cooresponding to the z-component of the refracted far-field
			 -# simulation_parameters.dat : an ascii file with the applied simulation parameters
			 
			 Note:
			 - all the files are stored in a subdirectory called (in MATLAB format) [idir '_rad_lam_' + num2str(lambda) + '/']
			 - the first five files contain one dimensional real data
			 - the next three files contain three dimensional complex data as a function of (x,y,z)
			 - the next sixfiles contain two dimensional complex data as a function of (\f$\theta\f$, \f$\phi\f$) 
			 - each complex number A+iB is represented as (A,B)
			 - the class uses the blitz++ arrays in column major format and hence the matrix A(i,j,k) is stored with the index k as the fastest running one and the index i being the slowest running one.
			 
			 \param idir - prefix to the directory to which the results have to be saved 
			 \return void.			 
			*/
		void save_results(string idir);
		//------------------------------------
		
	private:
		dipole_parameters param;

		double 	Dx, Dy, Dz, Dtheta, Dphi;
		stringstream outputstream, errorstream, namestream, outfilename;;
		ofstream fout;
		ifstream fin;
		
		complex<double> I;
		
		Array<double, 1> x, y, z, theta, phi;
		Array<double, 2> theta2, phi2;
		Array<double, 3> x3, y3, z3;
		Array<double, 3> simpson_13_weights;

		Array<complex<double>, 2> E_far_x, E_far_y, E_far_z;
		Array<complex<double>, 2> E_refracted_x, E_refracted_y, E_refracted_z; 
		Array<complex<double>, 3> E_dip_x, E_dip_y, E_dip_z;
		
		void load_parameters(dipole_parameters d);
		bool test_parameters(dipole_parameters d);
		void populate_x();
		void populate_y();
		void populate_z();
		void populate_theta();
		void populate_phi();
		void populate_far_field();
		void populate_refracted_field();
		void calculate_simpson_13_weights();
		void initialize();
		void initialize_E_dipole();
		void show_outfile_error(string filename);
		void show_output(string output_str);
		void show_error(string error_str);		
		void calculate_radiation();
		void calculate_refraction();		
};

#endif
