fixed bug if components with different grain sizes were used includin…

…g gaps in the distribution
polaris-MCRT · Sep 21, 2023 · 844c606 · 844c606
1 parent a136607
commit 844c606
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 18 deletions.
diff --git a/src/Dust.cpp b/src/Dust.cpp
@@ -2366,7 +2366,7 @@ bool CDustComponent::writeComponent(string path_data, string path_plot)
         // Init text file writer for scattering matrix
         ofstream scat_writer(path_scat.c_str());
 
-        uint nr_of_scat_theta = 2 * NANG - 1;
+        uint nr_of_scat_theta_tmp = 2 * NANG - 1;
 
         // Error message if the write does not work
         if(scat_writer.fail())
@@ -2391,9 +2391,9 @@ bool CDustComponent::writeComponent(string path_data, string path_plot)
         for(uint w = wavelength_offset; w < nr_of_wavelength; w++)
         {
             uint wID = w - wavelength_offset;
-            S11[wID] = new double[nr_of_scat_theta];
-            S12[wID] = new double[nr_of_scat_theta];
-            for(uint sth = 0; sth < nr_of_scat_theta; sth++)
+            S11[wID] = new double[nr_of_scat_theta_tmp];
+            S12[wID] = new double[nr_of_scat_theta_tmp];
+            for(uint sth = 0; sth < nr_of_scat_theta_tmp; sth++)
             {
                 // Init and reset variables
                 double sum = 0;
@@ -2402,7 +2402,7 @@ bool CDustComponent::writeComponent(string path_data, string path_plot)
                 double * S12_tmp = new double[nr_of_dust_species];
                 for(uint a = 0; a < nr_of_dust_species; a++)
                 {
-                    if(sizeIndexUsed(a))
+                    if(sizeIndexUsed(a) && nr_of_scat_theta[a][w] != 0)
                     {
                         double Csca_tmp = getCscaMean(a, w);
                         sum += Csca_tmp;
@@ -2446,7 +2446,7 @@ bool CDustComponent::writeComponent(string path_data, string path_plot)
         scat_writer << "set multiplot layout 2,1 rowsfirst" << endl;
 
         if(nr_of_wavelength > 1)
-            scat_writer << "set xrange[" << 0 << ":" << nr_of_scat_theta << "]" << endl;
+            scat_writer << "set xrange[" << 0 << ":" << nr_of_scat_theta_tmp << "]" << endl;
         scat_writer << "set yrange[" << S11min << ":" << S11max << "]" << endl;
         scat_writer << "set format x \"%.1f\"" << endl;
         scat_writer << "set format y \"%.1te%02T\"" << endl;
@@ -2467,13 +2467,13 @@ bool CDustComponent::writeComponent(string path_data, string path_plot)
         {
             uint wID = w - wavelength_offset;
 
-            for(uint sth = 0; sth < nr_of_scat_theta; sth++)
+            for(uint sth = 0; sth < nr_of_scat_theta_tmp; sth++)
                 scat_writer << sth << "\t" << S11[wID][sth] << endl;
             scat_writer << "e" << endl;
         }
 
         if(nr_of_wavelength > 1)
-            scat_writer << "set xrange[" << 0 << ":" << nr_of_scat_theta << "]" << endl;
+            scat_writer << "set xrange[" << 0 << ":" << nr_of_scat_theta_tmp << "]" << endl;
         scat_writer << "set yrange[" << S12min << ":" << S12max << "]" << endl;
         scat_writer << "set format x \"%.1f\"" << endl;
         scat_writer << "set format y \"%.1te%02T\"" << endl;
@@ -2494,7 +2494,7 @@ bool CDustComponent::writeComponent(string path_data, string path_plot)
         {
             uint wID = w - wavelength_offset;
 
-            for(uint sth = 0; sth < nr_of_scat_theta; sth++)
+            for(uint sth = 0; sth < nr_of_scat_theta_tmp; sth++)
                 scat_writer << sth << "\t" << S12[wID][sth] << endl;
             scat_writer << "e" << endl;
         }

diff --git a/src/Dust.h b/src/Dust.h
@@ -205,6 +205,9 @@ class CDustComponent
                 delete[] enthalpy[a];
             delete[] enthalpy;
         }
+        if(sca_mat != 0)
+            cleanScatteringData();
+
         if(nr_of_scat_theta != 0)
         {
             for(uint a = 0; a < nr_of_dust_species; a++)
@@ -280,9 +283,6 @@ class CDustComponent
         if(abs_prob != 0)
             delete[] abs_prob;
 
-        if(sca_mat != 0)
-            cleanScatteringData();
-
         if(a_eff != 0)
             delete[] a_eff;
         if(a_eff_1_5 != 0)
@@ -1118,13 +1118,16 @@ class CDustComponent
                 {
                     for(uint w = 0; w < nr_of_wavelength; w++)
                     {
-                        for(uint inc = 0; inc < nr_of_incident_angles; inc++)
+                        if(nr_of_scat_theta[a][w] != 0)
                         {
-                            for(uint sph = 0; sph < nr_of_scat_phi; sph++)
-                                delete[] sca_mat[a][w][inc][sph];
-                            delete[] sca_mat[a][w][inc];
+                            for(uint inc = 0; inc < nr_of_incident_angles; inc++)
+                            {
+                                for(uint sph = 0; sph < nr_of_scat_phi; sph++)
+                                    delete[] sca_mat[a][w][inc][sph];
+                                delete[] sca_mat[a][w][inc];
+                            }
+                            delete[] sca_mat[a][w];
                         }
-                        delete[] sca_mat[a][w];
                     }
                     delete[] sca_mat[a];
                 }

diff --git a/src/MathFunctions.h b/src/MathFunctions.h
@@ -1753,7 +1753,13 @@ class CMathFunctions
             integ_spline.setValue(0, 0);
             for(uint i = 1; i < N; i++)
             {
-                res += (x[i] - x[i - 1]) * y[i - 1] + 0.5 * (x[i] - x[i - 1]) * (y[i] - y[i - 1]);
+                if(y[i - 1] > 0)
+                {
+                    if(y[i] == 0)
+                        res += 0.5 * (y[i + 1] - y[i - 1]) / (x[i + 1] - x[i - 1]) * (x[i] - x[i - 1]) * (x[i] - x[i - 1]);
+                    else
+                        res += (x[i] - x[i - 1]) * y[i - 1] + 0.5 * (x[i] - x[i - 1]) * (y[i] - y[i - 1]);
+                }
                 integ_spline.setValue(i, res);
             }
         }