scatt/src/nmie.hpp

#ifndef SRC_NMIE_HPP_
#define SRC_NMIE_HPP_
//**********************************************************************************//
//    Copyright (C) 2009-2018  Ovidio Pena <ovidio@bytesfall.com>                   //
//    Copyright (C) 2013-2018  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 at least one of the following references:                      //
//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
//        a multilayered sphere," Computer Physics Communications,                  //
//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
//    [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie              //
//        calculation of electromagnetic near-field for a multilayered              //
//        sphere," Computer Physics Communications, vol. 214, May 2017,             //
//        pp. 225-230.                                                              //
//                                                                                  //
//    You should have received a copy of the GNU General Public License             //
//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
//**********************************************************************************//

#define VERSION "2.2"  //Compare with Makefile and setup.py
#include <array>
#include <complex>
#include <cstdlib>
#include <iostream>
#include <vector>
#include "nmie-precision.hpp"
#ifdef MULTI_PRECISION
#include <boost/math/constants/constants.hpp>
#endif
namespace nmie {
  int ScattCoeffs(const unsigned int L, const int pl,
                  const std::vector<double> &x, const std::vector<std::complex<double> > &m,
                  const int nmax,
                  std::vector<std::complex<double> > &an,
                  std::vector<std::complex<double> > &bn);

  int ExpanCoeffs(const unsigned int L, const int pl,
                  const std::vector<double> &x, const std::vector<std::complex<double> > &m,
                  const int nmax,
                  std::vector<std::vector<std::complex<double> > > &an,
                  std::vector<std::vector<std::complex<double> > > &bn,
                  std::vector<std::vector<std::complex<double> > > &cn,
                  std::vector<std::vector<std::complex<double> > > &dn);

//helper functions
template <typename FloatType>
double eval_delta(const int steps, const double from_value, const double to_value);


template<class T> inline T pow2(const T value) {return value*value;}

template<class T> inline T cabs(const std::complex<T> value)
{return nmm::sqrt(pow2(value.real()) + pow2(value.imag()));}

template <typename FloatType>
int newround(FloatType x) {
  return x >= 0 ? static_cast<int>(x + 0.5):static_cast<int>(x - 0.5);
  //return x >= 0 ? (x + 0.5).convert_to<int>():(x - 0.5).convert_to<int>();
}
template<typename T>
inline std::complex<T> my_exp(const std::complex<T> &x) {
  using std::exp; // use ADL
  T const &r = exp(x.real());
  return std::polar(r, x.imag());
}

// pl, nmax, mode_n, mode_type
    int nMie(const unsigned int L,
           const int pl,
           std::vector<double> &x, std::vector<std::complex<double> > &m,
           const unsigned int nTheta, std::vector<double> &Theta,
           const int nmax,
           double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
           double *g, double *Albedo,
           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2,
           int mode_n, int mode_type);
  // pl and nmax
    int nMie(const unsigned int L,
           const int pl,
           std::vector<double> &x, std::vector<std::complex<double> > &m,
           const unsigned int nTheta, std::vector<double> &Theta,
           const int nmax,
           double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
           double *g, double *Albedo,
           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
  // no pl and nmax
  int nMie(const unsigned int L,
           std::vector<double> &x, std::vector<std::complex<double> > &m,
           const unsigned int nTheta, std::vector<double> &Theta,
           double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
           double *g, double *Albedo,
           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
  // pl
  int nMie(const unsigned int L,
           const int pl,
           std::vector<double> &x, std::vector<std::complex<double> > &m,
           const unsigned int nTheta, std::vector<double> &Theta,
           double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
           double *g, double *Albedo,
           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
  // nmax
  int nMie(const unsigned int L,
           std::vector<double> &x, std::vector<std::complex<double> > &m,
           const unsigned int nTheta, std::vector<double> &Theta,
           const int nmax,
           double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr,
           double *g, double *Albedo,
           std::vector<std::complex<double> > &S1, std::vector<std::complex<double> > &S2);
  int nField(const unsigned int L, const int pl,
             const std::vector<double> &x, const std::vector<std::complex<double> > &m, const int nmax,
             const int mode_n, const int mode_type, const unsigned int ncoord,
             const std::vector<double> &Xp, const std::vector<double> &Yp, const std::vector<double> &Zp,
             std::vector<std::vector<std::complex<double> > > &E, std::vector<std::vector<std::complex<double> > > &H);

  // constants for per mode evaluation
  enum Modes {kAll = -1, kElectric = 0, kMagnetic = 1};

  template <typename FloatType = double>
  class MultiLayerMie {
   public:
    const FloatType PI_=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354886230577456498035593634568174324112515076069479451096596094025228879710893145669136867228748940560101503308;
    // light speed [m s-1]
    const double cc_ = 2.99792458e8;
    // assume non-magnetic (MU=MU0=const) [N A-2]
    const FloatType mu_ = 4.0*PI_*1.0e-7;
    // Run calculation
    void RunMieCalculation();
    void RunFieldCalculation();
    void RunFieldCalculationPolar(const int outer_arc_points = 1,
                                  const int radius_points=1,
                                  const double from_Rho=0, const double to_Rho=static_cast<double>(1.),
                                  const double from_Theta=0, const double to_Theta=static_cast<double>(3.14159265358979323),
                                  const double from_Phi=0, const double to_Phi=static_cast<double>(3.14159265358979323),
                                  const bool isIgnoreAvailableNmax = true); // TODO change to false for production

    void calcScattCoeffs();
    void calcExpanCoeffs();

    // Return calculation results
    FloatType GetQext();
    FloatType GetQsca();
    FloatType GetQabs();
    FloatType GetQbk();
    FloatType GetQpr();
    FloatType GetAsymmetryFactor();
    FloatType GetAlbedo();
    std::vector<std::complex<FloatType> > GetS1();
    std::vector<std::complex<FloatType> > GetS2();

    std::vector<std::complex<FloatType> > GetAn(){return an_;};
    std::vector<std::complex<FloatType> > GetBn(){return bn_;};

    std::vector< std::vector<std::complex<FloatType> > > GetLayerAn(){return aln_;};
    std::vector< std::vector<std::complex<FloatType> > > GetLayerBn(){return bln_;};
    std::vector< std::vector<std::complex<FloatType> > > GetLayerCn(){return cln_;};
    std::vector< std::vector<std::complex<FloatType> > > GetLayerDn(){return dln_;};

    // Problem definition
    // Modify size of all layers
    void SetLayersSize(const std::vector<FloatType> &layer_size);
    // Modify refractive index of all layers
    void SetLayersIndex(const std::vector< std::complex<FloatType> > &index);
    void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml, unsigned int &l);
    void GetIndexAtRadius(const FloatType Rho, std::complex<FloatType> &ml);
    // Modify scattering (theta) angles
    void SetAngles(const std::vector<FloatType> &angles);
    // Modify coordinates for field calculation
    void SetFieldCoords(const std::vector< std::vector<FloatType> > &coords);
    // Modify index of PEC layer
    void SetPECLayer(int layer_position = 0);
    // Modify the mode taking into account for evaluation of output variables
    void SetModeNmaxAndType(int mode_n, int mode_type){mode_n_ = mode_n; mode_type_ = mode_type;};

    // Set a fixed value for the maximun number of terms
    void SetMaxTerms(int nmax);
    // Get maximum number of terms
    int GetMaxTerms() {return nmax_;};

    bool isMieCalculated(){return isMieCalculated_;};
    // Clear layer information
    void ClearLayers();
    void MarkUncalculated();

    // Read parameters
    // Get total size parameter of particle
    FloatType GetSizeParameter();
    // Returns size of all layers
    std::vector<FloatType> GetLayersSize(){return size_param_;};
    // Returns refractive index of all layers
    std::vector<std::complex<FloatType> > GetLayersIndex(){return refractive_index_;};
    // Returns scattering (theta) angles
    std::vector<FloatType> GetAngles(){return theta_;};
    // Returns coordinates used for field calculation
    std::vector<std::vector<FloatType> > GetFieldCoords(){return coords_;};
    // Returns index of PEC layer
    int GetPECLayer(){return PEC_layer_position_;};

    std::vector<std::vector< std::complex<FloatType> > > GetFieldE(){return E_;};   // {X[], Y[], Z[]}
    std::vector<std::vector< std::complex<FloatType> > > GetFieldH(){return H_;};
    // Get fields in spherical coordinates.
    std::vector<std::vector< std::complex<FloatType> > > GetFieldEs(){return E_;};   // {rho[], teha[], phi[]}
    std::vector<std::vector< std::complex<FloatType> > > GetFieldHs(){return H_;};

  protected:
    // Size parameter for all layers
    std::vector<FloatType> size_param_;
    // Refractive index for all layers
    std::vector< std::complex<FloatType> > refractive_index_;
    // Scattering coefficients
    std::vector<std::complex<FloatType> > an_, bn_;
    std::vector< std::vector<std::complex<FloatType> > > aln_, bln_, cln_, dln_;
    // Points for field evaluation
    std::vector< std::vector<FloatType> > coords_;
    std::vector< std::vector<FloatType> > coords_polar_;

  private:
    unsigned int calcNstop(FloatType xL = -1);
    unsigned int calcNmax(FloatType xL = -1);

    std::complex<FloatType> calc_an(int n, FloatType XL, std::complex<FloatType> Ha, std::complex<FloatType> mL,
                                 std::complex<FloatType> PsiXL, std::complex<FloatType> ZetaXL,
                                 std::complex<FloatType> PsiXLM1, std::complex<FloatType> ZetaXLM1);
    std::complex<FloatType> calc_bn(int n, FloatType XL, std::complex<FloatType> Hb, std::complex<FloatType> mL,
                                 std::complex<FloatType> PsiXL, std::complex<FloatType> ZetaXL,
                                 std::complex<FloatType> PsiXLM1, std::complex<FloatType> ZetaXLM1);
    std::complex<FloatType> calc_S1(int n, std::complex<FloatType> an, std::complex<FloatType> bn,
                                 FloatType Pi, FloatType Tau);
    std::complex<FloatType> calc_S2(int n, std::complex<FloatType> an, std::complex<FloatType> bn,
                                 FloatType Pi, FloatType Tau);
    void calcD1D3(std::complex<FloatType> z,
                  std::vector<std::complex<FloatType> > &D1,
                  std::vector<std::complex<FloatType> > &D3);
    void calcPsiZeta(std::complex<FloatType> x,
                     std::vector<std::complex<FloatType> > &Psi,
                     std::vector<std::complex<FloatType> > &Zeta);
    void calcPiTau(const FloatType &costheta,
                   std::vector<FloatType> &Pi, std::vector<FloatType> &Tau);
    void calcSpherHarm(const std::complex<FloatType> Rho, const FloatType Theta, const FloatType Phi,
                       const std::complex<FloatType> &rn, const std::complex<FloatType> &Dn,
                       const FloatType &Pi, const FloatType &Tau, const FloatType &n,
                       std::vector<std::complex<FloatType> > &Mo1n, std::vector<std::complex<FloatType> > &Me1n,
                       std::vector<std::complex<FloatType> > &No1n, std::vector<std::complex<FloatType> > &Ne1n);

    void calcFieldByComponents(const FloatType Rho, const FloatType Theta, const FloatType Phi,
                               const std::vector<std::complex<FloatType> > &Psi,
                               const std::vector<std::complex<FloatType> > &D1n,
                               const std::vector<std::complex<FloatType> > &Zeta,
                               const std::vector<std::complex<FloatType> > &D3n,
                               const std::vector<FloatType> &Pi,
                               const std::vector<FloatType> &Tau,
                               std::vector<std::complex<FloatType> > &E,
                               std::vector<std::complex<FloatType> > &H);

    bool isExpCoeffsCalc_ = false;
    bool isScaCoeffsCalc_ = false;
    bool isMieCalculated_ = false;

    // Scattering angles for scattering pattern in radians
    std::vector<FloatType> theta_;
    // Should be -1 if there is no PEC.
    int PEC_layer_position_ = -1;

    int mode_n_ = Modes::kAll;
    int mode_type_ = Modes::kAll;

    // with calcNmax(int first_layer);
    int nmax_ = -1;
    int nmax_preset_ = -1;
    int available_maximal_nmax_ = -1;
    /// Store result
    FloatType Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
    std::vector<std::vector< std::complex<FloatType> > > E_, H_;  // {X[], Y[], Z[]}
    std::vector<std::vector< std::complex<FloatType> > > Es_, Hs_;  // {X[], Y[], Z[]}
    std::vector<std::complex<FloatType> > S1_, S2_;
    void calcMieSeriesNeededToConverge(const FloatType Rho);
    void calcPiTauAllTheta(const double from_Theta,
                           const double to_Theta,
                           std::vector<std::vector<FloatType>> &Pi,
                           std::vector<std::vector<FloatType>> &Tau);
    void calcRadialOnlyDependantFunctions(const FloatType from_Rho,
                                          const FloatType to_Rho,
                                          const bool isIgnoreAvailableNmax,
                                          std::vector<std::vector<std::complex<FloatType>>> &Psi,
                                          std::vector<std::vector<std::complex<FloatType>>> &D1n,
                                          std::vector<std::vector<std::complex<FloatType>>> &Zeta,
                                          std::vector<std::vector<std::complex<FloatType>>> &D3n);
  };  // end of class MultiLayerMie

}  // end of namespace nmie
#endif  // SRC_NMIE_HPP_