#ifndef _SUPERDIRECTIVITY_H
#define _SUPERDIRECTIVITY_H

#include <complex>
#include <vector>
#include <iostream>

#include "sph_bessel.h"

#define _USE_MATH_DEFINES
#include <cmath>

Real get_directivity_sphml_dZx_epsRe(std::vector<Real> &param, void *cls);
Real get_directivity_sphml_dZx_epsRe_parallel(std::vector<Real> &param, void *cls);

const int _degree = 30;

class sphere_ml {
public:
	sphere_ml() {
		deg_max = _degree;
		ampl_in = new Complex [4*deg_max];
		ampl_out = new Complex [4*deg_max];
		rt = new Complex [8*deg_max];
		rtml_in = new Complex [8*deg_max];
		rtml_out = new Complex [8*deg_max];
		pi1 = new Real [deg_max];
		tau0 = new Real [deg_max];
		tau1 = new Real [deg_max];
		theta = phi = Real(0);
		set_directivity_angles(theta, phi);
	}
	sphere_ml(int deg_max_) {
		deg_max = deg_max_;
		ampl_in = new Complex [4*deg_max];
		ampl_out = new Complex [4*deg_max];
		rt = new Complex [8*deg_max];
		rtml_in = new Complex [8*deg_max];
		rtml_out = new Complex [8*deg_max];
		pi1 = new Real [deg_max];
		tau0 = new Real [deg_max];
		tau1 = new Real [deg_max];
		theta = phi = Real(0);
		set_directivity_angles(theta, phi);
	}
	sphere_ml(const sphere_ml &sphml) {
		delete [] ampl_in; delete [] ampl_out;
		delete [] rt; delete [] rtml_in; delete [] rtml_out;
		delete [] pi1; delete [] tau0; delete [] tau1;
		deg_max = sphml.deg_max;
		ampl_in = new Complex [4*deg_max];
		ampl_out = new Complex [4*deg_max];
		rt = new Complex [8*deg_max];
		rtml_in = new Complex [8*deg_max];
		rtml_out = new Complex [8*deg_max];
		pi1 = new Real [deg_max];
		tau0 = new Real [deg_max];
		tau1 = new Real [deg_max];
		theta = sphml.theta; phi = sphml.phi;
		set_directivity_angles(theta, phi);
	}
	~sphere_ml() {
		delete [] ampl_in; delete [] ampl_out;
		delete [] rt; delete [] rtml_in; delete [] rtml_out;
		delete [] pi1; delete [] tau0; delete [] tau1;
	}

	int get_degree() const {return deg_max;}

		// pre-define direction in which the directivity should be evaluated (0,0 by default)
	void set_directivity_angles(Real theta, Real phi);
		// return directivity for an electric dipole placed at -kzd on axis Z
		// in the presence of a spherical multilayer: kR - vector of radii of
		// spherical interfaces; eps - vector of permittivities (should be one element larger than kR)
	Real get_directivity_edipole_x(Real kzd, const std::vector<Real> &kR, const std::vector<Complex> &eps);

		// return the amplitudes of the emitted waves, which were calculated at the last directivity evaluation
	std::vector<Complex> get_last_ampl() const;
		// calculate the radiation pattern
	std::vector<Complex> ampl_far(Real th_rad, Real ph_rad);

private:
	int deg_max; // maximum degree of spherical hamonics
	Real theta, phi; // directivity angles
	Real *pi1, *tau0, *tau1; // angular function arrays
	Complex *ampl_in, *ampl_out; // amplitude arrays
	Complex *rt, *rtml_in, *rtml_out; // s-matrix coefficient arrays

		// electric dipole is located at (-kz) and has amplitude 1 along x:
	void get_edipole_amp_x(Complex kz, Complex *ampl, bool in);
		// elextric dipole is located at (-kz) and has amplitude 1 along z:
	void get_edipole_amp_z(Complex kz, Complex *ampl, bool in);
		// elextric dipole is located at (-kz) and has amplitude 1 directed at angle th relative to z:
	void get_edipole_amp_xz(Complex kz, Complex *ampl, Real th, bool in);

		// reflection and transmission coefficients at spherical interface
	void sphere_sm(Real kr, Complex eps_in, Complex eps_out, Complex *res);
		// T-matrix of a spherical interface
	void sphere_tm(Real kr, Complex eps_in, Complex eps_out, Complex *res);
		// composition of S-matrices
	void compose_sm(Complex *in, Complex *out, Complex *res);
		// composition of T-matrices
	void compose_tm(Complex *in, Complex *out, Complex *res);

		// directivity of a dipole located at axiz Z
	Real directivity_dipole_z(Complex *ampl) const;
};

//#################################################################################
// class for parallel calculation:

class sphere_ml_parallel {
public:
	sphere_ml_parallel() {
		deg_max = _degree;
		pi1 = new Real [_degree];
		tau0 = new Real [_degree];
		tau1 = new Real [_degree];
		theta = phi = Real(0);
		set_directivity_angles(theta, phi);
	}

	int get_degree() const {return deg_max;}

		// pre-define direction in which the directivity should be evaluated (0,0 by default)
	void set_directivity_angles(Real theta, Real phi);
		// return directivity for an electric dipole placed at -kzd on axis Z
		// in the presence of a spherical multilayer: kR - vector of radii of
		// spherical interfaces; eps - vector of permittivities (should be one element larger than kR)
	Real get_directivity_edipole_x(Real kzd, const std::vector<Real> &kR, const std::vector<Complex> &eps);

private:
	int deg_max; // maximum degree of spherical hamonics
	Real theta, phi; // directivity angles
	Real *pi1, *tau0, *tau1; // angular function arrays

		// electric dipole is located at (-kz) and has amplitude 1 along x:
	void get_edipole_amp_x(Complex kz, Complex *ampl, bool in);
		// elextric dipole is located at (-kz) and has amplitude 1 along z:

		// reflection and transmission coefficients at spherical interface
	void sphere_sm(Real kr, Complex eps_in, Complex eps_out, Complex *res);
		// composition of S-matrices
	void compose_sm(Complex *in, Complex *out, Complex *res);

		// directivity of a dipole located at axiz Z
	Real directivity_dipole_z(Complex *ampl) const;
};

#endif // _SUPERDIRECTIVITY_H