/*! \file focus3d.h 
  	\author Vishnu Krishnamachari
		\date June 2006 - May 2008
		\brief Defines focus3d class.
		
		Defines a class focus3d. It is used for calculating the vectorial field distribution near the focus of a high numerical aperture objective. This code is based on the mathematical formulation of the problem delineated in "Principles of Nano-optics" by Novotny and Hecht. Determination of the field at at point (x,y,z) in the focal volume consists of integrating all the amplitude contributions from the spherical surface of the lens. If the surface of the lens is described in (f,\f$\theta\f$,\f$\phi\f$) coordinates, where f is the focal length of the lens, the amplitude at (x,y,z) is just a two-dimensional integration over the variables \f$\theta\f$ and \f$\phi\f$ ranging from 0 to \f$\theta_{max}\f$ and 0 to \f$2\pi\f$ respectively. \f$\theta_{max}\f$ corresponds to the angle of acceptance of the lens or the numerical aperture.
*/

#ifndef __FOCUS3D_H__
#define __FOCUS3D_H__

#include "precompile.h"
#include "focus3d_parameters.h"
#include "utilities.h"

		//! A class.
		/*!
				Use this class to calculate the three dimensional distribution of the focal field. This class uses the arrays from the blitz++ library. Hence to use this class and the dipole_radiation 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 three dimensional vectorial focal field distribution and saves the results:
				- \#include "focus3d.h"
				- int main() 
				-	{ ...
					-	focus3d_parameters input_params("parameters.dat"); // parameters.dat is a valid parameters file (see focus3d_parameters) 
					- focus3d focus(input_params);
					- focus.calculate_focus();
					- focus.save_results("directory_prefix");
				-	}
		*/
class focus3d
	{
public:
	//-----------------------------
			//! A Constructor.
			/*! 
			Takes the input parameters through an object of focus3d_parameters class and initializes all the internal matrices. It calls internally the focus3d_parameters::test() memeber function to test the sanity of the parameters. If insane, the program quits. Internally, this constructor calls the reset_parameters() method.
			
			\param q - A focus3d_parameters object
			*/
	focus3d(focus3d_parameters q);
	//-----------------------------
	
	//-----------------------------
	
			//! A Destructor.
			/*! 
			Releases the the memory associated with all the allocated internal matrices.
			*/
	~focus3d();
	//-----------------------------
	
	//-----------------------------
			//! A method for calculating focal field distribution.
			/*! 
			 Calculates the focus field distributions for all the three components of polarization.
			 
			 \return void.
			*/
	void calculate_focus();
	//-----------------------------
	
	//-----------------------------
			//! A method for resetting parameters.
			/*! 
			 Reloads the new parameters and redefines the sizes of the internal matrices.
			 \param q an instance of focus3d_parameters class. Note the effect of using this function is that it resets the effects of previous calls to set_incident_field and set_angular_weights.
			 
			 \param q - A focus3d_parameters object
			 \return void.
			*/
	void reset_parameters(focus3d_parameters q);
	//-----------------------------
	
	//-----------------------------
			//! A method for calculating focal field distributon and outputting it.
			/*! 
			 Calculates the focus field distributions for all the three components saves the results in the variables called by reference. xcomp, ycomp and z_comp are matrices of size (Nx,Ny,Nz) allocated before passing to this function by reference. 
			 
			 Note: if the focal field is desired at just one point (x,y,z) then set the parameter variables Nx = Ny = Nz = 1 and xmin =x, ymin = y and zmin = z.
			 
			 \param[out] x_comp - x-component of the focal field
			 \param[out] y_comp - y-component of the focal field
			 \param[out] z_comp - z-component of the focal field
			*/
	void calculate_focus(Array<complex<double>,3>& x_comp, 
												Array<complex<double>,3>& y_comp, 
												Array<complex<double>,3>& z_comp);
	//-----------------------------
	
	//-----------------------------
			//! 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_focus_x.bar 				: blitz array cooresponding to the x-component of the focal field
			 -# E_focus_y.bar 				: blitz array cooresponding to the y-component of the focal field
			 -# E_focus_z.bar 				: blitz array cooresponding to the z-component of the focal field
			 -# E_incident_x.bar 			: blitz array cooresponding to the x-component of the incident field
			 -# E_incident_y.bar 			: blitz array cooresponding to the y-component of the incident field
			 -# E_incident_z.bar 			: blitz array cooresponding to the z-component of the incident field
			 -# angular_weights.bar 	: blitz array cooresponding to the angular weights / scattering coeffs.
			 -# 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 '_focus_lam_' + num2str(lambda) + '/']
			 - the first five files contain one dimensional real data
			 - the next three files contain three dimensional complex data 
			 - the next four files 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);
	//-----------------------------
	
	//-----------------------------
			//! A method for loading the incident field from arrays.
			/*! 
			 	The incident field is the field that is incident on the flat surface tangent to the lens curvature. 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 (Ntheta, Nphi), i.e. E_x = E_x(\f$\theta\f$,\f$\phi\f$), E_z = E_y(\f$\theta\f$,\f$\phi\f$) and E_z = E_z(\f$\theta\f$,\f$\phi\f$). If the input is just an x-polarized plane parallel beam, then you might set E_x = 1; E_y =0; E_z = 0. This is the default case. 
				
				\param[in] E_x - E_x(\f$\theta\f$,\f$\phi\f$)
				\param[in] E_y - E_y(\f$\theta\f$,\f$\phi\f$)
				\param[in] E_z - E_z(\f$\theta\f$,\f$\phi\f$)
				\return void.		
			*/
	void set_incident_field(Array<complex<double>,2> E_x,
														Array<complex<double>,2> E_y,
														Array<complex<double>,2> E_z);
	//-----------------------------
	
	//-----------------------------
			//! A method for loading the incident field from scalars.
			/*! 
			 	The incident field is the field that is incident on the flat surface tangent to the lens curvature. 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 these three matrices of the size of (Ntheta, Nphi), i.e. E_x = E_x(\f$q\f$,\f$f\f$), E_z = E_y(\f$q\f$,\f$f\f$) and E_z = E_z(\f$q\f$,\f$f\f$). But if the matrices are constant matrices, then this function can be used input the the three complex numbers corresponding to the three polarizations. For example: if the input is just an x-polarized plane parallel beam, then you might set E_x = 1; E_y =0; E_z = 0. This is the default case. 
				
				\param[in] E_x - value of x-polarization of the field
				\param[in] E_y - value of y-polarization of the field
				\param[in] E_z - value of z-polarization of the field
				\return void.		
			*/
	void set_incident_field(complex<double> E_x,
                                  complex<double> E_y,
                                  complex<double> E_z);
	//-----------------------------
	
	//-----------------------------
			//! A method for loading the incident field from files.
			/*! 
			 	The incident field is the field that is incident on the flat surface tangent to the lens curvature. It has three components corresponding to x, y and z polarizations. 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] filename_xcomp - name of the file containing x-component of the incident field
				\param[in] filename_ycomp - name of the file containing y-component of the incident field
				\param[in] filename_zcomp - name of the file containing z-component of the incident field	
				\return void.		
			*/
	void set_incident_field(string filename_xcomp, string filename_ycomp, string filename_zcomp);
	//-----------------------------
	
	//-----------------------------
			//! A method for loading the angular weights from an array.
			/*!
				This function is used to include the effects of scattering on the focal amplitude. It populates an internal variable called angular_weights which effectively describes how much the rays traveling at angle (\f$\theta\f$,\f$\phi\f$) gets affected by the medium. If the medium is non-scattering then the input should be a uniform matrix of value 1.0 which is the case by default.		  
				\param[in] aweights - matrix of size (Ntheta,Nphi) which effectively describes the scattering function of the medium.
				\return void.		
				
			*/
	void set_angular_weights(Array<complex<double>,2> aweights);
	//-----------------------------
	
	//-----------------------------
			//! A method for loading the angular weights from files.
			/*!
				This function is used to include the effects of scattering on the focal amplitude. It populates an internal variable called angular_weights which effectively describes how much the rays traveling at angle (\f$\theta\f$,\f$\phi\f$) gets affected by the medium.
				
				\param[in] ampl_coeff_file - a std::string containing the name of the file with polynomial coefficients that fits the amplitude of the scattering of the medium.
				\param[in] phase_coeff_file - a std::string containing the name of the file with polynomial coefficients that fits the phase of the scattering of the medium.
				\return void.		
				
			*/
	void set_angular_weights(string ampl_coeff_file, string phase_coeff_file);
	
	//-----------------------------
			//! A method for loading the angular weights from a file.
			/*!
				This function is used to include the effects of scattering on the focal amplitude. It populates an internal variable called angular_weights which effectively describes how much the rays traveling at angle (\f$\theta\f$,\f$\phi\f$) gets affected by the medium.
				
				\param[in] g_of_theta_file - std::string with filename of a file containing three colums corresponding to \f$ \theta \f$, \f$ |G(\theta)| \f$, and angle(\f$ G(\theta) \f$).
				\return void.		
				
			*/
	void set_angular_weights(string g_of_theta_file);
	//-----------------------------
	
private:
	
	focus3d_parameters param;
	
	double Dx, Dy, Dz, Dtheta, Dphi;

	complex<double> I;		/*!< imaginary number i*/
	complex<double> prefactor;
	
	Array<double,1> x, y, z, theta, phi;
	Array<double,2> simpson_13_weights;
	Array<double,2> theta2;
	Array<double,2> phi2;
	Array<complex<double>,2> E_inc_x, E_inc_y, E_inc_z;
	Array<complex<double>,2> integrand_amplitude_x;
	Array<complex<double>,2> integrand_amplitude_y;
	Array<complex<double>,2> integrand_amplitude_z;
	Array<complex<double>,2> angular_weights; 				// G(k) function
	Array<complex<double>,2> exp_arg1;
	Array<complex<double>,2> exp_arg2;
	Array<complex<double>,2> exp_arg3;
	Array<complex<double>,2> integrand_exp_term;
	Array<complex<double>,3> E_focus_x, E_focus_y, E_focus_z;

	stringstream outputstream, errorstream, namestream, outfilename;;
	ofstream fout;
	ifstream fin;

	bool test_parameters(focus3d_parameters p);
	void load_parameters(focus3d_parameters p);
	void calculate_simpson_13_weights();
	void initialize_E_incs();
	void initialize_angular_weights();
	void calculate_integrand_amplitudes();
	void calculate_integrand_exp_args();
	void initialize();
	void populate_x();
	void populate_y();
	void populate_z();
	void populate_theta();
	void populate_phi();
	void show_outfile_error(string filename);
	void show_output(string output_str);
	void show_error(string error_str);
	};
//---------------------------



#endif
