moved expansion coefficents declaration to the protected part of the class

examples/example-get-Mie.cc

@@ -30,7 +30,7 @@
 #include "../src/nmie-applied.h"
 int main(int argc, char *argv[]) {
   try {
-    nmie::MultiLayerMieApplied multi_layer_mie_;  
+    nmie::MultiLayerMieApplied 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);
@@ -42,19 +42,29 @@ int main(int argc, char *argv[]) {
     //bool isSiAgSi = true;
     bool isSiAgSi = false;
     if (isSiAgSi) {
-      multi_layer_mie_.AddTargetLayer(core_width, index_Si);
-      multi_layer_mie_.AddTargetLayer(inner_width, index_Ag);
-      multi_layer_mie_.AddTargetLayer(outer_width, index_Si);
+      multi_layer_mie.AddTargetLayer(core_width, index_Si);
+      multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
+      multi_layer_mie.AddTargetLayer(outer_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, index_Si);
+      multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
+      multi_layer_mie.AddTargetLayer(outer_width, index_Si);
     }
-    multi_layer_mie_.SetWavelength(WL);
-    multi_layer_mie_.RunMieCalculation();
-    double Qabs = multi_layer_mie_.GetQabs();
-    printf("Qabs = %g", Qabs);
+    multi_layer_mie.SetWavelength(WL);
+    multi_layer_mie.RunMieCalculation();
+    double Qabs = multi_layer_mie.GetQabs();
+    printf("Qabs = %g\n", Qabs);
+    std::vector< std::vector<std::complex<double> > > aln, bln, cln, dln;
+    multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+    for (int l = 0; l<aln.size(); ++l) {
+    int n = 0;
+    printf("aln[%i][%i] = %g, %gi)\n", l, n, aln[l][n].real(), aln[l][n].imag());
+    printf("bln[%i][%i] = %g, %gi)\n", l, n, bln[l][n].real(), bln[l][n].imag());
+    printf("cln[%i][%i] = %g, %gi)\n", l, n, cln[l][n].real(), cln[l][n].imag());
+    printf("dln[%i][%i] = %g, %gi)\n", l, n, dln[l][n].real(), dln[l][n].imag());
+  }
+
 
 
   } catch( const std::invalid_argument& ia ) {

src/nmie-applied.cc

@@ -444,4 +444,20 @@ c    MM + 1  and - 1, alternately
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void MultiLayerMieApplied::GetExpanCoeffs( std::vector< std::vector<std::complex<double> > >& aln, std::vector< std::vector<std::complex<double> > >& bln, std::vector< std::vector<std::complex<double> > >& cln, std::vector< std::vector<std::complex<double> > >& dln) {
+    ConvertToSP();  // Case of call before running full Mie calculation.
+    // Calculate scattering coefficients an_ and bn_
+    calcScattCoeffs();
+    // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
+    calcExpanCoeffs();
+    aln = aln_;
+    bln = bln_;
+    cln = cln_;
+    dln = dln_;
+    
+  }  // end of void MultiLayerMieApplied::GetExpanCoeffs( ...)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+
 }  // end of namespace nmie

src/nmie-applied.h

@@ -120,7 +120,13 @@ namespace nmie {
     std::vector<double> GetPatternEkSP();
     std::vector<double> GetPatternHkSP();
     std::vector<double> GetPatternUnpolarizedSP();
-    
+
+    void GetExpanCoeffs
+      (std::vector< std::vector<std::complex<double> > >& aln,
+       std::vector< std::vector<std::complex<double> > >& bln,
+       std::vector< std::vector<std::complex<double> > >& cln,
+       std::vector< std::vector<std::complex<double> > >& dln);
+
 
     // Output results (data file + python script to plot it with matplotlib)
     void PlotSpectra();

src/nmie.cc

@@ -1045,7 +1045,7 @@ namespace nmie {
   //**********************************************************************************//
   void MultiLayerMie::calcExpanCoeffs() {
     if (!isScaCoeffsCalc_)
-      throw std::invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
+      throw std::invalid_argument("(calcExpanCoeffs) You should calculate external coefficients first!");
 
     isExpCoeffsCalc_ = false;
 

src/nmie.h

@@ -113,7 +113,10 @@ namespace nmie {
     std::vector<double> size_param_;
     // Refractive index for all layers
     std::vector< std::complex<double> > refractive_index_;
-    // Scattering angles for scattering pattern in radians
+    // Scattering coefficients
+    std::vector<std::complex<double> > an_, bn_;
+    std::vector< std::vector<std::complex<double> > > aln_, bln_, cln_, dln_;
+    void calcExpanCoeffs();
 
   private:
     void calcNstop();
@@ -142,7 +145,6 @@ namespace nmie {
                        const double& Pi, const double& Tau, const double& n,
                        std::vector<std::complex<double> >& Mo1n, std::vector<std::complex<double> >& Me1n, 
                        std::vector<std::complex<double> >& No1n, std::vector<std::complex<double> >& Ne1n);
-    void calcExpanCoeffs();
 
     void calcField(const double Rho, const double Theta, const double Phi,
                    std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
@@ -151,6 +153,7 @@ namespace nmie {
     bool isScaCoeffsCalc_ = false;
     bool isMieCalculated_ = false;
 
+    // Scattering angles for scattering pattern in radians
     std::vector<double> theta_;
     // Should be -1 if there is no PEC.
     int PEC_layer_position_ = -1;
@@ -158,10 +161,7 @@ namespace nmie {
     // with calcNmax(int first_layer);
     int nmax_ = -1;
     int nmax_preset_ = -1;
-    // Scattering coefficients
-    std::vector<std::complex<double> > an_, bn_;
     std::vector< std::vector<double> > coords_;
-    std::vector< std::vector<std::complex<double> > > aln_, bln_, cln_, dln_;
     /// Store result
     double 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<double> > > E_, H_;  // {X[], Y[], Z[]}