Functions/objfuns/rad_objfun.cpp

   1 #include <cmath>
   2
   3 #include "rad_objfun.h"
   4
   5 using namespace std;
   6
   7 /*Continuous radiator function from Jose study*/
   8
   9 extern "C"
  10   {
  11     double real_radiator_extern_(double *, const int &);
  12   }
  13
  14 double radiator( const vector<double> &x )
  15 {
  16         const int n = x.size();
  17         double J_result=0;
  18         const int number_layers=C_NUMBER_LAYERS;
  19         double* array = new double[number_layers*2];
  20         vector <double> structure;
  21
  22
  23         structure = var2struct(x);
  24
  25         for (unsigned int i=0;i<structure.size();i++){
  26             array[i] = structure[i];
  27         }
  28
  29         J_result = real_radiator_extern_(array, number_layers);
  30
  31         delete[] array;
  32
  33             return J_result;
  34 }
  35
  36 vector<double> var2struct(const vector<double> &x){
  37         const int n = x.size();
  38         int number_layers=C_NUMBER_LAYERS;
  39         int number_active_layers=C_NUMBER_ACTIVE_LAYERS;
  40         int layer_block,rest_layers,number_blocks;
  41         vector <double> array(2*number_layers,0);
  42
  43         number_blocks=(number_active_layers+1);
  44         layer_block=(number_layers-number_active_layers)/number_blocks;
  45         rest_layers=number_layers-(layer_block*number_blocks+number_active_layers);
  46         if(rest_layers <0) {
  47             cout << "Error assigning values to array\n";
  48             exit(1);
  49         }
  50
  51         int active_layer=0; // Active_layer and j have always the same values, but it is easier to follow the algorithm in this way
  52         int k=0;
  53         for (int j=0; j<number_blocks;j++){
  54
  55         for (int i=j*layer_block+active_layer;i<(j+1)*layer_block+active_layer;i++){
  56             array[i] = x[k];
  57             k++;
  58         }
  59         if(active_layer < number_active_layers) array[(j+1)*layer_block+active_layer]=-1.0;
  60         active_layer++;
  61         }
  62         //Rest of layers
  63         if (rest_layers != 0) {
  64             for (int i=number_blocks*layer_block+number_active_layers; i<number_blocks*layer_block+number_active_layers+rest_layers; i++){
  65                 array[i] = x[k];
  66                 k++;
  67
  68             }
  69         }
  70
  71         for (int i=number_layers;i<2*number_layers;i++){
  72             array[i] = x[k];
  73             k++;
  74         }
  75
  76         return array;
  77 }