Merge pull request #1 from MartinAlvarezSergio/Pull-requests

Pull requests
MartinAlvarezSergio · Jul 26, 2023 · 16b8f53 · 16b8f53
2 parents dc080fa + a675da1
commit 16b8f53
Show file tree

Hide file tree

Showing 7 changed files with 31 additions and 12 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -17,6 +17,7 @@ jobs:
       matrix:
         compiler: ["gcc", "clang++"]
         generator: ["make", "ninja"]
+      fail-fast: false
     name: "Ubuntu compile.sh :: ${{ matrix.compiler }} :: ${{ matrix.generator }}"
     runs-on: ubuntu-latest
     steps:
@@ -40,6 +41,7 @@ jobs:
             build_type: "Debug"
             exe_linker_flags: "-fsanitize=address"
             cxx_flags: "-fsanitize=address"
+      fail-fast: false
     name: "Ubuntu :: ${{ matrix.compiler }} :: ${{ matrix.build_type }} :: ${{ matrix.cxx_flags }} :: system ccfits cfitsio"
     runs-on: ubuntu-latest
     steps:

diff --git a/manual.pdf b/manual.pdf
diff --git a/src/Dust.cpp b/src/Dust.cpp
@@ -2754,11 +2754,12 @@ void CDustComponent::preCalcMieScatteringProb()
 #pragma omp parallel for
     for(int a = 0; a < int(nr_of_dust_species); a++)
     {
+        // Init arrays of interp
+        avg_scattering_frac[a] = new interp[nr_of_wavelength];
+        phase_pdf[a] = new interp[nr_of_wavelength];
+
         if(sizeIndexUsed(a))
         {
-            // Init arrays of interp
-            avg_scattering_frac[a] = new interp[nr_of_wavelength];
-            phase_pdf[a] = new interp[nr_of_wavelength];
             for(uint w = 0; w < nr_of_wavelength; w++)
             {
                 // Init pointer arrays
@@ -4748,7 +4749,11 @@ void CDustComponent::getEscapePhotonMie(CGridBasic * grid,
     }
     else
     {
-        cout << "\nHINT: Photon package intensity or first scattering matrix element is zero!\n" << endl;
+        if(tmp_stokes.I() <= 0.0)
+            cout << "\nERROR: Photon package intensity is zero or negative!\n" << endl;
+        if(mat_sca(0, 0) <= 0.0)
+            cout << "\nERROR: First scattering matrix element is zero or negative!\n" << endl;
+
         tmp_stokes.clear();
         pp_escape->setStokesVector(tmp_stokes);
         return;
@@ -4926,7 +4931,11 @@ void CDustComponent::miesca(photon_package * pp, uint a, CRandomGenerator * rand
     }
     else
     {
-        cout << "\nHINT: Photon package intensity or first scattering matrix element is zero!\n" << endl;
+        if(tmp_stokes.I() <= 0.0)
+            cout << "\nERROR: Photon package intensity is zero or negative!\n" << endl;
+        if(mat_sca(0, 0) <= 0.0)
+            cout << "\nERROR: First scattering matrix element is zero or negative!\n" << endl;
+
         tmp_stokes.clear();
         pp->setStokesVector(tmp_stokes);
         return;
@@ -4958,7 +4967,7 @@ void CDustComponent::miesca(photon_package * pp, uint a, CRandomGenerator * rand
         run_counter++;
     }
     if(run_counter == 1000 || abs(root_phi) > 1e-10)
-        cout << "\nERROR: No phi found\n" << endl;
+        cout << "\nERROR: No scattering angle phi found!\n" << endl;
 
     phi = PI - gamma + phi;
 

diff --git a/src/Dust.h b/src/Dust.h
@@ -124,7 +124,7 @@ class CDustComponent
         nr_of_mixtures = 0;
         calorimetry_type = 0;
         alignment = ALIG_PA;
-        phID = 1;
+        phID = 0;
         nr_of_dust_species = 0;
         nr_of_incident_angles = 0;
         nr_of_scat_theta = 0;

diff --git a/src/Parameters.h b/src/Parameters.h
@@ -682,10 +682,11 @@ class parameters
 
     uint getPhaseFunctionID(uint i) const
     {
+        if(phIDs.size() == 1)
+            return phIDs[0];
         if(i < phIDs.size())
             return phIDs[i];
-        else
-            return PH_ISO;
+        return PH_ISO;
     }
 
     double getFHighJ() const

diff --git a/src/Synchrotron.cpp b/src/Synchrotron.cpp
@@ -240,9 +240,15 @@ syn_param CSynchrotron::get_Power_Law_Parameter(double n_e,
 
     // additional correction for g_min>1 (see Reissl et al. 2018)
     double gmin_den = pow(g_min, 1 - p);
+<<<<<<< HEAD
     res.j_I *= gmin_den; // res.j_I /= (g_min * g_min); // SMA: old was hardcoded for p = 3.0
     res.j_Q *= gmin_den; // res.j_Q /= (g_min * g_min); // SMA: old was hardcoded for p = 3.0
     res.j_V *= gmin_den; // res.j_V /= (g_min * g_min); // SMA: old was hardcoded for p = 3.0
+=======
+    res.j_I *= gmin_den; 
+    res.j_Q *= gmin_den; 
+    res.j_V *= gmin_den; 
+>>>>>>> tmp
 
     double du = sqrt(res.j_I * res.j_I - res.j_Q * res.j_Q);
 
@@ -273,9 +279,9 @@ syn_param CSynchrotron::get_Power_Law_Parameter(double n_e,
                        // * log(g_min)/tan_theta;
 
     // additional correction for g_min>1 (see Reissl et al. 2018)
-    res.alpha_I *= gmin_den; //res.alpha_I /= (g_min * g_min); // SMA: old was hardcoded for p = 3.0
-    res.alpha_Q *= gmin_den; //res.alpha_Q /= (g_min * g_min); // SMA: old was hardcoded for p = 3.0
-    res.alpha_V *= gmin_den; //res.alpha_V /= (g_min * g_min); // SMA: old was hardcoded for p = 3.0
+    res.alpha_I *= gmin_den; 
+    res.alpha_Q *= gmin_den; 
+    res.alpha_V *= gmin_den; 
 
     // converting back into SI;
     res.scale();

diff --git a/src/Synchrotron.h b/src/Synchrotron.h
@@ -495,6 +495,7 @@ class CSynchrotron
         return Gamma((3. * p + 12.) / 12.) * Gamma((3. * p + 22.) / 12.) /
                (4. * (pow(g_min, 1. - p) - pow(g_max, 1. - p)));
     }
+
 
     double getI_Q_p(double p)
     {