//**********************************************************************************//
//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
//                                                                                  //
//    This file is part of scattnlay                                                //
//                                                                                  //
//    This program 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.                                           //
//                                                                                  //
//    This program 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.                                  //
//                                                                                  //
//    The only additional remark is that we expect that all publications            //
//    describing work using this software, or all commercial products               //
//    using it, cite the following reference:                                       //
//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
//        a multilayered sphere," Computer Physics Communications,                  //
//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
//                                                                                  //
//    You should have received a copy of the GNU General Public License             //
//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
//**********************************************************************************//
//   This program returns expansion coefficents of Mie series
#include "../src/nmie.hpp"
#include "../src/nmie.hpp"
#include "../src/nmie-precision.hpp"
#include <complex>
#include <cstdio>
#include <string>
#include "../src/nmie-applied.hpp"
#include "../src/nmie-applied-impl.hpp"
#include "../src/nmie-precision.hpp"
#include "./read-spectra.h"

// template<class T> inline T pow2(const T value) {return value*value;}
int main(int argc, char *argv[]) {
  using namespace nmie ;
  try {
    //read_spectra::ReadSpectra Si_index;
    //read_spectra::ReadSpectra plot_core_index_, plot_TiN_;
   // std::string core_filename("Si-int.txt");
    //std::string core_filename("Ag.txt");
    //std::string TiN_filename("TiN.txt");
    //std::string TiN_filename("Ag-int.txt");
    //std::string TiN_filename("Si.txt");
    //std::string shell_filename(core_filename);

    nmie::MultiLayerMieApplied<nmie::FloatType> multi_layer_mie;  
    const std::complex<double> epsilon_Si(16, 0);
    //const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
    const std::complex<double> index_Si = std::sqrt(epsilon_Si);
    //const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
    double WL=550; //nm
    //double core_width = 102; //nm Si // radius
    double core_width = 196/2.0; //nm Si // radius
    //double inner_width = 8.22; //nm Ag
    //double outer_width = 67.91; //nm  Si
    //bool isSiAgSi = true;
    //double delta_width = 25.0;
    // //bool isSiAgSi = false;
    // if (isSiAgSi) {
      // core_width = 5.27; //nm Si
      // inner_width = 8.22; //nm Ag
      // outer_width = 67.91; //nm  Si
    multi_layer_mie.AddTargetLayer(core_width, index_Si);
      // multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
      // multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
    // } else {
    //   inner_width = 31.93; //nm Ag
    //   outer_width = 4.06; //nm  Si
    //   multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
    //   multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
    // }

    // for (int dd = 0; dd<50; ++dd) {
    //   delta_width = dd;
    // FILE *fp;
    // std::string fname = "absorb-layered-spectra-d"+std::to_string(dd)+".dat";
    // fp = fopen(fname.c_str(), "w");

    multi_layer_mie.SetWavelength(WL);
    multi_layer_mie.RunMieCalculation();
    
    double Qsca = static_cast<double>(multi_layer_mie.GetQsca());
    printf("Qsca = %g\n", Qsca);//*3.14159*core_width*core_width*1e-6);
    std::vector< std::vector<std::complex<nmie::FloatType> > > aln, bln, cln, dln;
    multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);






 //      nmie::ConvertComplexVectorVector<double>(aln);
 //    std::string str = std::string("#WL ");
 //    for (int l = 0; l<d_aln.size(); ++l) {
 //      for (int n = 0; n<3; ++n) {
	// str+="|a|^2+|d|^2_ln"+std::to_string(l)+std::to_string(n)+" "
	//   + "|b|^2+|c|^2_ln"+std::to_string(l)+std::to_string(n)+" ";
 //      }
 //    // }
 //    str+="\n";
 //    fprintf(fp, "%s", str.c_str());
 //    fprintf(fp, "# |a|+|d|");
 //    str.clear();

 //    double from_WL = 400;
 //    double to_WL = 600;
 //    int total_points = 401;
 //    double delta_WL = std::abs(to_WL - from_WL)/(total_points-1.0);
 //    Si_index.ReadFromFile(core_filename).ResizeToComplex(from_WL, to_WL, total_points)
 //      .ToIndex();
 //    Ag_index.ReadFromFile(TiN_filename).ResizeToComplex(from_WL, to_WL, total_points)
 //      .ToIndex();
 //    auto Si_data = Si_index.GetIndex();
 //    auto Ag_data = Ag_index.GetIndex();
 //    for (int i=0; i < Si_data.size(); ++i) {
 //      const double& WL = Si_data[i].first;
 //      const std::complex<double>& Si = Si_data[i].second;
 //      const std::complex<double>& Ag = Ag_data[i].second;
 //      str+=std::to_string(WL);
 //      multi_layer_mie.ClearTarget();      
 //      if (isSiAgSi) {
	// multi_layer_mie.AddTargetLayer(core_width, Si);
	// multi_layer_mie.AddTargetLayer(inner_width, Ag);
	// multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
 //      } else {
	// inner_width = 31.93; //nm Ag
	// outer_width = 4.06; //nm  Si
	// multi_layer_mie.AddTargetLayer(inner_width, Ag);
	// multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
 //      }
 //      multi_layer_mie.SetWavelength(WL);
 //      multi_layer_mie.RunMieCalculation();
 //      multi_layer_mie.GetQabs();
 //      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
 //      for (int l = 0; l<aln.size(); ++l) {
	// for (int n = 0; n<3; ++n) {
	//   str+=" "+std::to_string(static_cast<double>(pow2(std::abs(aln[l][n]))+
 //                                                      pow2(std::abs(dln[l][n]))))
	//     + " "
	//     + std::to_string(static_cast<double>(pow2(std::abs(bln[l][n]))
 //                                                 + pow2(std::abs(cln[l][n])) ));

	//   // str+=" "+std::to_string(aln[l][n].real() - pow2(std::abs(aln[l][n]))
	//   // 			  +dln[l][n].real() - pow2(std::abs(dln[l][n])))
	//   //   + " "
	//   //   + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n]))
	//   // 			  +cln[l][n].real() - pow2(std::abs(cln[l][n])) );
	// }
 //      }
 //      str+="\n";
 //      fprintf(fp, "%s", str.c_str());
 //      str.clear();
 //    }

 //    fclose(fp);
    // }

    // WL = 500;
    // multi_layer_mie.SetWavelength(WL);
    // multi_layer_mie.RunMieCalculation();
    // multi_layer_mie.GetQabs();
    // multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);

    printf("\n Scattering");
    for (int l = 0; l<aln.size(); ++l) {
      int n = 0;
      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
      n = 1;
      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
       n = 2;
      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
             n = 3;
      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
             n = 4;
      printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
      printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
      printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
      printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
    }
    // printf("\n Absorbtion\n");
    // for (int l = 0; l<aln.size(); ++l) {
    //   if (l == aln.size()-1) printf(" Total ");
    //   printf("===== l=%i   =====\n", l);
    //   int n = 0;
    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
    //   printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
    //   n = 1;
    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
    //   printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
    //   // n = 2;
    //   // printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
    //   // printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
    //   // printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
    //   // printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
    // }


  } catch( const std::invalid_argument& ia ) {
    // Will catch if  multi_layer_mie fails or other errors.
    std::cerr << "Invalid argument: " << ia.what() << std::endl;
    return -1;
  }  
    return 0;
}