Merge pull request #13 from MartinAlvarezSergio/Pull-request

Suggested bug fixes to synchrotron

src/RadiativeTransfer.cpp

@@ -2186,9 +2186,8 @@ void CRadiativeTransfer::getSyncIntensity(photon_package * pp,
             // Set current index in photon package
             pp->setSpectralID(i_wave + i_extra * nr_used_wavelengths);
 
-            // Get wavelength/frequency of the photon package
-            double mult = 1e+26 * subpixel_fraction * tracer[i_det]->getDistanceFactor() * con_c /
-                          (pp->getFrequency() * pp->getFrequency());
+            // Convert W/m2/Hz/sr to Jy
+            double mult = 1e+26 * subpixel_fraction * tracer[i_det]->getDistanceFactor();
 
             // Include foreground extinction if necessary
             mult *= dust->getForegroundExtinction(tracer[i_det]->getWavelength(pp->getWavelength()));

src/Synchrotron.cpp

@@ -104,7 +104,7 @@ syn_param CSynchrotron::get_Thermal_Parameter(double n_e, double T_e, double l,
 
     n_e *= 1e-6; // n_e in 1/ccm
 
-    B *= 1e4; // B in mu Gauss
+    B *= 1e4; // B in Gauss
 
     double sin_theta = sin(theta);
     double cos_theta = cos(theta);
@@ -239,9 +239,10 @@ syn_param CSynchrotron::get_Power_Law_Parameter(double n_e,
     res.j_V = res.j_I / sqrt(nu / (3. * nu_c * sin_theta)) * getI_V_p(p, tan_theta);
 
     // additional correction for g_min>1 (see Reissl et al. 2018)
-    res.j_I /= (g_min * g_min);
-    res.j_Q /= (g_min * g_min);
-    res.j_V /= (g_min * g_min);
+    double gmin_den = pow(g_min, 1 - p);
+    res.j_I *= gmin_den; 
+    res.j_Q *= gmin_den; 
+    res.j_V *= gmin_den; 
 
     double du = sqrt(res.j_I * res.j_I - res.j_Q * res.j_Q);
 
@@ -272,9 +273,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 /= (g_min * g_min);
-    res.alpha_Q /= (g_min * g_min);
-    res.alpha_V /= (g_min * g_min);
+    res.alpha_I *= gmin_den; 
+    res.alpha_Q *= gmin_den; 
+    res.alpha_V *= gmin_den; 
 
     // converting back into SI;
     res.scale();

src/Synchrotron.h

@@ -11,10 +11,12 @@
 #define syn_e 4.80320680e-10 // electron charge    [Statcoulomb]
 #define syn_h 6.6260693e-27  // Planck constant    [erg/s]
 #define syn_kB 1.380662e-16  // Boltzmann constant [erg/K]
-#define syn_c 2.99792458e10  // speed of light     [m/s]
+#define syn_c 2.99792458e10  // speed of light     [cm/s]
 
-// conversion factor 1/cm -> 1/m
-#define syn_SI 100.0
+// conversion factor 100.0 (1/cm -> 1/m) * 1e-3 (erg s^-1 cm^-2 Hz^-1 -> W m^-2 Hz^-1) 
+#define syn_SI 0.1
+// conversion factor 100.0 (1/cm -> 1/m) 
+#define syn_SI_abs 100.0
 
 // container class for the parameters of sync. RT
 class syn_param
@@ -74,9 +76,9 @@ class syn_param
         j_Q *= syn_SI;
         j_V *= syn_SI;
 
-        alpha_I *= syn_SI;
-        alpha_Q *= syn_SI;
-        alpha_V *= syn_SI;
+        alpha_I *= syn_SI_abs;
+        alpha_Q *= syn_SI_abs;
+        alpha_V *= syn_SI_abs;
 
         kappa_Q *= syn_SI;
         kappa_V *= syn_SI;
@@ -490,7 +492,7 @@ class CSynchrotron
 
     double Gamma_A_p(double g_min, double g_max, double p)
     {
-        return Gamma((3. * p + 2.) / 12.) * Gamma((3. * p + 22.) / 12.) /
+        return Gamma((3. * p + 12.) / 12.) * Gamma((3. * p + 22.) / 12.) /
                (4. * (pow(g_min, 1. - p) - pow(g_max, 1. - p)));
     }