Add kappa_gamma estimators

gammapkt.cc

@@ -647,6 +647,16 @@ constexpr auto meanf_sigma(const double x) -> double {
   return tot;
 }
 
+// get the gamma-ray opacity (with the expected energy loss per interaction factor included)
+auto get_kappa(const Packet &pkt) -> double {
+  const double xx = H * pkt.nu_cmf / ME / CLIGHT / CLIGHT;
+  const int mgi = grid::get_cell_modelgridindex(pkt.where);
+  const double chi = ((meanf_sigma(xx) * grid::get_nnetot(mgi)) + get_chi_photo_electric_rf(pkt) +
+                      (sigma_pair_prod_rf(pkt) * (1. - (2.46636e+20 / pkt.nu_cmf))));
+  const double kappa = 1 / grid::get_rho(mgi) * chi;
+  return kappa;
+}
+
 // Subroutine to record the heating rate in a cell due to gamma rays.
 // By heating rate I mean, for now, really the rate at which the code is making
 // k-packets in that cell which will then convert into r-packets. This is (going
@@ -677,6 +687,9 @@ void update_gamma_dep(const Packet &pkt, const double dist, const int mgi, const
   assert_testmodeonly(heating_cont >= 0.);
   assert_testmodeonly(std::isfinite(heating_cont));
   atomicadd(globals::dep_estimator_gamma[nonemptymgi], heating_cont);
+
+  const double kappa_E_cont = heating_cont * get_kappa(pkt);
+  atomicadd(globals::estimator_gamma_kappa_integrated, kappa_E_cont);
 }
 
 // handle physical pair production event
@@ -1079,6 +1092,9 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) {
   // printout("initialise pol state of packet %g, %g, %g, %g,
   // %g\n",pkt.stokes_qu[0],pkt.stokes_qu[1],pkt.pol_dir[0],pkt.pol_dir[1],pkt.pol_dir[2]);
   // printout("pkt direction %g, %g, %g\n",pkt.dir[0],pkt.dir[1],pkt.dir[2]);
+
+  atomicadd(globals::estimator_gamma_kappa_decayspec, get_kappa(pkt) * pkt.e_cmf);
+  atomicadd(globals::estimator_gamma_nu_cmf_decayspec, pkt.nu_cmf * pkt.e_cmf);
 }
 
 __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {
@@ -1096,6 +1112,7 @@ __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {
 
   if (pkt.type != TYPE_GAMMA && pkt.type != TYPE_ESCAPE) {
     atomicadd(globals::timesteps[nts].gamma_dep_discrete, pkt.e_cmf);
+    atomicadd(globals::estimator_gamma_kappa_absorbedspec, get_kappa(pkt) * pkt.e_cmf);
 
     if constexpr (GAMMA_THERMALISATION_SCHEME != ThermalisationScheme::DETAILED) {
       // no transport, so the path-based gamma deposition estimator won't get updated unless we do it here

globals.h

@@ -212,6 +212,11 @@ inline std::vector<double> dep_estimator_positron;
 inline std::vector<double> dep_estimator_electron;
 inline std::vector<double> dep_estimator_alpha;
 
+inline double estimator_gamma_nu_cmf_decayspec;
+inline double estimator_gamma_kappa_integrated;
+inline double estimator_gamma_kappa_decayspec;
+inline double estimator_gamma_kappa_absorbedspec;
+
 inline int bfestimcount{0};
 
 // for USE_LUT_PHOTOION = true

sn3d.cc

@@ -617,6 +617,11 @@ void zero_estimators() {
   std::ranges::fill(globals::dep_estimator_electron, 0.);
   std::ranges::fill(globals::dep_estimator_alpha, 0.);
 
+  globals::estimator_gamma_nu_cmf_decayspec = 0.;
+  globals::estimator_gamma_kappa_integrated = 0.;
+  globals::estimator_gamma_kappa_decayspec = 0.;
+  globals::estimator_gamma_kappa_absorbedspec = 0.;
+
   if constexpr (USE_LUT_PHOTOION) {
     if (globals::nbfcontinua_ground > 0) {
       std::ranges::fill(globals::gammaestimator, 0.);
@@ -643,6 +648,21 @@ void normalise_deposition_estimators(int nts) {
   globals::timesteps[nts].electron_dep = 0.;
   globals::timesteps[nts].alpha_dep = 0.;
 
+  double gamma_dep_erg = 0.;
+  for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) {
+    gamma_dep_erg += (globals::dep_estimator_gamma[nonemptymgi] / nprocs);
+  }
+  globals::estimator_gamma_kappa_integrated /= gamma_dep_erg;
+  globals::estimator_gamma_kappa_decayspec /= globals::timesteps[nts].gamma_emission;
+  globals::estimator_gamma_kappa_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete;
+  globals::estimator_gamma_nu_cmf_decayspec /= globals::timesteps[nts].gamma_emission;
+  const double mean_en_mev = H * globals::estimator_gamma_nu_cmf_decayspec / MEV;
+
+  printout(
+      "timestep %d kappa_eff [cm2/g]: path-integrated %g decay_spec %g deposit_spec %g decay gamma mean E[MeV] %g\n",
+      nts, globals::estimator_gamma_kappa_integrated, globals::estimator_gamma_kappa_decayspec,
+      globals::estimator_gamma_kappa_absorbedspec, mean_en_mev);
+
   for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) {
     const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);