Update gammapkt.cc

gammapkt.cc

@@ -27,6 +27,8 @@
 namespace gammapkt {
 // Code for handing gamma rays - creation and propagation
 
+constexpr bool USE_XCOM_GAMMAPHOTOION = false;
+
 struct gamma_line {
   double energy{};  // in erg
   double probability{};
@@ -573,52 +575,70 @@ static auto get_chi_photo_electric_rf(const struct packet *pkt_ptr) -> double {
 
   if (globals::gamma_kappagrey < 0) {
     chi_cmf = 0.;
-    // Cross sections from Equation 2 of Ambwani & Sutherland (1988), attributed to Veigele (1973)
-
-    // 2.41326e19 Hz = 100 keV / H
-    // const double hnu_over_100kev = pkt_ptr->nu_cmf / 2.41326e+19;
-    // for interpolation compute frequency in MeV
-    const double hnu_over_1MeV = pkt_ptr->nu_cmf / 2.41326e+20;
-    const double log10_hnu_over_1MeV = log10(hnu_over_1MeV);
-
-    // double sigma_cmf_cno = 0.0448e-24 * pow(hnu_over_100kev, -3.2);
-    const int numb_elements = get_nelements();
-    for (int i = 0; i < numb_elements; i++) {
-      // determine charge number:
-      int Z = get_atomicnumber(i);
-      double n_i = grid::get_elem_numberdens(mgi, i);  // number density in the current cell
-      if (n_i == 0) {
-        continue;
-      }
-      // get indices of lower and upper boundary
-      int E_gtr_idx = -2;  //
-      auto numb_energies = std::ssize(photoion_data[Z - 1]);
-      for (int j = 0; j < numb_energies; j++) {
-        if (photoion_data[Z - 1][j].energy > hnu_over_1MeV) {
-          E_gtr_idx = j;
-          break;
+    if (!USE_XCOM_GAMMAPHOTOION) {
+      // Cross sections from Equation 2 of Ambwani & Sutherland (1988), attributed to Veigele (1973)
+
+      // 2.41326e19 Hz = 100 keV / H
+      const double hnu_over_100kev = pkt_ptr->nu_cmf / 2.41326e+19;
+
+      // double sigma_cmf_cno = 0.0448e-24 * pow(hnu_over_100kev, -3.2);
+
+      const double sigma_cmf_si = 1.16e-24 * pow(hnu_over_100kev, -3.13);
+
+      const double sigma_cmf_fe = 25.7e-24 * pow(hnu_over_100kev, -3.0);
+
+      // Now need to multiply by the particle number density.
+
+      const double chi_cmf_si = sigma_cmf_si * (rho / MH / 28);
+      // Assumes Z = 14. So mass = 28.
+
+      const double chi_cmf_fe = sigma_cmf_fe * (rho / MH / 56);
+      // Assumes Z = 28. So mass = 56.
+
+      const double f_fe = grid::get_ffegrp(mgi);
+
+      chi_cmf = (chi_cmf_fe * f_fe) + (chi_cmf_si * (1. - f_fe));
+    } else {
+      const double hnu_over_1MeV = pkt_ptr->nu_cmf / 2.41326e+20;
+      const double log10_hnu_over_1MeV = log10(hnu_over_1MeV);
+      const int numb_elements = get_nelements();
+      for (int i = 0; i < numb_elements; i++) {
+        // determine charge number:
+        int Z = get_atomicnumber(i);
+        double n_i = grid::get_elem_numberdens(mgi, i);  // number density in the current cell
+        if (n_i == 0) {
+          continue;
         }
+        // get indices of lower and upper boundary
+        int E_gtr_idx = -2;  //
+        auto numb_energies = std::ssize(photoion_data[Z - 1]);
+        for (int j = 0; j < numb_energies; j++) {
+          if (photoion_data[Z - 1][j].energy > hnu_over_1MeV) {
+            E_gtr_idx = j;
+            break;
+          }
+        }
+        if (E_gtr_idx == -1) {  // packet energy smaller than all tabulated values
+          chi_cmf += photoion_data[i][0].sigma_xcom * n_i;
+          continue;
+        }
+        if (E_gtr_idx == -2) {  // packet energy greater than all tabulated values
+          chi_cmf += photoion_data[i][numb_energies - 1].sigma_xcom * n_i;
+          continue;
+        }
+        int E_smaller_idx = E_gtr_idx - 1;
+        double log10_E = log10_hnu_over_1MeV;
+        double log10_E_gtr = log10(photoion_data[Z - 1][E_gtr_idx].energy);
+        double log10_E_smaller = log10(photoion_data[Z - 1][E_smaller_idx].energy);
+        double log10_sigma_lower = log10(photoion_data[Z - 1][E_smaller_idx].sigma_xcom);
+        double log10_sigma_gtr = log10(photoion_data[Z - 1][E_gtr_idx].sigma_xcom);
+        // interpolate or extrapolate, both linear in log10-log10 space
+        double log10_intpol = log10_E_smaller + (log10_sigma_gtr - log10_sigma_lower) /
+                                                    (log10_E_gtr - log10_E_smaller) * (log10_E - log10_E_smaller);
+        double sigma_intpol = pow(10., log10_intpol) * 1.0e-24;  // now in cm^2
+        double chi_cmf_contrib = sigma_intpol * n_i;
+        chi_cmf += chi_cmf_contrib;
       }
-      if (E_gtr_idx == -1) {  // packet energy smaller than all tabulated values
-        chi_cmf += photoion_data[i][0].sigma_xcom * n_i;
-        continue;
-      }
-      if (E_gtr_idx == -2) {  // packet energy greater than all tabulated values
-        chi_cmf += photoion_data[i][numb_energies - 1].sigma_xcom * n_i;
-        continue;
-      }
-      int E_smaller_idx = E_gtr_idx - 1;
-      double log10_E = log10_hnu_over_1MeV;
-      double log10_E_gtr = log10(photoion_data[Z - 1][E_gtr_idx].energy);
-      double log10_E_smaller = log10(photoion_data[Z - 1][E_smaller_idx].energy);
-      double log10_sigma_lower = log10(photoion_data[Z - 1][E_smaller_idx].sigma_xcom);
-      double log10_sigma_gtr = log10(photoion_data[Z - 1][E_gtr_idx].sigma_xcom);
-      // interpolate or extrapolate, both linear in log10-log10 space
-      double log10_intpol = log10_E_smaller + (log10_sigma_gtr - log10_sigma_lower) / (log10_E_gtr - log10_E_smaller) *
-                                                  (log10_E - log10_E_smaller);
-      double sigma_intpol = pow(10., log10_intpol) * 1.0e-24;  // now in cm^2
-      double chi_cmf_contrib = sigma_intpol * n_i;
-      chi_cmf += chi_cmf_contrib;
     }
   } else {
     chi_cmf = globals::gamma_kappagrey * rho;