//**********************************************************************************//
//    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 <complex>
#include <cstdio>
#include <string>
#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, Ag_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(18.4631066585, 0.6259727805);
    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 = 500;             // nm
    double core_width = 5.27;    // nm Si
    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 Qabs = static_cast<double>(multi_layer_mie.GetQabs());
      printf("Qabs = %g\n", Qabs);
      std::vector<std::vector<std::complex<nmie::FloatType> > > aln, bln, cln,
          dln;
      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
      std::vector<std::vector<std::complex<double> > > d_aln =
          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());
    // }
    // 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;
}