text: #include "openmc/bremsstrahlung.h"
  
  #include "openmc/constants.h"
  #include "openmc/material.h"
  #include "openmc/random_lcg.h"
  #include "openmc/search.h"
  #include "openmc/settings.h"
  
  #include "xtensor/xmath.hpp"
  
  namespace openmc {
  
  //==============================================================================
  // Global variables
  //==============================================================================
  
  namespace data {
  
  xt::xtensor<double, 1> ttb_e_grid;
  xt::xtensor<double, 1> ttb_k_grid;
  vector<Bremsstrahlung> ttb;
  
  } // namespace data
  
  //==============================================================================
  // Non-member functions
  //==============================================================================
  
  void thick_target_bremsstrahlung(Particle& p, double* E_lost)
  {
    if (p.material() == MATERIAL_VOID)
      return;
  
    int photon = static_cast<int>(ParticleType::photon);
    if (p.E() < settings::energy_cutoff[photon])
      return;
  
    // Get bremsstrahlung data for this material and particle type
    BremsstrahlungData* mat;
    if (p.type() == ParticleType::positron) {
      mat = &model::materials[p.material()]->ttb_->positron;
    } else {
      mat = &model::materials[p.material()]->ttb_->electron;
    }
  
    double e = std::log(p.E());
    auto n_e = data::ttb_e_grid.size();
  
    // Find the lower bounding index of the incident electron energy
    size_t j =
      lower_bound_index(data::ttb_e_grid.cbegin(), data::ttb_e_grid.cend(), e);
    if (j == n_e - 1)
      --j;
  
    // Get the interpolation bounds
    double e_l = data::ttb_e_grid(j);
    double e_r = data::ttb_e_grid(j + 1);
    double y_l = mat->yield(j);
    double y_r = mat->yield(j + 1);
  
    // Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF
    // interpolated in log energy, which can be interpreted as the probability
    // of index j+1
    double f = (e - e_l) / (e_r - e_l);
  
    // Get the photon number yield for the given energy using linear
    // interpolation on a log-log scale
    double y = std::exp(y_l + (y_r - y_l) * f);
  
    // Sample number of secondary bremsstrahlung photons
    int n = y + prn(p.current_seed());
  
    *E_lost = 0.0;
    if (n == 0)
      return;
  
    // Sample index of the tabulated PDF in the energy grid, j or j+1
    double c_max;
    int i_e;
    if (prn(p.current_seed()) <= f || j == 0) {
      i_e = j + 1;
  
      // Interpolate the maximum value of the CDF at the incoming particle
      // energy on a log-log scale
      double p_l = mat->pdf(i_e, i_e - 1);
      double p_r = mat->pdf(i_e, i_e);
      double c_l = mat->cdf(i_e, i_e - 1);
      double a = std::log(p_r / p_l) / (e_r - e_l) + 1.0;
      c_max = c_l + std::exp(e_l) * p_l / a * (std::exp(a * (e - e_l)) - 1.0);
    } else {
      i_e = j;
  
      // Maximum value of the CDF
      c_max = mat->cdf(i_e, i_e);
    }
  
    // Sample the energies of the emitted photons
    for (int i = 0; i < n; ++i) {
      // Generate a random number r and determine the index i for which
      // cdf(i) <= r*cdf,max <= cdf(i+1)
      double c = prn(p.current_seed()) * c_max;
      int i_w = lower_bound_index(&mat->cdf(i_e, 0), &mat->cdf(i_e, 0) + i_e, c);
  
      // Sample the photon energy
      double w_l = data::ttb_e_grid(i_w);
      double w_r = data::ttb_e_grid(i_w + 1);
      double p_l = mat->pdf(i_e, i_w);
      double p_r = mat->pdf(i_e, i_w + 1);
      double c_l = mat->cdf(i_e, i_w);
      double a = std::log(p_r / p_l) / (w_r - w_l) + 1.0;
      double w = std::exp(w_l) *
                 std::pow(a * (c - c_l) / (std::exp(w_l) * p_l) + 1.0, 1.0 / a);
  
      if (w > settings::energy_cutoff[photon]) {
        // Create secondary photon
        p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon);
        *E_lost += w;
      }
    }
  }
  
  } // namespace openmc
  