#include <iostream>
#include <fstream>
#include <limits>
#include <mpi.h>

#include "superdirectivity.h"
#include "sph_bessel.h"
#include "jade.h"

using namespace std;

int main(int argc, char* argv[]) {
		// directivity calculation
/**/
	{
		cout.precision(15);

		int ni = 3, ND = 30;
		Real kz, wl, wv, dir;
		vector<Real> kR;
		vector<Complex> eps;

		wl = 0.455; wv = 2*_pi/wl;

		sphere_ml sph(ND);

//		kR.push_back(2*_pi/0.455*0.07985);
//		eps.push_back(pow(Real(6.325),2));
//		eps.push_back(Real(1));
/**
		kR.push_back(2*_pi*0.0486056995345798/0.455);
		kR.push_back(2*_pi*0.0824637171512022/0.455);
		kR.push_back(2*_pi*0.08645/0.455);
		eps.push_back(pow(Real(4.98055106546034),2));
		eps.push_back(pow(Real(5.0),2));
		eps.push_back(pow(Real(3.34491667760117),2));
		eps.push_back(Real(1));
/**/
		kR.push_back(0.023666594336234 * wv);
		kR.push_back(0.136411859930056 * wv);
		kR.push_back(0.1365 * wv);
		eps.push_back(pow(Real(1.0),2));
		eps.push_back(pow(Real(7.97246729336805),2));
		eps.push_back(pow(Real(5.17006034408423),2));
		eps.push_back(Real(1)); // premittivity of the supprounding medium

			// dependence of the directivity from the dipole position
		ofstream out_dir("e:/Alexey/03 - code/data/directivity/dir.dat");
		for (kz=0.01*kR[2]; kz<1.5*kR[2]; kz+=.01*kR[2]) {
			out_dir << kz << " " << sph.get_directivity_edipole_x(kz,kR,eps) << endl;
		}
			// single directivity evaluation
		kz = 2*_pi*0.299; //0.5 * (kR[0] + kR[1]);
		dir = sph.get_directivity_edipole_x(kz,kR,eps);
		cout << "directivity at position kz = " << kz << " equals to " << dir << endl;
			
			// get amplitude vector from the last directivity calculation:
		vector<Complex> v_amp = sph.get_last_ampl();
		cout << "amplitudes of e-waves (m = 1):\n";
		for (int n=0; n<ND; ++n) cout << n+1 << " " << (v_amp[3*n]) << endl;
		cout << "amplitudes of h-waves (m = 1):\n";
		for (int n=0; n<ND; ++n) cout << n+1 << " " << (v_amp[3*n+1]) << endl;
		cout << "amplitudes of h-waves (m = 0):\n";
		for (int n=0; n<ND; ++n) cout << n+1 << " " << (v_amp[3*n+2]) << endl;
			
			// radiation pattern based on the last directivity calculation
		vector<Complex> E_par, E_perp;
		ofstream out_pat("e:/Alexey/03 - code/data/directivity/pat.dat");
		for (Real th=0; th<361; th+=1) {
			E_par = sph.ampl_far(th*_pi/180, Real(0));
			E_perp = sph.ampl_far(th*_pi/180, 0.5*_pi);
			out_pat << th+90 << " " << (norm(E_par[0])) << " " << (norm(E_perp[1])) << endl;
		}

		return 0;
	}
/**/
/**/
/**/
		// directivity optimization section
/**/
	{
		int ni = 3; // number of spherical interfaces
		Real kr_max = Real(1); // maximum kr
		Real eps_max = Real(225); // maximum epsilon
		sphere_ml_parallel cls_sph; // class for calculation of directivity
		vector<Real> vb_min, vb_max; // parameter boundaries

			// boundaries for all optimization parameters:
		vb_min.push_back(Real(0.01)); vb_max.push_back(Real(1.5)*kr_max); // dipole position
		for (int i=0; i<ni; ++i)
			vb_min.push_back(Real(0.01)); vb_max.push_back(Real(1.5)*kr_max); // radii
		for (int i=0; i<ni; ++i)
			vb_min.push_back(Real(1)); vb_max.push_back(eps_max); // real permittivities

		ofstream out_opt("e:/Alexey/03 - code/data/directivity/opt.dat");

		MPI_Init(&argc, &argv);

		int rank;
		MPI_Comm_rank(MPI_COMM_WORLD, &rank);

			// Settings for optimization algorithm;
		int dimension = 2*ni+1; /// Number of parameters to optimize.
		int total_population = 3*dimension;  /// Total number of individuals in population.
		jade::SubPopulation sub_population;
		sub_population.Init(total_population, dimension);
		sub_population.FitnessFunction_p = &get_directivity_sphml_dZx_epsRe_parallel;
			/// Low and upper bound for all dimensions;
			//sub_population.SetAllBoundsVectors({0.005, 0.005}, {0.05, 0.02});
		sub_population.SetAllBoundsVectors(vb_min, vb_max);
		//sub_population.SetTargetToMinimum();
		sub_population.SetTargetToMaximum();
		sub_population.SetTotalGenerationsMax(10000);
		sub_population.SetBestShareP(0.05);
		sub_population.SetAdapitonFrequencyC(0.1);
		sub_population.SetDistributionLevel(0);
		//sub_population.PrintParameters("f1");
		//sub_population.RunOptimization(out_opt, reinterpret_cast<void*>(&cls_sph));
		sub_population.RunOptimization(cout, reinterpret_cast<void*>(&cls_sph));
		sub_population.PrintResult("directivity function ");
		MPI_Finalize();

		return 0;
	}
/**/
}