Добавьте файлы проекта.

jade.cpp

@@ -0,0 +1,834 @@
+///
+/// @file   jade.cc
+/// @author Ladutenko Konstantin <kostyfisik at gmail (.) com>
+/// @date   Thu Aug 15 19:26:53 2013
+/// @copyright 2013 Ladutenko Konstantin
+/// @section LICENSE
+/// This file is part of JADE++.
+///
+/// JADE++ is free software: you can redistribute it and/or modify
+/// it under the terms of the GNU General Public License as published by
+/// the Free Software Foundation, either version 3 of the License, or
+/// (at your option) any later version.
+///
+/// JADE++ is distributed in the hope that it will be useful,
+/// but WITHOUT ANY WARRANTY; without even the implied warranty of
+/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+/// GNU General Public License for more details.
+///
+/// You should have received a copy of the GNU General Public License
+/// along with JADE++.  If not, see <http://www.gnu.org/licenses/>.
+///
+/// @brief JADE++ is a free (GPLv3+) high performance implementation of
+/// adaptive differential evolution optimization algorithm from
+/// Jingqiao Zhang and Arthur C. Sanderson book 'Adaptive Differential
+/// Evolution. A Robust Approach to Multimodal Problem Optimization'
+/// Springer, 2009.  Crossover rate was patched according to PMCRADE
+/// approach supposed by Jie Li, Wujie Zhu, Mengjun Zhou, and Hua Wang
+/// in 'Power Mean Based Crossover Rate Adaptive Differential
+/// Evolution' in H. Deng et al. (Eds.): AICI 2011, Part II, LNAI
+/// 7003, pp. 34–41, 2011
+#include "./jade.h"
+#include <mpi.h>
+#include <random>
+#include <cstdio>
+#include <cmath>
+#include <algorithm>
+#include <map>
+#include <iterator>
+#include <string>
+#include <vector>
+#include <iostream>
+
+namespace jade {
+  /// @todo Replace all simple kError returns with something meangfull. TODO change kError to throw exception
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::Selection(std::vector<Real> crossover_u, long i)  {
+    bool is_evaluated = false;
+    Real f_current = 0.0;
+    for (auto f : evaluated_fitness_for_current_vectors_) {
+      if (f.second == i) {
+        f_current = f.first;
+        is_evaluated = true;
+      }
+    }  // end of searching of pre-evaluated fitness for current individual
+    if (!is_evaluated) error_status_ = kError;
+    Real f_best = evaluated_fitness_for_current_vectors_.front().first;
+    Real f_crossover_u = FitnessFunction(crossover_u);
+    bool is_success = f_crossover_u > f_current
+        || f_crossover_u == f_best;  //Selected for maxima search
+    if (is_find_minimum_) is_success = !is_success;
+    if (!is_success) {
+      // Case of current x and f were new for current generation.
+      x_vectors_next_generation_[i] = x_vectors_current_[i];
+      for (auto &f : evaluated_fitness_for_next_generation_) {
+        if (f.second == i) f.first = f_current;
+      }  // end of saving old fitness value in new generation.
+    } else {  // if is_success == true
+      x_vectors_next_generation_[i] = crossover_u;
+      for (auto &f : evaluated_fitness_for_next_generation_) 
+        if (f.second == i) f.first = f_crossover_u;
+      to_be_archived_best_A_.push_back(x_vectors_current_[i]);
+      successful_mutation_parameters_S_F_.push_back(mutation_F_[i]);
+      successful_crossover_parameters_S_CR_.push_back(crossover_CR_[i]);
+      // if (process_rank_ == kOutput)
+      //   printf("n%li f_new=%4.2f\n",i,f_crossover_u);
+      //PrintSingleVector(crossover_u);
+    }  // end of dealing with success crossover
+    return kDone;
+  } // end of int SubPopulation::Selection(std::vector<Real> crossover_u);
+
+  int SubPopulation::Selection(std::vector<Real> crossover_u, void* param, long i)  {
+    bool is_evaluated = false;
+    Real f_current = 0.0;
+    for (auto f : evaluated_fitness_for_current_vectors_) {
+      if (f.second == i) {
+        f_current = f.first;
+        is_evaluated = true;
+      }
+    }  // end of searching of pre-evaluated fitness for current individual
+    if (!is_evaluated) error_status_ = kError;
+    Real f_best = evaluated_fitness_for_current_vectors_.front().first;
+    Real f_crossover_u = FitnessFunction_p(crossover_u, param);
+    bool is_success = f_crossover_u > f_current
+      || f_crossover_u == f_best;  //Selected for maxima search
+    if (is_find_minimum_) is_success = !is_success;
+    if (!is_success) {
+      // Case of current x and f were new for current generation.
+      x_vectors_next_generation_[i] = x_vectors_current_[i];
+      for (auto &f : evaluated_fitness_for_next_generation_) {
+        if (f.second == i) f.first = f_current;
+      }  // end of saving old fitness value in new generation.
+    } else {  // if is_success == true
+      x_vectors_next_generation_[i] = crossover_u;
+      for (auto &f : evaluated_fitness_for_next_generation_) 
+        if (f.second == i) f.first = f_crossover_u;
+      to_be_archived_best_A_.push_back(x_vectors_current_[i]);
+      successful_mutation_parameters_S_F_.push_back(mutation_F_[i]);
+      successful_crossover_parameters_S_CR_.push_back(crossover_CR_[i]);
+      // if (process_rank_ == kOutput)
+      //   printf("n%li f_new=%4.2f\n",i,f_crossover_u);
+      //PrintSingleVector(crossover_u);
+    }  // end of dealing with success crossover
+    return kDone;
+  } // end of int SubPopulation::Selection(std::vector<Real> crossover_u);
+    
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::ArchiveCleanUp() {
+    auto archived_last=archived_best_A_.end();
+    archived_best_A_.splice(archived_last, to_be_archived_best_A_);
+    auto size_A = archived_best_A_.size();
+    long initial_diff = size_A - subpopulation_;
+    // if (process_rank_ == kOutput)
+    //   printf("diff = %li size_A=%li subpop=%li \n ", initial_diff, size_A, subpopulation_);
+    // if (initial_diff < 1) return kDone;
+    // if (process_rank_ == kOutput)
+    //   printf("diff = %li size_A=%li \n ", initial_diff, size_A);
+    for (long i = 0; i < initial_diff; ++i) {
+      long index_to_remove = randint(0, size_A - 1);
+      auto element_A = archived_best_A_.begin();
+      std::advance(element_A, index_to_remove);
+      archived_best_A_.erase(element_A);
+      --size_A;
+    }
+    const unsigned long new_size_A = archived_best_A_.size();
+    if (new_size_A > subpopulation_) error_status_ = kError; //TODO change kError to throw exception
+    return kDone;
+  } // end of int SubPopulation:: ArchiveCleanUp();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::Adaption()  {    
+    long elements = 0;
+    Real sum = 0.0;
+    for (auto CR : successful_crossover_parameters_S_CR_) {
+      sum += CR;
+      ++elements;
+    }  // end of collecting data for mean CR
+    if (elements == 0) return kDone;
+    // Should never be reached for symmetric filling of S_CR and S_F.
+    if (successful_mutation_parameters_S_F_.size() == 0) return kDone;  
+    Real mean_a_CR = sum / static_cast<Real>(elements);
+    // Original JADE adaption of mu_CR.
+    if (!isPMCRADE_) {
+      adaptor_crossover_mu_CR_ =
+        (1 - adaptation_frequency_c_) *  adaptor_crossover_mu_CR_
+        + adaptation_frequency_c_ * mean_a_CR;
+    } else {
+      // PMCRADE patch for mu_CR
+      Real std_S_CR = 0;
+      for (auto CR : successful_crossover_parameters_S_CR_)
+        std_S_CR += pow2(CR - mean_a_CR);
+      std_S_CR = sqrt(std_S_CR / static_cast<Real>(elements));
+      const Real PMCRADE_const = 0.07;
+      if (std_S_CR < PMCRADE_const) {
+        adaptor_crossover_mu_CR_ =
+          (1 - adaptation_frequency_c_) *  adaptor_crossover_mu_CR_
+          + adaptation_frequency_c_ * mean_a_CR;
+      } else {
+        Real mean_pow2_CR = 0.0;
+        for (auto CR : successful_crossover_parameters_S_CR_)
+          mean_pow2_CR += pow2(CR);
+        mean_pow2_CR = sqrt(mean_pow2_CR/static_cast<Real>(elements));
+        adaptor_crossover_mu_CR_ =
+          (1 - adaptation_frequency_c_) *  adaptor_crossover_mu_CR_
+          + adaptation_frequency_c_ * mean_pow2_CR;
+      }
+      // end of PMCRADE patch
+    }  // end of if isPMCRADE_
+    Real sum_F = 0.0, sum_F2 = 0.0;
+    for (auto F : successful_mutation_parameters_S_F_) {
+      sum_F += F;
+      sum_F2 += F*F;
+    }  // end of collection data for Lehmer mean F
+    Real mean_l_F = sum_F2 / sum_F;
+    adaptor_mutation_mu_F_ =
+      (1 - adaptation_frequency_c_) *  adaptor_mutation_mu_F_
+      + adaptation_frequency_c_ * mean_l_F;
+    return kDone;
+  } // end of int SubPopulation::Adaption();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::RunOptimization() {
+    //if (process_rank_ == kOutput) printf("Start optimization..\n");
+    if (error_status_) return error_status_;
+    adaptor_mutation_mu_F_ = 0.5;
+    adaptor_crossover_mu_CR_ = 0.5;
+    archived_best_A_.clear();
+    CreateInitialPopulation();
+    x_vectors_next_generation_ = x_vectors_current_;
+    EvaluateCurrentVectors();
+    evaluated_fitness_for_next_generation_ =
+      evaluated_fitness_for_current_vectors_;
+    for (long g = 0; g < total_generations_max_; ++g) {
+      if (process_rank_ == kOutput && g%100 == 0) printf("%li\n",g);      
+      to_be_archived_best_A_.clear();
+      successful_mutation_parameters_S_F_.clear();
+      successful_crossover_parameters_S_CR_.clear();        
+      // //debug section
+      // if (process_rank_ == kOutput)
+      //   printf("==============  Generation %li =============\n", g);
+      PrintPopulation();      
+      //PrintEvaluated();
+			std::cin.get();
+      //end of debug section
+      for (unsigned long i = 0; i < subpopulation_; ++i) {
+        SetCRiFi(i);
+        std::vector<Real> mutated_v, crossover_u;
+        mutated_v = Mutation(i);
+        crossover_u = Crossover(mutated_v, i);
+        Selection(crossover_u, i);
+      }  // end of for all individuals in subpopulation
+      ArchiveCleanUp();
+      Adaption();
+      x_vectors_current_.swap(x_vectors_next_generation_);
+      evaluated_fitness_for_current_vectors_
+        .swap(evaluated_fitness_for_next_generation_);
+      SortEvaluatedCurrent();
+      if (error_status_) return error_status_;
+    }  // end of stepping generations
+    PrintPopulation();      
+    PrintEvaluated();
+    return kDone;
+  }  // end of int SubPopulation::RunOptimization()
+
+	int SubPopulation::RunOptimization(std::ostream &out, void* param) {
+    //if (process_rank_ == kOutput) printf("Start optimization..\n");
+		if (error_status_) return error_status_;
+		adaptor_mutation_mu_F_ = 0.5;
+		adaptor_crossover_mu_CR_ = 0.5;
+		archived_best_A_.clear();
+		CreateInitialPopulation();
+		x_vectors_next_generation_ = x_vectors_current_;
+		EvaluateCurrentVectors(param);
+		evaluated_fitness_for_next_generation_ = evaluated_fitness_for_current_vectors_;
+		for (long g = 0; g < total_generations_max_; ++g) {
+			if (process_rank_ == kOutput && g%100 == 0) printf("%li\n",g);      
+			to_be_archived_best_A_.clear();
+			successful_mutation_parameters_S_F_.clear();
+			successful_crossover_parameters_S_CR_.clear();        
+
+      PrintPopulation(g, out);
+
+      for (unsigned long i = 0; i < subpopulation_; ++i) {
+				SetCRiFi(i);
+				std::vector<Real> mutated_v, crossover_u;
+				mutated_v = Mutation(i);
+				crossover_u = Crossover(mutated_v, i);
+				Selection(crossover_u, param, i);
+      }  // end of for all individuals in subpopulation
+			ArchiveCleanUp();
+			Adaption();
+			x_vectors_current_.swap(x_vectors_next_generation_);
+			evaluated_fitness_for_current_vectors_
+				.swap(evaluated_fitness_for_next_generation_);
+			SortEvaluatedCurrent();
+			if (error_status_) return error_status_;
+		}  // end of stepping generations
+		PrintPopulation();   
+		PrintEvaluated();
+		return kDone;
+	}  // end of int SubPopulation::RunOptimization()
+
+	 // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::Crossover(std::vector<Real> mutation_v, long i) {
+    const Real CR_i = crossover_CR_[i];
+    std::vector<Real> crossover_u, x_current;
+    crossover_u.resize(dimension_);
+    x_current = x_vectors_current_.at(i);
+    unsigned long j_rand = randint(0, dimension_ - 1);
+    for (unsigned long c = 0; c < dimension_; ++c) {
+      if (c == j_rand || rand(0,1) < CR_i)
+        crossover_u[c] = mutation_v[c];
+      else
+        crossover_u[c] = x_current[c];
+    }
+    // //debug section
+    // if (process_rank_ == kOutput)
+    //   printf("x -> v -> u with CR_i=%4.2f j_rand=%li\n", CR_i, j_rand);
+    // PrintSingleVector(mutation_v);
+    // PrintSingleVector(x_current);
+    // PrintSingleVector(crossover_u);
+    // //end of debug section
+    return crossover_u;
+  } // end of  std::vector<Real> SubPopulation::Crossover();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::Mutation(long i) {
+    std::vector<Real> mutation_v, x_best_current, x_random_current,
+      x_random_archive_and_current, x_current;    
+    x_current = x_vectors_current_.at(i);
+    x_best_current = GetXpBestCurrent();
+    // //debug
+    // if (process_rank_ == kOutput) printf("x_best: ");
+    // PrintSingleVector(x_best_current);
+    long index_of_random_current = -1;
+    x_random_current = GetXRandomCurrent(&index_of_random_current, i);
+    // //debug
+    // if (process_rank_ == kOutput) printf("x_random: ");
+    // PrintSingleVector(x_random_current);    
+    x_random_archive_and_current = GetXRandomArchiveAndCurrent(index_of_random_current, i);
+    // //debug
+    // if (process_rank_ == kOutput) printf("x_random with archive: ");
+    // PrintSingleVector(x_random_archive_and_current);
+    mutation_v.resize(dimension_);
+    Real F_i = mutation_F_[i];
+    for (unsigned long c = 0; c < dimension_; ++c) {
+      // Mutation
+      mutation_v[c] = x_current[c]
+        + F_i * (x_best_current[c] - x_current[c])
+        + F_i * (x_random_current[c]
+                 - x_random_archive_and_current[c]);
+      // Bounds control
+      if (mutation_v[c] > x_ubound_[c])
+        mutation_v[c] = (x_ubound_[c] + x_current[c])/2;
+      if (mutation_v[c] < x_lbound_[c])
+        mutation_v[c] = (x_lbound_[c] + x_current[c])/2;
+    }
+    // //debug section
+    // int isSame = 888, isSame2 = 7777, isSame3 =11111;
+    // for (long c = 0; c < dimension_; ++c) {
+    //   Real norm = std::abs(mutation_v[c] - x_current[c]);
+    //   Real norm2 = std::abs(x_random_current[c] - x_current[c]);
+    //   Real norm3 = std::abs(x_random_current[c]
+    //                           - x_random_archive_and_current[c]);
+    //   if ( norm  > 0.0001) isSame = 0;
+    //   if ( norm2  > 0.0001) isSame2 = 0;
+    //   if ( norm3  > 0.0001) isSame3 = 0;
+    // }
+    // if (process_rank_ == kOutput) printf("mutation_v%i%i%i:  ",isSame, isSame2, isSame3);
+    // PrintSingleVector(mutation_v);
+    // if (process_rank_ == kOutput) printf("current _v%i%i%i:  ",isSame, isSame2, isSame3);
+    // PrintSingleVector(x_current);    
+    // if (process_rank_ == kOutput)
+    //   printf("  -> f = %4.2f                                    F_i=%4.2f\n",
+    //          FitnessFunction(mutation_v), F_i);
+    // //end of debug section
+    return mutation_v;
+  } // end of std::vector<Real> SubPopulation::Mutation();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::GetXpBestCurrent() {
+    const unsigned long n_best_total = static_cast<long>
+      (floor(subpopulation_ * best_share_p_ ));
+    if (n_best_total == subpopulation_) error_status_ = kError; //TODO change kError to throw exception
+    long best_n = randint(0, n_best_total);
+    long best_n_index = -1, i = 0;
+    for (auto x : evaluated_fitness_for_current_vectors_) {
+      if (i == best_n) {
+        best_n_index = x.second;
+        break;
+      }
+      ++i;
+    }
+    return x_vectors_current_.at(best_n_index);
+  }  // end of std::vector<Real> SubPopulation::GetXpBestCurrent();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::GetXRandomCurrent(long *index, long forbidden_index) {
+    long random_n = randint(0, subpopulation_-1);
+    while (random_n == forbidden_index) random_n = randint(0, subpopulation_-1);
+    (*index) = random_n;
+    return x_vectors_current_.at(random_n);
+  }  // end of std::vector<Real> SubPopulation::GetXRandomCurrent()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::GetXRandomArchiveAndCurrent
+  (unsigned long forbidden_index1, unsigned long forbidden_index2) {
+    long archive_size = archived_best_A_.size();
+    unsigned long random_n = randint(0, subpopulation_ + archive_size - 1);
+    while (random_n == forbidden_index1 || random_n == forbidden_index2)
+      random_n = randint(0, subpopulation_ + archive_size - 1);
+    if (random_n < subpopulation_) return x_vectors_current_.at(random_n);
+    random_n -= subpopulation_;
+    unsigned long i = 0;
+    for (auto x : archived_best_A_) {
+      if (i == random_n) {
+        // //debug
+        // if (process_rank_ == kOutput) printf("Using Archive!!\n");
+        return x;
+      }
+      ++i;
+    }  // end of selecting from archive
+    error_status_ = kError; //TODO change kError to throw exception
+    std::vector<Real> x;
+    return x;    
+  }  // end of std::vector<Real> SubPopulation::GetXRandomArchiveAndCurrent()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::SetCRiFi(long i) {
+    long k = 0;
+    while (1) {
+      mutation_F_[i] = randc(adaptor_mutation_mu_F_, 0.1);
+      if (mutation_F_[i] > 1) {
+        mutation_F_[i] = 1;
+        break;
+      }
+      if (mutation_F_[i] > 0) break;
+      ++k;
+      if (k > 10) {
+        mutation_F_[i] = 0.001;
+        break;
+      }
+      if (k > 100) printf("k");
+    }
+    crossover_CR_[i] = randn(adaptor_crossover_mu_CR_,0.1);
+    if (crossover_CR_[i] > 1) crossover_CR_[i] = 1;
+    if (crossover_CR_[i] < 0) crossover_CR_[i] = 0;    
+    return kDone;
+  }  // end of int SubPopulation::SetCRiFi(long i)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  /// %todo Change returned kError to meangfull error code.
+  int SubPopulation::Init(long total_population, long dimension) {// NOLINT
+    total_population_  = total_population;
+    if (total_population_ < 1)
+      throw std::invalid_argument("You should set population > 0!");
+    dimension_ = dimension;
+    if (dimension_ < 1) 
+      throw std::invalid_argument("You should set dimension > 0!");
+    MPI_Comm_rank(MPI_COMM_WORLD, &process_rank_);
+    MPI_Comm_size(MPI_COMM_WORLD, &number_of_processes_);
+    if (process_rank_ < 0)
+      throw std::invalid_argument("MPI problem: process_rank_ < 0!");
+    if (number_of_processes_ < 1)
+      throw std::invalid_argument("MPI problem: number_of_processes_ < 1!");
+    std::random_device rd;
+    generator_.seed(rd());
+
+    // //debug
+    // CheckRandom();
+
+    subpopulation_ = total_population;
+
+    current_generation_ = 0;
+    x_vectors_current_.resize(subpopulation_);
+    for (auto &x : x_vectors_current_) x.resize(dimension_);
+    x_vectors_next_generation_.resize(subpopulation_);
+    for (auto &x : x_vectors_next_generation_) x.resize(dimension_);
+    evaluated_fitness_for_current_vectors_.resize(subpopulation_);
+    // //debug
+    // if (process_rank_ == kOutput) printf("%i, x1 size = %li \n", process_rank_, x_vectors_current_.size());
+    x_lbound_.resize(dimension_);
+    x_ubound_.resize(dimension_);
+    mutation_F_.resize(subpopulation_);
+    crossover_CR_.resize(subpopulation_);
+    return kDone;
+  }  // end of void SubPopulation::Init()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::PrintPopulation() {
+		if (process_rank_ == kOutput) {
+			printf("\n");	
+			for (auto x : evaluated_fitness_for_current_vectors_) {
+				long n = x.second;
+				Real fitness = x.first;
+				printf("%6.2f:% 3li||", fitness, n);
+				for (unsigned long c = 0; c < dimension_; ++c) 
+					printf(" %+7.2f ", x_vectors_current_[n][c]);
+				printf("\n");      
+			}
+      // for (long i = 0; i < subpopulation_; ++i) {
+      //   printf("n%li:", i);
+      //   for (long c = 0; c < dimension_; ++c) {
+      //     printf(" %5.2f ", x_vectors_current_[i][c]);
+      //   }
+      //   printf("\n");
+      // }  // end of for each individual
+    }  // end of if output
+    return kDone;
+  }  // end of int SubPupulation::PrintPopulation()
+
+	int SubPopulation::PrintPopulation(long niter, std::ostream &out) {
+		if (process_rank_ == kOutput) {
+			//printf("\n");
+			auto p = evaluated_fitness_for_current_vectors_.begin();
+			out << niter << " " << (*p).first;
+			for (unsigned long c = 0; c < dimension_; ++c) out << " " << x_vectors_current_[(*p).second][c];
+			out << std::endl;
+		}  // end of if output
+		return kDone;
+	}  // end of int SubPupulation::PrintPopulation()
+		 
+		 // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::PrintEvaluated() {
+    if (process_rank_ == kOutput) {
+      for (auto x : evaluated_fitness_for_current_vectors_)
+        printf("%li:%4.2f  ", x.second, x.first);
+      printf("\n");
+    }  // end of if output
+    return kDone;
+  }  // end of int SubPopulation::PrintEvaluated()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::PrintSingleVector(std::vector<Real> x) {
+    if (process_rank_ == kOutput) {
+      for (auto c : x) printf("%5.2f ", c);
+      printf("\n");
+    }  // end of output
+    return kDone;
+  }  // end of int SubPopulation::PrintSingleVector(std::vector<Real> x)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void SubPopulation::SetFeed(std::vector<std::vector<Real> > x_feed_vectors) {
+    isFeed_  = true;
+    if (dimension_ < 1) 
+      throw std::invalid_argument("You should set dimension before feed!");
+    for (auto x : x_feed_vectors)
+      if (x.size() != dimension_) 
+	throw std::invalid_argument
+	  ("Feed and optimization dimenstion should be the same size!");
+    x_feed_vectors_.clear();
+    for (auto &x : x_feed_vectors)
+      x_feed_vectors_.push_back(x);
+  }  // end of void SubPopulation::SetFeed()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::CreateInitialPopulation() {
+    if ((subpopulation_ - x_feed_vectors_.size()) <1)
+      throw std::invalid_argument("Too large feed!");
+    x_vectors_current_.resize(subpopulation_ - x_feed_vectors_.size());
+    for (auto &x : x_vectors_current_)
+      for (unsigned long i = 0; i < dimension_; ++i) {
+        if (x_lbound_[i] > x_ubound_[i])
+	  throw std::invalid_argument("Wrong order of bounds!");
+        x[i] = rand(x_lbound_[i], x_ubound_[i]);                            // NOLINT
+      }  // end of for each dimension
+    // //debug
+    // for (auto x : x_vectors_current_[0]) if (process_rank_ == kOutput) printf("%g ",x);
+    for (auto x: x_feed_vectors_) {
+      x_vectors_current_.push_back(x);
+        if (process_rank_ == kOutput) {
+	  printf("--=-- Feed:\n");
+	  for (auto index:x) printf(" %+7.2f", index);
+	  printf("\n");
+	}	
+    }
+    if (x_vectors_current_.size() != subpopulation_)
+      throw std::invalid_argument("Population is not full after feed!");	
+    x_feed_vectors_.clear();
+    // if (process_rank_ == kOutput) {
+    //   printf("==--== Initial population:\n");    
+    //   EvaluateCurrentVectors();
+    //   PrintPopulation();
+    // }
+    return kDone;
+  }  // end of int SubPopulation::CreateInitialPopulation()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::SortEvaluatedCurrent() {
+    evaluated_fitness_for_current_vectors_
+      .sort([=](const std::pair<Real, long>& a,                     // NOLINT
+               const std::pair<Real, long>& b) {                   // NOLINT
+              bool cmp = a.first < b.first;
+              if (is_find_minimum_) return cmp;
+              return !cmp;
+            });
+    return kDone;
+  }  // end of int SubPopulation::SortEvaluatedCurrent()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::EvaluateCurrentVectors() {
+    evaluated_fitness_for_current_vectors_.clear();
+    for (unsigned long i = 0; i < subpopulation_; ++i) {                        // NOLINT
+      auto tmp = std::make_pair(FitnessFunction(x_vectors_current_[i]), i);
+      evaluated_fitness_for_current_vectors_.push_back(tmp);
+    }
+    SortEvaluatedCurrent();
+    // //debug
+    // if (process_rank_ == kOutput) printf("\n After ");
+    // for (auto val : evaluated_fitness_for_current_vectors_)
+    //   if (process_rank_ == kOutput) printf("%g ", val.first);
+    return kDone;
+  }  // end of int SubPopulation::EvaluateCurrentVectors()
+  
+  int SubPopulation::EvaluateCurrentVectors(void* param) {
+    evaluated_fitness_for_current_vectors_.clear();
+    for (unsigned long i = 0; i < subpopulation_; ++i) {                        // NOLINT
+      auto tmp = std::make_pair(FitnessFunction_p(x_vectors_current_[i], param), i);
+      evaluated_fitness_for_current_vectors_.push_back(tmp);
+    }
+    SortEvaluatedCurrent();
+    // //debug
+    // if (process_rank_ == kOutput) printf("\n After ");
+    // for (auto val : evaluated_fitness_for_current_vectors_)
+    //   if (process_rank_ == kOutput) printf("%g ", val.first);
+    return kDone;
+  }  // end of int SubPopulation::EvaluateCurrentVectors()
+     // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void SubPopulation::SetAllBoundsVectors
+  (std::vector<Real> lbound, std::vector<Real> ubound) {
+    x_lbound_.clear();
+    x_ubound_.clear();    
+    for (auto x : lbound) x_lbound_.push_back(x);
+    for (auto x : ubound) x_ubound_.push_back(x);
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::SetAllBounds(Real lbound, Real ubound) {
+    if (lbound >= ubound) {
+      error_status_ = kError;
+      return kError; //TODO change kError to throw exception
+    }
+    for (auto &x : x_lbound_) x = lbound;
+    for (auto &x : x_ubound_) x = ubound;
+    return kDone;
+  }  // end of int SubPopulation::SetAllBounds(Real lbound, Real ubound)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  Real SubPopulation::randn(Real mean, Real stddev) {
+    std::normal_distribution<Real> distribution(mean, stddev);
+    return distribution(generator_);
+  }  // end of Real SubPopulation::randn(Real mean, Real stddev)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  Real SubPopulation::randc(Real location, Real scale) {
+    std::cauchy_distribution<Real> distribution(location, scale);
+    return distribution(generator_);
+  }  // end of Real SubPopulation::randc(Real location, Real scale)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  long SubPopulation::randint(long lbound, long ubound) {    // NOLINT
+    std::uniform_int_distribution<long> distribution(lbound, ubound);  // NOLINT
+    return distribution(generator_);
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  Real SubPopulation::rand(Real lbound, Real ubound) {                // NOLINT
+    std::uniform_real_distribution<Real> distribution(lbound, ubound);
+    return distribution(generator_);
+  }  // end of Real rand(Real lbound, Real ubound)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void SubPopulation::CheckRandom() {
+    std::map<int, int> hist_n, hist_c, hist_int, hist;
+    int factor = 1000;
+    for (int n = 0; n < 100*factor; ++n) {
+      ++hist_n[std::round(randn(0, 2))];
+      ++hist_c[std::round(randc(0, 2))];
+      ++hist_int[std::round(randint(0, 10))];
+      ++hist[std::round(rand(0, 15))];                                      // NOLINT 
+    }
+    if (process_rank_ == kOutput) {
+      printf("Normal (0,2)\n");
+      for (auto p : hist_n) {
+        if (p.second > factor)
+          printf("%i: % 4i %s\n", process_rank_, p.first,
+                 std::string(p.second/factor, '*').c_str());
+      }  // end of for p in hist_n
+      printf("Cauchy (0,2)\n");
+      for (auto p : hist_c) {
+        if (p.second > factor)
+          printf("%i: % 4i %s\n", process_rank_, p.first,
+                 std::string(p.second/factor, '*').c_str());
+      }  // end of for p in hist_c
+      printf("Uniform int (0,10)\n");
+      for (auto p : hist_int) {
+        if (p.second > factor)
+          printf("%i: % 4i %s\n", process_rank_, p.first,
+                 std::string(p.second/factor, '*').c_str());
+      }  // end of for p in hist_int
+      printf("Uniform Real (0,15)\n");
+      for (auto p : hist) {
+        if (p.second > factor)
+          printf("%i: % 4i %s\n", process_rank_, p.first,
+                 std::string(p.second/factor, '*').c_str());
+      }  // end of for p in hist
+    }  //  end of if current process_rank_ == kOutput
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::SetBestShareP(Real p) {
+    if (p < 0 || p > 1) {
+      error_status_ = kError;
+      return kError; //TODO change kError to throw exception
+    }
+    best_share_p_ = p;
+    return kDone;
+  }  // end of int SubPopulation::SetBestShareP(Real p)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::SetAdapitonFrequencyC(Real c) {
+    if (c < 0 || c > 1) {
+      error_status_ = kError;
+      return kError; //TODO change kError to throw exception
+    }
+    adaptation_frequency_c_ = c;
+    return kDone;
+  }  // end of int SubPopulation::SetAdapitonFrequency(Real c)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::SetDistributionLevel(int level) {
+    distribution_level_ = level;
+    return kDone;
+  }
+  // end of int SubPopulation::SetDistributionLevel(int level)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::PrintParameters(std::string comment) {
+    if (process_rank_ == 0) {
+      printf("#%s dim=%li NP=%li(of %li) p=%4.2f c=%4.2f generation=%li\n",
+             comment.c_str(),dimension_, subpopulation_, total_population_,
+             best_share_p_, adaptation_frequency_c_, total_generations_max_
+             );
+      fflush(stdout);
+    }  // end of output
+    return kDone;
+  }  // end of int SubPopulation::PrintParameters(std::string comment)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::PrintResult(std::string comment) {    
+    if (distribution_level_ == 0) {      
+      auto x = evaluated_fitness_for_current_vectors_.front();
+      std::vector<Real> to_send {x.first};
+      //printf("%8.5g\n  ", x.first);
+      AllGatherVectorDouble(to_send);
+      if (process_rank_ == 0) {
+        Real sum = 0;
+        Real size = static_cast<Real>(recieve_Real_.size());
+        for (auto x : recieve_Real_) sum += x;
+        Real mean = sum/size;
+        Real sigma = 0;
+        for (auto x : recieve_Real_) sigma += pow2(x - mean);
+        sigma = sqrt(sigma/size);
+        printf("%s gen%li, mean %4.1e (stddev %4.1e = %3.2g %%) runs(%g)\n",
+               comment.c_str(), total_generations_max_,  mean,sigma,
+               sigma*100.0/mean,size);
+        // for (auto x : recieve_Real_)
+        //   printf("%18.15g\n", x);
+      }
+    }
+    return kDone;
+  }  // end of int SubPopulation::PrintResult()
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::GetFinalFitness() {
+    recieve_Real_.clear();    
+    if (distribution_level_ == 0) {      
+      auto x = evaluated_fitness_for_current_vectors_.front();
+      std::vector<Real> to_send {x.first};
+      AllGatherVectorDouble(to_send);      
+    }
+    return recieve_Real_;
+  }  // end of sdt::vector<Real> SubPopulation::GetFinalFitness();
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::GetBest(Real *best_fitness) {
+    auto x = evaluated_fitness_for_current_vectors_.front();
+    (*best_fitness) = x.first;
+    return x_vectors_current_[x.second];
+  }  // end of std::vector<Real> SubPopulation::GetBest(Real *best_fitness)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  std::vector<Real> SubPopulation::GetWorst(Real *worst_fitness) {
+    auto x = evaluated_fitness_for_current_vectors_.back();
+    (*worst_fitness) = x.first;
+    return x_vectors_current_[x.second];
+  }  // end of std::vector<Real> SubPopulation::GetWorst(Real *worst_fitness)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::AllGatherVectorDouble(std::vector<Real> to_send) {
+    long size_single = to_send.size();
+    long size_all = size_single * number_of_processes_;
+    recieve_Real_.clear();
+    recieve_Real_.resize(size_all);
+    MPI_Allgather(&to_send.front(), size_single, MPI_DOUBLE,
+                  &recieve_Real_.front(), size_single, MPI_DOUBLE,
+                  MPI_COMM_WORLD);
+    return kDone;
+  }  // end of int SubPopulation::AllGatherVectorDouble(std::vector<Real> to_send);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  int SubPopulation::AllGatherVectorLong(std::vector<long> to_send) {
+    long size_single = to_send.size();
+    long size_all = size_single * number_of_processes_;
+    recieve_long_.clear();
+    recieve_long_.resize(size_all);
+    MPI_Allgather(&to_send.front(), size_single, MPI_LONG,
+                  &recieve_long_.front(), size_single, MPI_LONG,
+                  MPI_COMM_WORLD);
+    return kDone;
+  }  // end of int SubPopulation::AllGatherVectorLong(std::vector<long> to_send);
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+}  // end of namespace jade

jade.h

@@ -0,0 +1,231 @@
+#ifndef SRC_JADE_H_
+#define SRC_JADE_H_
+///
+/// @file   jade.h
+/// @author Ladutenko Konstantin <kostyfisik at gmail (.) com>
+/// @date   Thu Aug 15 19:21:57 2013
+/// @copyright 2013 Ladutenko Konstantin
+/// @section LICENSE
+/// This file is part of JADE++.
+///
+/// JADE++ is free software: you can redistribute it and/or modify
+/// it under the terms of the GNU General Public License as published by
+/// the Free Software Foundation, either version 3 of the License, or
+/// (at your option) any later version.
+///
+/// JADE++ is distributed in the hope that it will be useful,
+/// but WITHOUT ANY WARRANTY; without even the implied warranty of
+/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+/// GNU General Public License for more details.
+///
+/// You should have received a copy of the GNU General Public License
+/// along with JADE++.  If not, see <http://www.gnu.org/licenses/>.
+
+/// @brief JADE++ is a free (GPLv3+) high performance implementation of
+/// adaptive differential evolution optimization algorithm from
+/// Jingqiao Zhang and Arthur C. Sanderson book 'Adaptive Differential
+/// Evolution. A Robust Approach to Multimodal Problem Optimization'
+/// Springer, 2009. Crossover rate was patched according to PMCRADE
+/// approach supposed by Jie Li, Wujie Zhu, Mengjun Zhou, and Hua Wang
+/// in 'Power Mean Based Crossover Rate Adaptive Differential
+/// Evolution' in H. Deng et al. (Eds.): AICI 2011, Part II, LNAI
+/// 7003, pp. 34–41, 2011
+
+#include <random>
+#include <utility>
+#include <list>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include "sph_bessel.h"
+
+namespace jade {
+    /// @brief Population controlled by single MPI process.
+  class SubPopulation {
+  public:
+      /// @brief Externaly defined fitness function, used by pointer.
+    Real (*FitnessFunction)(std::vector<Real> x) = nullptr;
+    Real (*FitnessFunction_p)(std::vector<Real>&, void*) = nullptr;
+    /// @brief Class initialization.
+    int Init(long total_population, long dimension);              // NOLINT
+      /// @brief Vizualize used random distributions (to do manual check).
+    void SetFeed(std::vector<std::vector<Real> > x_feed_vectors);
+    void CheckRandom();
+      /// @brief Find optimum value of fitness function.
+	  int RunOptimization();
+    int RunOptimization(std::ostream&, void* param);
+      /// @brief Set maximum number of generations used for optimization.
+    void SetTotalGenerationsMax(long gen) {total_generations_max_ = gen;} // NOLINT
+      /// @brief Select if to find global minimum or maximum of fitness function.
+    void SetTargetToMinimum() {is_find_minimum_ = true;}
+    void SetTargetToMaximum() {is_find_minimum_ = false;}
+      /// @brief Set adaption parameters.
+    int SetBestShareP(Real p);
+    int SetAdapitonFrequencyC(Real c);
+      /// @brief Set level of algorithm distribution.
+      /// 0 - no distribution, each MPI process acts independantly.
+    int SetDistributionLevel(int level);
+      /// @brief Set same search bounds for all components of fitness
+      /// function input vector.
+    int SetAllBounds(Real lbound, Real ubound);
+    void SetAllBoundsVectors(std::vector<Real> lbound, std::vector<Real> ubound);
+      /// @brief Print Optimization parameters.
+    int PrintParameters(std::string comment);
+      /// @brief Print final result
+    int PrintResult(std::string comment);
+    std::vector<Real> GetFinalFitness();
+    std::vector<Real> GetBest(Real *best_fitness);
+    std::vector<Real> GetWorst(Real *worst_fitness);
+    int ErrorStatus() {return error_status_;};
+    void SwitchOffPMCRADE(){isPMCRADE_ = false;};
+  
+  private:
+    bool isPMCRADE_ = true;
+    bool isFeed_ = false;
+    int CreateInitialPopulation();
+    int PrintPopulation();
+		int PrintPopulation(long, std::ostream&); 
+    int PrintEvaluated();
+    int PrintSingleVector(std::vector<Real> x);
+    int SortEvaluatedCurrent();
+      /// @brief Apply fitness function to current population.
+    int EvaluateCurrentVectors();
+    int EvaluateCurrentVectors(void* param);
+      /// @brief Generate crossover and mutation factors for current individual
+    int SetCRiFi(long i);
+      /// @name Main algorithm steps.
+      // @{
+    int Selection(std::vector<Real> crossover_u, long individual_index);
+    int Selection(std::vector<Real> crossover_u, void* param, long individual_index);
+    int ArchiveCleanUp();
+    int Adaption();
+    std::vector<Real> Mutation(long individual_index);
+    std::vector<Real> Crossover(std::vector<Real> mutated_v,
+                                  long individual_index);
+    // @}
+    /// @name Other algorithm steps.
+    // @{
+    std::vector<Real> GetXpBestCurrent();
+    /// @brief Returns random vector from current population and
+    /// vector`s index.
+    std::vector<Real> GetXRandomCurrent(long *index,
+                                          long forbidden_index);
+    std::vector<Real> GetXRandomArchiveAndCurrent(
+        unsigned long forbidden_index1, unsigned long forbidden_index2);
+    // @}
+    /// @name Population, individuals and algorithm .
+    // @{
+    /// @brief Search minimum or maximum of fitness function.
+    bool is_find_minimum_ = true;
+    /// @brief Maximum number of generations used for optimization.
+    long total_generations_max_ = 0;                                   // NOLINT
+    /// @brief Total number of individuals in all subpopulations.
+    long total_population_ = 0;                                        // NOLINT
+    /// @brief Number of individuals in subpopulation
+    unsigned long subpopulation_ = 0;                                  // NOLINT
+    /* /// @brief All individuals are indexed. First and last index of */
+    /* /// individuals in subpopulations. */
+    /* long index_first_ = -1, index_last_ = -1;                          // NOLINT */
+    /// @brief Dimension of the optimization task (number of variables
+    /// to optimize).
+    unsigned long dimension_ = 0;                                              // NOLINT
+    /// @brief Current generation of evalution process;
+    long current_generation_ = -1;                                     // NOLINT
+    /// @brief Several feed vectors.
+    std::vector<std::vector<Real> > x_feed_vectors_;
+    /// @brief Current state vectors of all individuals in subpopulation.
+    std::vector<std::vector<Real> > x_vectors_current_;
+    /// @brief State vectors of all individuals in subpopulation in
+    /// new generation.
+    std::vector<std::vector<Real> > x_vectors_next_generation_;
+    /// @brief Sometimes sorted list of evaluated fitness function.
+    std::list<std::pair<Real, long> >                                // NOLINT
+        evaluated_fitness_for_current_vectors_;
+    /// @brief Sometimes sorted list of evaluated fitness function for
+    /// next generation.
+    std::list<std::pair<Real, long> >                                // NOLINT
+        evaluated_fitness_for_next_generation_;
+    /// @brief Archived best solutions (state vactors)
+    std::list<std::vector<Real> > archived_best_A_;
+    std::list<std::vector<Real> > to_be_archived_best_A_;
+    /// @brief Sometimes sorted list of evaluated fitness function for
+    /// best vectors.
+    std::list<std::pair<Real, long> >                                // NOLINT
+        evaluated_fitness_for_archived_best_;
+    /// @brief Low and upper bounds for x vectors.
+    std::vector<Real> x_lbound_;
+    std::vector<Real> x_ubound_;
+    /// @brief JADE+ adaption parameter for mutation factor
+    Real adaptor_mutation_mu_F_ = 0.5;
+    /// @brief JADE+ adaption parameter for crossover probability
+    Real adaptor_crossover_mu_CR_ = 0.5;
+    /// @brief Individual mutation and crossover parameters for each individual.
+    std::vector<Real> mutation_F_, crossover_CR_;
+    std::list<Real> successful_mutation_parameters_S_F_;
+    std::list<Real> successful_crossover_parameters_S_CR_;
+    /// @brief Share of all individuals in current population to be
+    /// the best, recomended value range 0.05-0.2
+    //const Real best_share_p_ = 0.12;
+    Real best_share_p_ = 0.05;
+
+    /* //debug Change it back before production!!!! */
+    /* const Real best_share_p_ = 0.3; */
+
+    /// @brief 1/c - number of generations accounted for parameter
+    /// adaption, recommended value 5 to 20 generation;
+    //const Real adaptation_frequency_c_ = 1.0/20.0;    
+    Real adaptation_frequency_c_ = 0.1;    
+    // @}
+    /// @name Random generation
+    /// Names are in notation from Jingqiao Zhang and Arthur C. Sanderson book.
+    // @{
+    /// @todo Select random generator enginge for best results in DE!
+    //std::mt19937_64 generator_;
+    std::ranlux48 generator_;
+    /// @brief randn(&mu;, &sigma^2; ) denotes a random value from a normal
+    /// distribution of mean &mu; and variance &sigma^2;
+    Real randn(Real mean, Real stddev);
+    /// @brief randc(&mu;, &delta; ) a random value from a Cauchy distribution
+    /// with location and scale parameters &mu; and &delta;
+    Real randc(Real location, Real scale);
+    /// @brief randint(1, D) is an integer randomly chosen from 1 to D
+    long randint(long lbound, long ubound);                  // NOLINT
+    /// @brief rand(a, b) is an uniform random number chosen from a to b
+    Real rand(Real lbound, Real ubound);                              // NOLINT
+    // @}
+    /// @name MPI section
+    // @{
+    int process_rank_;
+    int number_of_processes_;
+    int AllGatherVectorDouble(std::vector<Real> to_send);
+    std::vector<Real> recieve_Real_;
+    int AllGatherVectorLong(std::vector<long> to_send);
+    std::vector<long> recieve_long_;
+       
+    // @}
+    /// @brief Subpopulation status. If non-zero than some error has appeared.
+    int error_status_ = 0;  // @todo move to exceptions!
+    int distribution_level_ = 0;
+  };  // end of class SubPopulation
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  /// @brief Error codes
+  ///  
+  /// Error codes used with jade
+  /// @todo move to exceptions!
+  enum Errors {
+    /// no error
+    kDone = 0,
+    /// Unspecified (pending to be described).
+    kError
+  };  // end of enum Errors
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  const int kOutput = 0; /// Process rank to do output with printf
+  template<class T> inline T pow2(const T value) {return value*value;}
+}  // end of namespace jade
+
+#endif  // SRC_JADE_H_

main.cpp

@@ -0,0 +1,128 @@
+
+#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;
+	}
+/**/
+}

sph_bessel.cpp

@@ -0,0 +1,264 @@
+
+#include "sph_bessel.h"
+#include <limits>
+
+using namespace std;
+using namespace sph_bessel;
+
+#define DG 12
+
+Complex sph_bessel::besj(Complex z, int n) {
+	if (n == 0) return besj0(z);
+	else if (n == 1) return besj1(z);
+	//	else if (n < 0) return (n%2) ? besy(z,-n-1) : -besy(z,-n-1);
+
+	if (abs(z) < 0.6) {
+		int i, ip;
+		Real tvv;
+		Complex tc, tcc;
+
+		tc = 1.; for (i=1; i<n+1; i++) tc *= z/Real(2*i+1); i = 1; tcc = tc;
+		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
+		do {tcc += ( tc *= -Real(0.5)*z*z/Real(i*(2*(n+i)+1)) ); i++;} while (abs(tc) > tvv);
+
+		return tcc;
+	}
+	else if (abs(z)/Real(n) < 1.) {
+		int nn, i, pw;
+		Real tv1, tv2, tv3, tx, tx1, tx2, ty1, ty2;
+		Complex tc, tc1, tc2, tcc;
+
+		if (n > int(abs(z))) pw = abs(int( (n*log(0.5*_e*abs(z)/n) - 0.5*log(n*abs(z)) - _ln2)/_ln10 )) + DG;
+		else pw = DG;
+
+		tv1 = abs(z); tv2 = 2./_e/tv1; tv3 = pw*_ln10 - _ln2;
+		tx2 = (tx1 = n) + 2;
+		ty1 = tx1*log(tv2*tx1) + 0.5*log(tv1*tx1) - tv3;
+		ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
+		do {
+			tx = tx2 - ty2*(tx2-tx1)/(ty2-ty1); tx1 = tx2; tx2 = tx;
+			ty1 = ty2; ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
+		} while (fabs(tx2-tx1) > 0.5);
+		nn = int(tx2);
+
+		tc1 = 0.; tc2 = exp(-pw*_ln10);
+		for (i=nn; i>0; i--) {
+			tc = Real(2*i+1)/z*tc2 - tc1; tc1 = tc2; tc2 = tc;
+			if (i == n+1) tcc = tc2;
+		}
+		return (abs(besj1(z)) > abs(besj0(z))) ? tcc*besj1(z)/tc1 : tcc*besj0(z)/tc2;
+	}
+	else {
+		int i, tm;
+		Real tvv;
+		Complex tc, tp, tq, pp, qq;
+		i = -int(0.4343*log(abs(z)/n)) - DG; tvv = (i < -100) ? 1.e-100 : exp(2.3*i);
+		tp = pp = 1.; tm = (2*n+1)*(2*n+1); tq = qq = Real(tm-1)*(tc = Real(0.125)/z); tc *= tc; i = 1;
+		do {
+			pp += ( tp *= -tc*Real((tm - (4*i-1)*(4*i-1))*(tm - (4*i-3)*(4*i-3)))/Real(2*i*(2*i-1)) );
+			qq += ( tq *= -tc*Real((tm - (4*i+1)*(4*i+1))*(tm - (4*i-1)*(4*i-1)))/Real(2*i*(2*i+1)) );
+			i++;
+		} while ((abs(tp) > tvv) && (abs(tq) > tvv));
+		return (pp*cos(z-_pi_2*(n+1)) - qq*sin(z-_pi_2*(n+1)))/z;
+	}
+}
+
+Complex sph_bessel::besjd(Complex z, int n) {
+	if (n == 0) return besj0d(z);
+	else if (n == 1) return besj1d(z);
+	//	else if (n < 0) return (n%2) ? besy(z,-n-1) : -besy(z,-n-1);
+
+	if (abs(z) < 0.6) {
+		int i, ip;
+		Real tv, tvv;
+		Complex tc, tcc;
+
+		tc = 1/3.; for (i=2; i<n+1; i++) tc *= z/Real(2*i+1); i = 1; tcc = Real(n)*tc;
+		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
+		do {tv = abs( tc *= -Real(0.5)*z*z/Real(i*(2*(n+i)+1)) ); tcc += Real(n+2*i)*tc; i++;} while (tv > tvv);
+
+		return tcc;
+	}
+	else if (abs(z)/n < 1.) {
+		int nn, i, pw;
+		Real tv1, tv2, tv3, tx, tx1, tx2, ty1, ty2;
+		Complex tc, tc1, tc2, tcc;
+
+		if (n > int(abs(z))) pw = abs(int( (n*log(0.5*_e*abs(z)/n) - 0.5*log(n*abs(z)) - _ln2)/_ln10 )) + DG;
+		else pw = DG;
+
+		tv1 = abs(z); tv2 = 2./_e/tv1; tv3 = pw*_ln10 - _ln2;
+		tx2 = (tx1 = n) + 2;
+		ty1 = tx1*log(tv2*tx1) + 0.5*log(tv1*tx1) - tv3;
+		ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
+		do {
+			tx = tx2 - ty2*(tx2-tx1)/(ty2-ty1); tx1 = tx2; tx2 = tx;
+			ty1 = ty2; ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
+		} while (fabs(tx2-tx1) > 0.5);
+		nn = int(tx2);
+
+		tc1 = 0.; tc2 = exp(-pw*_ln10);
+		for (i=nn; i>0; i--) {
+			tc = Real(2*i+1)/z*tc2 - tc1; tc1 = tc2; tc2 = tc;
+			if (i == n+1) tcc = Real(n)/z*tc2 - tc1;
+		}
+		return (abs(besj1(z)) > abs(besj0(z))) ? tcc*besj1(z)/tc1 : tcc*besj0(z)/tc2;
+	}
+	else {
+		return (Real(n)*besj(z,n-1) - Real(n+1)*besj(z,n+1))/Real(2*n+1);
+	}
+}
+
+Complex sph_bessel::besy(Complex z, int n) {
+	if (n == 0) return besy0(z);
+	else if (n == 1) return besy1(z);
+	//	else if (n < 0) return (abs(n)%2) ? besj(z,n) : -besj(z,n);
+
+	if (abs(z) < 0.6) {
+		int i, ip;
+		Real tv, tvv;
+		Complex tc, tcc;
+
+		tc = -Real(1)/z; for (i=1; i<n+1; i++) tc *= Real(2*i-1)/z; i = 1; tcc = tc;
+		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
+		do {tv = abs( tc *= Real(0.5)*z*z/Real(2*(n-i+1)-1)/Real(i) ); tcc += tc; i++;} while (tv > tvv);
+
+		return tcc;
+	}
+	else if (abs(z)/n < 1.) {
+		Complex tc, tc1 = besy0(z), tc2 = besy1(z);
+		for (int i=2; i<n+1; i++) {
+			tc = Real(2*i-1)*tc2/z - tc1; tc1 = tc2; tc2 = tc;
+		}
+		return tc2;
+	}
+	else {
+		int i, tm;
+		Real tvv;
+		Complex tc, tp, tq, pp, qq;
+		i = -int(0.4343*log(abs(z))) - DG; tvv = (i < -100) ? 1.e-100 : exp(2.3*i);
+		tp = pp = 1.; tm = (2*n+1)*(2*n+1); tq = qq = Real(tm-1)*(tc = Real(0.125)/z); tc *= tc; i = 1;
+		do {
+			pp += ( tp *= -tc*Real((tm - (4*i-1)*(4*i-1))*(tm - (4*i-3)*(4*i-3)))/Real(2*i*(2*i-1)) );
+			qq += ( tq *= -tc*Real((tm - (4*i+1)*(4*i+1))*(tm - (4*i-1)*(4*i-1)))/Real(2*i*(2*i+1)) );
+			i++;
+		} while ((abs(tp) > tvv) && (abs(tq) > tvv));
+		return (pp*sin(z-_pi_2*(n+1)) + qq*cos(z-_pi_2*(n+1)))/z;
+	}
+}
+
+Complex sph_bessel::besyd(Complex z, int n) {
+	if (n == 0) return besy0d(z);
+	else if (n == 1) return besy1d(z);
+	//	else if (n < 0) return (abs(n)%2) ? besj(z,n) : -besj(z,n);
+
+	if (abs(z) < 0.6) {
+		int i, ip;
+		Real tvv;
+		Complex tc, tcc;
+		tc = Real(1)/z/z; for (i=1; i<n+1; i++) tc *= Real(2*i-1)/z; i = 1; tcc = Real(n+1)*tc;
+		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
+		do {tcc -= Real(2*i-n-1)*( tc *= Real(0.5)*z*z/Real(2*(n-i+1)-1)/Real(i) ); i++;} while (abs(tc) > tvv);
+		return tcc;
+	}
+	else if (abs(z)/n < 1.) {
+		Complex tc, tc1 = besy0(z), tc2 = besy1(z);
+		for (int i=2; i<n+1; i++) {
+			tc = Real(2*i-1)*tc2/z - tc1; tc1 = tc2; tc2 = tc;
+		}
+		return (tc1 - Real(n+1)/z*tc2);
+	}
+	else return (Real(n)*besy(z,n-1) - Real(n+1)*besy(z,n+1))/Real(2*n+1);
+}
+
+Complex sph_bessel::besh1(Complex z, int n) {
+	if (n == 0) return besh10(z);
+	else if (n == 1) return besh11(z);
+	//	return (besj(z,n) + j_*besy(z,n));
+
+	if (z.imag() > -1.e-16*abs(z)) {
+		Complex tc, tc1 = besh10(z), tc2 = besh11(z); int i = 1;
+		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
+		return tc2;
+	}
+	else return (Real(2)*besj(z,n) - besh2(z,n));
+}
+
+Complex sph_bessel::besh2(Complex z, int n) {
+	if (n == 0) return besh20(z);
+	else if (n == 1) return besh21(z);
+	//	return (besj(z,n) - j_*besy(z,n));
+
+	if (z.imag() < 1.e-16*abs(z)) {
+		Complex tc, tc1 = besh20(z), tc2 = besh21(z); int i = 1;
+		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
+		return tc2;
+	}
+	else return (Real(2)*besj(z,n) + besh1(z,n));
+}
+
+Complex sph_bessel::besh1d(Complex z, int n) {
+	if (n == 0) return besh10d(z);
+	else if (n == 1) return besh11d(z);
+	//	return (besjd(z,n) + j_*besyd(z,n));
+
+	if (z.imag() > -1.e-16*abs(z)) {
+		Complex tc, tc1 = besh10(z), tc2 = besh11(z); int i = 1;
+		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
+		return (tc1 - Real(n+1)*tc2/z);
+	}
+	else return (Real(2)*besjd(z,n) - besh2d(z,n));
+}
+
+Complex sph_bessel::besh2d(Complex z, int n) {
+	if (n == 0) return besh20d(z);
+	else if (n == 1) return besh21d(z);
+	//	return (besjd(z,n) - j_*besyd(z,n));
+
+	if (z.imag() < 1.e-16*abs(z)) {
+		Complex tc, tc1 = besh20(z), tc2 = besh21(z); int i = 1;
+		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
+		return (tc1 - Real(n+1)*tc2/z);
+	}
+	else return (Real(2)*besjd(z,n) + besh1d(z,n));
+}
+
+//###################################################
+
+void get_pitau_m01(Real x, Real *pi1, Real *tau0, Real *tau1, int deg_max) {
+	if (fabs(x-1) < numeric_limits<Real>::epsilon()) {
+		memset(tau0,0,(deg_max-1)*sizeof(Real));
+		for (int n=0; n<deg_max; ++n) {
+			pi1[n] = tau1[n] = 0.5*sqrt(0.5*(n+1)*(n+2)*(2*n+3));
+		}
+	}
+	else if (fabs(x+1) < numeric_limits<Real>::epsilon()) {
+		memset(tau0,0,(deg_max-1)*sizeof(Real));
+		for (int n=0; n<deg_max; n+=2) {
+			tau1[n] = -(pi1[n] = 0.5*sqrt(0.5*(n+1)*(n+2)*(2*n+3)));
+		}
+		for (int n=1; n<deg_max; n+=2) {
+			pi1[n] = -(tau1[n] = 0.5*sqrt(0.5*(n+1)*(n+2)*(2*n+3)));
+		}
+	}
+	else {
+		Real pl1, pl2, tmp;
+		Real tsin2 = Real(1)-x*x, tsin = sqrt(tsin2);
+
+		pl1 = _1sq2;
+		pl2 = sqrt(1.5)*x;
+		pi1[0] = 0.5*sqrt(3);
+		tau0[0] = -sqrt(1.5)*tsin;
+		tau1[0] = 0.5*sqrt(3)*x;
+
+		for (int n=1; n<deg_max; ++n) {
+			tmp = sqrt(Real(2*n+3))/Real(n+1)*(x*pl2*sqrt(Real(2*n+1)) - Real(n)/sqrt(Real(2*n-1))*pl1);
+			pl1 = pl2; pl2 = tmp;
+			tmp = sqrt(Real(n+1)/Real(n+2))/tsin2;
+			tau1[n] = sqrt(Real(2*n+3)/Real(2*n+1))*pl1;
+			pi1[n] = -tmp*( tau0[n] = x*pl2 - tau1[n] );
+			tau0[n] *= Real(n+1)/tsin;
+			tau1[n] = tmp*( (Real(1)+(n+1)*tsin2)*pl2 - x*tau1[n] );
+		}
+	}
+}

sph_bessel.h

@@ -0,0 +1,65 @@
+
+#ifndef _SPH_BESSEL_H
+#define _SPH_BESSEL_H
+
+#include <iostream>
+#include <complex>
+
+typedef double Real;
+typedef std::complex<Real> Complex;
+
+const Complex im1(Real(0),Real(1));
+const Real _pi = 3.14159265358979323846;
+const Real _pi_2 = 1.57079632679489661923;
+const Real _e = 2.71828182845904523536;
+const Real _ln10 = 2.30258509299404568402;
+const Real _ln2 = 0.693147180559945309417;
+const Real _1sq2 = 0.707106781186547524401;
+
+inline Complex get_pow_i1(int p) {
+	int i = p%4;
+	if (i == 0) return Real(1);
+	else if (i == 1) return im1;
+	else if (i == 2) return Real(-1);
+	else return -im1;
+}
+
+namespace sph_bessel {
+	inline Complex besj0(Complex z) {if (abs(z) < 1.e-8) return 1.; else return sin(z)/z;};
+	inline Complex besj1(Complex z) {if (abs(z) < 1.e-15) return z/Real(3); else return (sin(z)-z*cos(z))/z/z;}
+	inline Complex besj0d(Complex z) {if (abs(z) < 1.e-15) return -z/Real(3); else return (z*cos(z)-sin(z))/z/z;}
+	inline Complex besj1d(Complex z) {
+		if (abs(z) < 1.e-8) return Real(1)/Real(3);
+		else return ((z*z - Real(2))*sin(z) + Real(2)*z*cos(z))/(z*z*z);
+	}
+	inline Complex besy0(Complex z) {return -cos(z)/z;}
+	inline Complex besy1(Complex z) {return (-cos(z) - z*sin(z))/z/z;}
+	inline Complex besy0d(Complex z) {return (cos(z) + z*sin(z))/z/z;;}
+	inline Complex besy1d(Complex z) {return -( (z*z - Real(2))*cos(z) - Real(2)*z*sin(z) )/z/z/z;}
+	inline Complex besh10(Complex z) {return -im1*exp(im1*z)/z;}
+	inline Complex besh11(Complex z) {return (-z-im1)*exp(im1*z)/z/z;}
+	inline Complex besh10d(Complex z) {return (z+im1)*exp(im1*z)/z/z;}
+	inline Complex besh11d(Complex z) {return (-im1 + Real(2)*(z+im1)/z/z )*exp(im1*z)/z;}
+	inline Complex besh20(Complex z) {return im1*exp(-im1*z)/z;}
+	inline Complex besh21(Complex z) {return (-z+im1)*exp(-im1*z)/z/z;}
+	inline Complex besh20d(Complex z) {return (z-im1)*exp(-im1*z)/z/z;}
+	inline Complex besh21d(Complex z) {return (im1 + Real(2)*(z-im1)/z/z)*exp(-im1*z)/z;}
+
+	Complex besj(Complex, int);
+	Complex besjd(Complex, int);
+	Complex besy(Complex, int);
+	Complex besyd(Complex, int);
+	Complex besh1(Complex, int);
+	Complex besh2(Complex, int);
+	Complex besh1d(Complex, int);
+	Complex besh2d(Complex, int);
+
+	inline Complex bes_dzj(Complex z, int n) {return besj(z,n)+z*besjd(z,n);};
+	inline Complex bes_dzy(Complex z, int n) {return besy(z,n)+z*besyd(z,n);};
+	inline Complex bes_dzh1(Complex z, int n) {return besh1(z,n)+z*besh1d(z,n);};
+	inline Complex bes_dzh2(Complex z, int n) {return besh2(z,n)+z*besh2d(z,n);};
+}
+
+void get_pitau_m01(Real x, Real *pi1, Real *tau0, Real *tau1, int deg_max);
+
+#endif // _SPH_BESSEL_H

superdirectivity.cpp

@@ -0,0 +1,435 @@
+
+#include "superdirectivity.h"
+#include "sph_bessel.h"
+
+#include <algorithm>
+
+using namespace std;
+using namespace sph_bessel;
+
+Real get_directivity_sphml_dZx_epsRe(vector<Real> &param, void *cls) {
+	int ni = (size(param)-1)/2; // 1, ni, ni
+	sphere_ml *sph = reinterpret_cast<sphere_ml*>(cls);
+	vector<Real> kR(param.begin()+1, param.begin()+ni+1);
+	sort(kR.begin(),kR.end());
+	vector<Complex> eps;
+	for (int i=0; i<ni; ++i)
+		eps.push_back(Complex(param[ni+1+i],Real(0)));
+	eps.push_back(Complex(Real(1),Real(0)));
+	return sph->get_directivity_edipole_x(param[0],kR,eps);
+}
+
+Real get_directivity_sphml_dZx_epsRe_parallel(vector<Real> &param, void *cls) {
+	int ni = (size(param)-1)/2;
+	sphere_ml_parallel *sph = reinterpret_cast<sphere_ml_parallel*>(cls);
+	vector<Real> kR(param.begin()+1, param.begin()+ni+1);
+	sort(kR.begin(),kR.end());
+	vector<Complex> eps;
+	for (int i=0; i<ni; ++i)
+		eps.push_back(Complex(param[ni+1+i],Real(0)));
+	eps.push_back(Complex(Real(1),Real(0)));
+	return sph->get_directivity_edipole_x(param[0],kR,eps);
+}
+
+//###############################################################################
+// class sphere_ml
+
+// public members
+
+void sphere_ml::set_directivity_angles(Real theta_, Real phi_) {
+	theta = theta_; phi = phi_;
+	get_pitau_m01(cos(theta),pi1,tau0,tau1,deg_max);
+}
+
+Real sphere_ml::get_directivity_edipole_x(Real kz, const std::vector<Real> &kR, const std::vector<Complex> &eps) {
+	int ni = size(kR);
+	if (size(eps)-1 != ni) {std::cout << "get_directivity_edipole_x error\n"; return Real(0);}
+	Real D = 0.;
+	Complex kzn;
+	if (kz <= kR[0]) {
+			// s-matrix of a the multilayer
+		sphere_sm(kR[0], eps[0], eps[1], rtml_out);
+		//cout << "sphere coeff:\n";
+		//for (int i=0; i<deg_max; ++i) {cout << i << " " << rtml_out[8*i+2] << " " << rtml_out[8*i+3] << endl;}
+		for (int k=1; k<ni; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_out, rt, rtml_out);
+		}
+			// dipole amplitudes
+		kzn = kz*sqrt(eps[0]);
+		if (arg(kzn) < -1e-14) kzn = -kzn;
+		get_edipole_amp_x(kzn, ampl_out, 1);
+			// output amplitudes
+		for (int n=0; n<deg_max; ++n) {
+			ampl_out[3*n] *= rtml_out[8*n+2]; // e, m = +-1
+			ampl_out[3*n+1] *= rtml_out[8*n+3]; // h, m = +-1
+			ampl_out[3*n+2] *= rtml_out[8*n+3]; // h, m = 0
+		}
+	}
+	else if (kz <= kR[kR.size()-1]) { // the dipole is located in a layer
+			// determine a dipole layer
+		int kd = 0;
+		while (kz > kR[++kd]); // determine a first interface with radius larger than the dipole position
+			// calculate s-matrices of multilayers inside and outside the dipole position:
+		sphere_sm(kR[0], eps[0], eps[1], rtml_in);
+		for (int k=1; k<kd; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_in, rt, rtml_in);
+		}
+		sphere_sm(kR[kd], eps[kd], eps[kd+1], rtml_out);
+		for (int k=kd+1; k<ni; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_out, rt, rtml_out);
+		}
+			// dipole amplitudes
+		kzn = kz*sqrt(eps[kd]);
+		if (arg(kzn) < -1e-14) kzn = -kzn;
+		get_edipole_amp_x(kzn, ampl_in, 0);
+		get_edipole_amp_x(kzn, ampl_out, 1);
+			// output self-consistent amplitudes
+		for (int n=0; n<deg_max; ++n) {
+			ampl_out[3*n] = ( ampl_in[3*n]*rtml_in[8*n+6] + ampl_out[3*n] ) * rtml_out[8*n+2]
+										/ (Real(1) - rtml_in[8*n+6]*rtml_out[8*n]); // e, m = 1
+			ampl_out[3*n+1] = ( ampl_in[3*n+1]*rtml_in[8*n+7] + ampl_out[3*n+1] ) * rtml_out[8*n+3]
+										/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 1
+			ampl_out[3*n+2] = ( ampl_in[3*n+2]*rtml_in[8*n+7] + ampl_out[3*n+2] ) * rtml_out[8*n+3]
+										/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 0
+		}
+	}
+	else { // the dipole is located outside the multilayer
+			// s-matrix of a the multilayer
+		sphere_sm(kR[0], eps[0], eps[1], rtml_out);
+		for (int k=1; k<ni; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_out, rt, rtml_out);
+		}
+			// dipole amplitudes
+		kzn = kz*sqrt(eps[ni]);
+		if (arg(kzn) < -1e-14) kzn = -kzn;
+		get_edipole_amp_x(kzn, ampl_in, 0);
+		get_edipole_amp_x(kzn, ampl_out, 1);
+			// output amplitudes
+		for (int n=0; n<deg_max; ++n) {
+			ampl_out[3*n] += ampl_in[3*n] * rtml_out[8*n+6]; // e, m = +-1
+			ampl_out[3*n+1] += ampl_in[3*n+1] * rtml_out[8*n+7]; // h, m = +-1
+			ampl_out[3*n+2] += ampl_in[3*n+2] * rtml_out[8*n+7]; // h, m = 0
+		}
+	}
+
+	return directivity_dipole_z(ampl_out);
+}
+
+void sphere_ml::get_edipole_amp_x(Complex kz, Complex *ampl, bool out) {
+	if (!out) {
+		Real tv = -0.25/sqrt(_pi), tvn, tv1, tv2;
+		Complex zf, zfd;
+		memset(ampl,0,3*deg_max*sizeof(Complex));
+		for (int n=0; n<deg_max; ++n) {
+			tvn = tv*sqrt(Real(2*n+3));
+			zf = besh1(kz,n+1); zfd = besh1d(kz,n+1);
+			ampl[3*n] = tvn*zf; // e-polarization, m = 1
+			ampl[3*n+1] = -im1*tvn*(zfd + zf/kz); // h-polarization m = 1
+			tv = -tv;
+		}
+	}
+	else {
+		Real tv = -0.25/sqrt(_pi), tvn, tv1, tv2;
+		Complex zf, zfd;
+		memset(ampl,0,3*deg_max*sizeof(Complex));
+		for (int n=0; n<deg_max; ++n) {
+			tvn = tv*sqrt(Real(2*n+3));
+			zf = besj(kz,n+1); zfd = besjd(kz,n+1);
+			ampl[3*n] = tvn*zf; // m = -1,1
+			ampl[3*n+1] = -im1*tvn*(zfd + zf/kz);
+			tv = -tv;
+		}
+	}
+}
+
+std::vector<Complex> sphere_ml::get_last_ampl() const {
+	std::vector<Complex> res;
+	res.resize(3*deg_max);
+	for (int i=0; i<3*deg_max; ++i) res[i] = ampl_out[i];
+	return res;
+}
+
+std::vector<Complex> sphere_ml::ampl_far(Real th_rad, Real ph_rad) {
+	Complex tci = -im1, tc, ts = Real(0);
+	std::vector<Complex> E_out = {Real(0),Real(0)};
+	get_pitau_m01(cos(th_rad),pi1,tau0,tau1,deg_max);
+	for (int n=0; n<deg_max; ++n) {
+		tc = tci/sqrt(Real((n+1)*(n+2)));
+		E_out[0] -= tc * ( ampl_out[3*n] * pi1[n] + ampl_out[3*n+1] * tau1[n] );
+		E_out[1] += tc * ( ampl_out[3*n] * tau1[n] + ampl_out[3*n+1] * pi1[n] );
+		ts += tc * ampl_out[3*n+2] * tau0[n];
+		tci *= -im1;
+	}
+	E_out[0] *= sqrt(2/_pi) * cos(ph_rad);
+	E_out[1] = sqrt(2/_pi) * ( sin(ph_rad)*E_out[1] + 0.5*ts );
+	return E_out;
+}
+
+// private members
+
+void sphere_ml::sphere_sm(Real kr, Complex eps_in, Complex eps_out, Complex *res) {
+	Complex t_besj_in, t_besdzj_in, t_besj_out, t_besdzj_out;
+	Complex t_besh_in, t_besdzh_in, t_besh_out, t_besdzh_out;
+	Complex kr_in, kr_out, te, tc, tcc;
+
+	kr_in = kr*sqrt(eps_in); if (arg(kr_in) < -1.e-8) kr_in = -kr_in;
+	kr_out = kr*sqrt(eps_out); if (arg(kr_out) < -1.e-8) kr_out = -kr_out;
+	te = eps_in/eps_out;
+
+	for (int n=0; n<deg_max; ++n) {
+		t_besj_in = besj(kr_in,n+1);
+		t_besdzj_in = bes_dzj(kr_in,n+1);
+		t_besj_out = besj(kr_out,n+1);
+		t_besdzj_out = bes_dzj(kr_out,n+1);
+
+		t_besh_in = besh1(kr_in,n+1);
+		t_besdzh_in = bes_dzh1(kr_in,n+1);
+		t_besh_out = besh1(kr_out,n+1);
+		t_besdzh_out = bes_dzh1(kr_out,n+1);
+
+		tc = Real(1) / (t_besj_in*t_besdzh_out - t_besh_out*t_besdzj_in);
+		res[0+8*n] = tc * (t_besh_out*t_besdzh_in - t_besh_in*t_besdzh_out); // 00e
+		res[2+8*n] = tc * im1 / kr_in; // 01e
+		res[6+8*n] = tc * (t_besj_out*t_besdzj_in - t_besj_in*t_besdzj_out); // 11e
+		res[4+8*n] = tc * im1 / kr_out; // 10e
+		tc = Real(1) / (te*t_besj_in*t_besdzh_out - t_besh_out*t_besdzj_in);
+		res[1+8*n] = tc * (t_besh_out*t_besdzh_in - te*t_besh_in*t_besdzh_out); // 00h
+		res[3+8*n] = tc * im1 / kr_out; // 01h
+		res[7+8*n] = tc * (t_besj_out*t_besdzj_in - te*t_besj_in*t_besdzj_out); // 11h
+		res[5+8*n] = tc * im1 * te / kr_in; // 10h
+	}
+}
+
+void sphere_ml::compose_sm(Complex *in, Complex *out, Complex *res) {
+	Complex tmp, tmp1;
+	for (int n=0; n<deg_max; ++n) {
+		tmp = Real(1) / (Real(1) - in[8*n+6]*out[8*n+0]);
+		tmp1 = out[8*n+2] * out[8*n+4];
+		res[8*n+0] = in[8*n+0] + tmp * in[8*n+2] * in[8*n+4] * out[8*n+0]; // 00e
+		res[8*n+2] = tmp * in[8*n+2] * out[8*n+2]; // 01e
+		res[8*n+4] = tmp * in[8*n+4] * out[8*n+4]; // 10e
+		res[8*n+6] = out[8*n+6] + tmp * tmp1 * in[8*n+6]; // 11e
+		
+		tmp = Real(1)/(Real(1) - in[8*n+7]*out[8*n+1]);
+		tmp1 = out[8*n+3] * out[8*n+5];
+		res[8*n+1] = in[8*n+1] + tmp * in[8*n+3] * in[8*n+5] * out[8*n+1]; // 00h
+		res[8*n+3] = tmp * in[8*n+3] * out[8*n+3]; // 01h
+		res[8*n+5] = tmp * in[8*n+5] * out[8*n+5]; // 10h
+		res[8*n+7] = out[8*n+7] + tmp * tmp1 * in[8*n+7]; // 11h
+	}
+}
+
+Real sphere_ml::directivity_dipole_z(Complex *ampl) const {
+	Real w_sca = 0;
+	Complex ti = -im1, tc;
+	Complex sum1, sum2, sum3;
+	
+	sum1 = sum2 = sum3 = Real(0);
+	
+	for (int n=0; n<deg_max; ++n) {
+		tc = ti/sqrt(Real((n+1)*(n+2)));
+		sum1 += tc * ( ampl[3*n]*pi1[n] + ampl[3*n+1]*tau1[n] );
+		sum2 += tc * ( ampl[3*n]*tau1[n] + ampl[3*n+1]*pi1[n] );
+		sum3 += tc * ampl[3*n+2]*tau0[n];
+		ti *= -im1;
+		w_sca += ( norm(ampl[3*n]) + norm(ampl[3*n+1]) ) + 0.5*norm(ampl[3*n+2]);
+	}
+	sum1 *= Real(2)*cos(phi);
+	sum2 *= Real(2)*im1*sin(phi);
+	
+//	return Real(2)*((sum1+sum3)*conj(sum1+sum3) + sum2*conj(sum2)).real() / w_sca;
+	return (norm(sum1+sum3) + norm(sum2)) / w_sca;
+}
+
+//############################################################################################
+
+void sphere_ml_parallel::set_directivity_angles(Real theta_, Real phi_) {
+	theta = theta_; phi = phi_;
+	get_pitau_m01(cos(theta),pi1,tau0,tau1,deg_max);
+}
+
+Real sphere_ml_parallel::get_directivity_edipole_x(Real kz, const std::vector<Real> &kR, const std::vector<Complex> &eps) {
+	int ni = size(kR);
+	if (size(eps)-1 != ni) {std::cout << "get_directivity_edipole_x error\n"; return Real(0);}
+	
+	Real D = 0.;
+	Complex kzn;
+	Complex rtml_in[8*_degree], rtml_out[8*_degree], rt[8*_degree];
+	Complex ampl_in[3*_degree], ampl_out[3*_degree];
+	
+	if (kz <= kR[0]) {
+		// s-matrix of a the multilayer
+		sphere_sm(kR[0], eps[0], eps[1], rtml_out);
+		for (int k=1; k<ni; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_out, rt, rtml_out);
+		}
+		// dipole amplitudes
+		kzn = kz*sqrt(eps[0]);
+		if (arg(kzn) < -1e-14) kzn = -kzn;
+		get_edipole_amp_x(kzn, ampl_out, 1);
+		// output amplitudes
+		for (int n=0; n<deg_max; ++n) {
+			ampl_out[3*n] *= rtml_out[8*n+2]; // e, m = +-1
+			ampl_out[3*n+1] *= rtml_out[8*n+3]; // h, m = +-1
+			ampl_out[3*n+2] *= rtml_out[8*n+3]; // h, m = 0
+		}
+	}
+	else if (kz <= kR[kR.size()-1]) { // the dipole is located in a layer
+																		// determine a dipole layer
+		int kd = 0;
+		while (kz > kR[++kd]); // determine a first interface with radius larger than the dipole position
+													 // calculate s-matrices of multilayers inside and outside the dipole position:
+		sphere_sm(kR[0], eps[0], eps[1], rtml_in);
+		for (int k=1; k<kd; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_in, rt, rtml_in);
+		}
+		sphere_sm(kR[kd], eps[kd], eps[kd+1], rtml_out);
+		for (int k=kd+1; k<ni; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_out, rt, rtml_out);
+		}
+		// dipole amplitudes
+		kzn = kz*sqrt(eps[kd]);
+		if (arg(kzn) < -1e-14) kzn = -kzn;
+		get_edipole_amp_x(kzn, ampl_in, 0);
+		get_edipole_amp_x(kzn, ampl_out, 1);
+		// output amplitudes
+		for (int n=0; n<deg_max; ++n) {
+			ampl_out[3*n] = ( ampl_in[3*n]*rtml_in[8*n+6] + ampl_out[3*n] ) * rtml_out[8*n+2]
+				/ (Real(1) - rtml_in[8*n+6]*rtml_out[8*n]); // e, m = 1
+			ampl_out[3*n+1] = ( ampl_in[3*n+1]*rtml_in[8*n+7] + ampl_out[3*n+1] ) * rtml_out[8*n+3]
+				/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 1
+			ampl_out[3*n+2] = ( ampl_in[3*n+2]*rtml_in[8*n+7] + ampl_out[3*n+2] ) * rtml_out[8*n+3]
+				/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 0
+		}
+	}
+	else { // the dipole is located outside the multilayer
+				 // s-matrix of a the multilayer
+		sphere_sm(kR[0], eps[0], eps[1], rtml_out);
+		for (int k=1; k<ni; ++k) {
+			sphere_sm(kR[k], eps[k], eps[k+1], rt);
+			compose_sm(rtml_out, rt, rtml_out);
+		}
+		// dipole amplitudes
+		kzn = kz*sqrt(eps[ni]);
+		if (arg(kzn) < -1e-14) kzn = -kzn;
+		get_edipole_amp_x(kzn, ampl_in, 0);
+		get_edipole_amp_x(kzn, ampl_out, 1);
+		// output amplitudes
+		for (int n=0; n<deg_max; ++n) {
+			ampl_out[3*n] += ampl_in[3*n] * rtml_out[8*n+6]; // e, m = +-1
+			ampl_out[3*n+1] += ampl_in[3*n+1] * rtml_out[8*n+7]; // h, m = +-1
+			ampl_out[3*n+2] += ampl_in[3*n+2] * rtml_out[8*n+7]; // h, m = 0
+		}
+	}
+
+	return directivity_dipole_z(ampl_out);
+}
+
+// private members
+
+void sphere_ml_parallel::get_edipole_amp_x(Complex kz, Complex *ampl, bool out) {
+	if (!out) {
+		Real tv = -0.25/sqrt(_pi), tvn, tv1, tv2;
+		Complex zf, zfd;
+		memset(ampl,0,3*deg_max*sizeof(Complex));
+		for (int n=0; n<deg_max; ++n) {
+			tvn = tv*sqrt(Real(2*n+3));
+			zf = besh1(kz,n+1); zfd = besh1d(kz,n+1);
+			ampl[3*n] = tvn*zf; // e-polarization, m = 1
+			ampl[3*n+1] = -im1*tvn*(zfd + zf/kz); // h-polarization m = 1
+			tv = -tv;
+		}
+	}
+	else {
+		Real tv = -0.25/sqrt(_pi), tvn, tv1, tv2;
+		Complex zf, zfd;
+		memset(ampl,0,3*deg_max*sizeof(Complex));
+		for (int n=0; n<deg_max; ++n) {
+			tvn = tv*sqrt(Real(2*n+3));
+			zf = besj(kz,n+1); zfd = besjd(kz,n+1);
+			ampl[3*n] = tvn*zf; // m = -1,1
+			ampl[3*n+1] = -im1*tvn*(zfd + zf/kz);
+			tv = -tv;
+		}
+	}
+}
+
+void sphere_ml_parallel::sphere_sm(Real kr, Complex eps_in, Complex eps_out, Complex *res) {
+	Complex t_besj_in, t_besdzj_in, t_besj_out, t_besdzj_out;
+	Complex t_besh_in, t_besdzh_in, t_besh_out, t_besdzh_out;
+	Complex kr_in, kr_out, te, tc, tcc;
+
+	kr_in = kr*sqrt(eps_in); if (arg(kr_in) < -1.e-8) kr_in = -kr_in;
+	kr_out = kr*sqrt(eps_out); if (arg(kr_out) < -1.e-8) kr_out = -kr_out;
+	te = eps_in/eps_out;
+
+	for (int n=0; n<deg_max; ++n) {
+		t_besj_in = besj(kr_in,n+1);
+		t_besdzj_in = bes_dzj(kr_in,n+1);
+		t_besj_out = besj(kr_out,n+1);
+		t_besdzj_out = bes_dzj(kr_out,n+1);
+
+		t_besh_in = besh1(kr_in,n+1);
+		t_besdzh_in = bes_dzh1(kr_in,n+1);
+		t_besh_out = besh1(kr_out,n+1);
+		t_besdzh_out = bes_dzh1(kr_out,n+1);
+
+		tc = Real(1) / (t_besj_in*t_besdzh_out - t_besh_out*t_besdzj_in);
+		res[0+8*n] = tc * (t_besh_out*t_besdzh_in - t_besh_in*t_besdzh_out); // 00e
+		res[2+8*n] = tc * im1 / kr_in; // 01e
+		res[6+8*n] = tc * (t_besj_out*t_besdzj_in - t_besj_in*t_besdzj_out); // 11e
+		res[4+8*n] = tc * im1 / kr_out; // 10e
+		tc = Real(1) / (te*t_besj_in*t_besdzh_out - t_besh_out*t_besdzj_in);
+		res[1+8*n] = tc * (t_besh_out*t_besdzh_in - te*t_besh_in*t_besdzh_out); // 00h
+		res[3+8*n] = tc * im1 / kr_out; // 01h
+		res[7+8*n] = tc * (t_besj_out*t_besdzj_in - te*t_besj_in*t_besdzj_out); // 11h
+		res[5+8*n] = tc * im1 * te / kr_in; // 10h
+	}
+}
+
+void sphere_ml_parallel::compose_sm(Complex *in, Complex *out, Complex *res) {
+	Complex tmp, tmp1;
+	for (int n=0; n<deg_max; ++n) {
+		tmp = Real(1) / (Real(1) - in[8*n+6]*out[8*n+0]);
+		tmp1 = out[8*n+2] * out[8*n+4];
+		res[8*n+0] = in[8*n+0] + tmp * in[8*n+2] * in[8*n+4] * out[8*n+0]; // 00e
+		res[8*n+2] = tmp * in[8*n+2] * out[8*n+2]; // 01e
+		res[8*n+4] = tmp * in[8*n+4] * out[8*n+4]; // 10e
+		res[8*n+6] = out[8*n+6] + tmp * tmp1 * in[8*n+6]; // 11e
+
+		tmp = Real(1)/(Real(1) - in[8*n+7]*out[8*n+1]);
+		tmp1 = out[8*n+3] * out[8*n+5];
+		res[8*n+1] = in[8*n+1] + tmp * in[8*n+3] * in[8*n+5] * out[8*n+1]; // 00h
+		res[8*n+3] = tmp * in[8*n+3] * out[8*n+3]; // 01h
+		res[8*n+5] = tmp * in[8*n+5] * out[8*n+5]; // 10h
+		res[8*n+7] = out[8*n+7] + tmp * tmp1 * in[8*n+7]; // 11h
+	}
+}
+
+Real sphere_ml_parallel::directivity_dipole_z(Complex *ampl) const {
+	Real w_sca = 0;
+	Complex ti = -im1, tc;
+	Complex sum1, sum2, sum3;
+
+	sum1 = sum2 = sum3 = Real(0);
+
+	for (int n=0; n<deg_max; ++n) {
+		tc = ti/sqrt(Real((n+1)*(n+2)));
+		sum1 += tc * ( ampl[3*n]*pi1[n] + ampl[3*n+1]*tau1[n] );
+		sum2 += tc * ( ampl[3*n]*tau1[n] + ampl[3*n+1]*pi1[n] );
+		sum3 += tc * ampl[3*n+2]*tau0[n];
+		ti *= -im1;
+		w_sca += ( norm(ampl[3*n]) + norm(ampl[3*n+1]) ) + 0.5*norm(ampl[3*n+2]);
+	}
+	sum1 *= Real(2)*cos(phi);
+	sum2 *= Real(2)*im1*sin(phi);
+
+	return (norm(sum1+sum3) + norm(sum2)) / w_sca;
+}

superdirectivity.h

@@ -0,0 +1,151 @@
+
+#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 = 15;
+
+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

superdirectivity.sln

@@ -0,0 +1,31 @@
+
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio Version 16
+VisualStudioVersion = 16.0.29709.97
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "superdirectivity", "superdirectivity.vcxproj", "{3200BA27-FDBD-444C-971B-D359DFBEC711}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|x64 = Debug|x64
+		Debug|x86 = Debug|x86
+		Release|x64 = Release|x64
+		Release|x86 = Release|x86
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Debug|x64.ActiveCfg = Debug|x64
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Debug|x64.Build.0 = Debug|x64
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Debug|x86.ActiveCfg = Debug|Win32
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Debug|x86.Build.0 = Debug|Win32
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Release|x64.ActiveCfg = Release|x64
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Release|x64.Build.0 = Release|x64
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Release|x86.ActiveCfg = Release|Win32
+		{3200BA27-FDBD-444C-971B-D359DFBEC711}.Release|x86.Build.0 = Release|Win32
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+	GlobalSection(ExtensibilityGlobals) = postSolution
+		SolutionGuid = {46A7CD6F-F439-4DC8-983E-4E1CFC2A5B45}
+	EndGlobalSection
+EndGlobal

superdirectivity.vcxproj

@@ -0,0 +1,160 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup Label="ProjectConfigurations">
+    <ProjectConfiguration Include="Debug|Win32">
+      <Configuration>Debug</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|Win32">
+      <Configuration>Release</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Debug|x64">
+      <Configuration>Debug</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|x64">
+      <Configuration>Release</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+  </ItemGroup>
+  <PropertyGroup Label="Globals">
+    <VCProjectVersion>16.0</VCProjectVersion>
+    <ProjectGuid>{3200BA27-FDBD-444C-971B-D359DFBEC711}</ProjectGuid>
+    <RootNamespace>superdirectivity</RootNamespace>
+    <WindowsTargetPlatformVersion>10.0</WindowsTargetPlatformVersion>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v142</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v142</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v142</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v142</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
+  <ImportGroup Label="ExtensionSettings">
+  </ImportGroup>
+  <ImportGroup Label="Shared">
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <PropertyGroup Label="UserMacros" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <LinkIncremental>true</LinkIncremental>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <LinkIncremental>true</LinkIncremental>
+    <IncludePath>C:\Program Files %28x86%29\Microsoft SDKs\MPI\Include;$(IncludePath)</IncludePath>
+    <LibraryPath>C:\Program Files %28x86%29\Microsoft SDKs\MPI\Lib\x64;$(LibraryPath)</LibraryPath>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <LinkIncremental>false</LinkIncremental>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <LinkIncremental>false</LinkIncremental>
+    <IncludePath>C:\Program Files %28x86%29\Microsoft SDKs\MPI\Include;$(IncludePath)</IncludePath>
+    <LibraryPath>C:\Program Files %28x86%29\Microsoft SDKs\MPI\Lib\x64;$(LibraryPath)</LibraryPath>
+  </PropertyGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalDependencies>msmpi.lib;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalDependencies>msmpi.lib;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemGroup>
+    <ClInclude Include="jade.h" />
+    <ClInclude Include="sph_bessel.h" />
+    <ClInclude Include="superdirectivity.h" />
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="jade.cpp" />
+    <ClCompile Include="main.cpp" />
+    <ClCompile Include="sph_bessel.cpp" />
+    <ClCompile Include="superdirectivity.cpp" />
+  </ItemGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
+  <ImportGroup Label="ExtensionTargets">
+  </ImportGroup>
+</Project>

superdirectivity.vcxproj.filters

@@ -0,0 +1,42 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup>
+    <Filter Include="Исходные файлы">
+      <UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
+      <Extensions>cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
+    </Filter>
+    <Filter Include="Файлы заголовков">
+      <UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
+      <Extensions>h;hh;hpp;hxx;hm;inl;inc;ipp;xsd</Extensions>
+    </Filter>
+    <Filter Include="Файлы ресурсов">
+      <UniqueIdentifier>{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}</UniqueIdentifier>
+      <Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
+    </Filter>
+  </ItemGroup>
+  <ItemGroup>
+    <ClInclude Include="superdirectivity.h">
+      <Filter>Файлы заголовков</Filter>
+    </ClInclude>
+    <ClInclude Include="sph_bessel.h">
+      <Filter>Файлы заголовков</Filter>
+    </ClInclude>
+    <ClInclude Include="jade.h">
+      <Filter>Файлы заголовков</Filter>
+    </ClInclude>
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="superdirectivity.cpp">
+      <Filter>Исходные файлы</Filter>
+    </ClCompile>
+    <ClCompile Include="sph_bessel.cpp">
+      <Filter>Исходные файлы</Filter>
+    </ClCompile>
+    <ClCompile Include="main.cpp">
+      <Filter>Исходные файлы</Filter>
+    </ClCompile>
+    <ClCompile Include="jade.cpp">
+      <Filter>Исходные файлы</Filter>
+    </ClCompile>
+  </ItemGroup>
+</Project>