From 844c6066e2dc70eb3b7a0e54752f3abfc732d624 Mon Sep 17 00:00:00 2001 From: mlietzow Date: Thu, 21 Sep 2023 16:04:27 +0200 Subject: [PATCH] fixed bug if components with different grain sizes were used including gaps in the distribution --- src/Dust.cpp | 18 +++++++++--------- src/Dust.h | 19 +++++++++++-------- src/MathFunctions.h | 8 +++++++- 3 files changed, 27 insertions(+), 18 deletions(-) diff --git a/src/Dust.cpp b/src/Dust.cpp index 928b67e..dd71785 100644 --- 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 index adebed9..732a062 100644 --- 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 index 291773b..9284d21 100644 --- 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); } }