#include "openmc/particle.h"

#include <algorithm> // copy, min
#include <cmath>     // log, abs

#include <fmt/core.h>

#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/collision_track.h"
#include "openmc/constants.h"
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/particle_data.h"
#include "openmc/photon.h"
#include "openmc/physics.h"
#include "openmc/physics_mg.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/surface.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/track_output.h"
#include "openmc/weight_windows.h"

#ifdef OPENMC_DAGMC_ENABLED
#include "DagMC.hpp"
#endif

namespace openmc {

//==============================================================================
// Particle implementation
//==============================================================================

double Particle::speed(double E) const
{
  // Determine mass in eV/c^2
  double mass = this->mass();
  // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
  return C_LIGHT * std::sqrt(E * (E + 2 * mass)) / (E + mass);
}

double Particle::speed() const
{
  if (settings::run_CE) {
    return speed(this->E());
  } else {
    auto mat = this->material();
    if (mat == MATERIAL_VOID)
      return 1.0 / data::mg.default_inverse_velocity_[this->g()];
    auto& macro_xs = data::mg.macro_xs_[mat];
    int macro_t = this->mg_xs_cache().t;
    int macro_a = macro_xs.get_angle_index(this->u());
    return 1.0 / macro_xs.get_xs(
                   MgxsType::INVERSE_VELOCITY, this->g(), macro_t, macro_a);
  }
}

double Particle::mass() const
{
  switch (type().pdg_number()) {
  case PDG_NEUTRON:
    return MASS_NEUTRON_EV;
  case PDG_ELECTRON:
  case PDG_POSITRON:
    return MASS_ELECTRON_EV;
  case PDG_PHOTON:
    return 0.0;
  default:
    return this->type().mass() * AMU_EV;
  }
}

bool Particle::create_secondary(
  double wgt, Direction u, double E, ParticleType type)
{
  // If energy is below cutoff for this particle, don't create secondary
  // particle
  int idx = type.transport_index();
  if (idx == C_NONE) {
    return false;
  }
  if (E < settings::energy_cutoff[idx]) {
    return false;
  }

  // Increment number of secondaries created (for ParticleProductionFilter)
  n_secondaries()++;

  SourceSite bank;
  bank.particle = type;
  bank.wgt = wgt;
  bank.r = r();
  bank.u = u;
  bank.E = settings::run_CE ? E : g();
  bank.time = time();
  bank_second_E() += bank.E;
  bank.parent_id = current_work();
  if (settings::use_shared_secondary_bank) {
    bank.progeny_id = n_progeny()++;
  }
  bank.wgt_born = wgt_born();
  bank.wgt_ww_born = wgt_ww_born();
  bank.n_split = n_split();

  local_secondary_bank().emplace_back(bank);
  return true;
}

void Particle::split(double wgt)
{
  SourceSite bank;
  bank.particle = type();
  bank.wgt = wgt;
  bank.r = r();
  bank.u = u();
  bank.E = settings::run_CE ? E() : g();
  bank.time = time();

  // Convert signed index to a signed surface ID
  if (surface() == SURFACE_NONE) {
    bank.surf_id = SURFACE_NONE;
  } else {
    int surf_id = model::surfaces[surface_index()]->id_;
    bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
  }

  bank.wgt_born = wgt_born();
  bank.wgt_ww_born = wgt_ww_born();
  bank.n_split = n_split();
  bank.n_collision = n_collision();
  bank.parent_id = current_work();
  if (settings::use_shared_secondary_bank) {
    bank.progeny_id = n_progeny()++;
  }

  local_secondary_bank().emplace_back(bank);
}

void Particle::from_source(const SourceSite* src)
{
  // Reset some attributes
  clear();
  surface() = SURFACE_NONE;
  cell_born() = C_NONE;
  material() = C_NONE;
  n_collision() = src->n_collision;
  fission() = false;
  zero_flux_derivs();
  lifetime() = 0.0;
#ifdef OPENMC_DAGMC_ENABLED
  history().reset();
#endif

  // Copy attributes from source bank site
  type() = src->particle;
  wgt() = src->wgt;
  wgt_last() = src->wgt;
  r() = src->r;
  u() = src->u;
  r_born() = src->r;
  r_last_current() = src->r;
  r_last() = src->r;
  u_last() = src->u;
  if (settings::run_CE) {
    E() = src->E;
    g() = 0;
  } else {
    g() = static_cast<int>(src->E);
    g_last() = static_cast<int>(src->E);
    E() = data::mg.energy_bin_avg_[g()];
  }
  E_last() = E();
  time() = src->time;
  time_last() = src->time;
  parent_nuclide() = src->parent_nuclide;
  delayed_group() = src->delayed_group;

  // Convert signed surface ID to signed index
  if (src->surf_id != SURFACE_NONE) {
    auto it = model::surface_map.find(std::abs(src->surf_id));
    if (it != model::surface_map.end()) {
      int index_plus_one = it->second + 1;
      surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
    }
  }

  wgt_born() = src->wgt_born;
  wgt_ww_born() = src->wgt_ww_born;
  n_split() = src->n_split;
}

void Particle::event_calculate_xs()
{
  // Set the random number stream
  stream() = STREAM_TRACKING;

  // Store pre-collision particle properties
  wgt_last() = wgt();
  E_last() = E();
  u_last() = u();
  r_last() = r();
  time_last() = time();

  // Reset event variables
  event() = TallyEvent::KILL;
  event_nuclide() = NUCLIDE_NONE;
  event_mt() = REACTION_NONE;

  // If the cell hasn't been determined based on the particle's location,
  // initiate a search for the current cell. This generally happens at the
  // beginning of the history and again for any secondary particles
  if (lowest_coord().cell() == C_NONE) {
    if (!exhaustive_find_cell(*this)) {
      mark_as_lost(
        "Could not find the cell containing particle " + std::to_string(id()));
      return;
    }

    // Set birth cell attribute
    if (cell_born() == C_NONE)
      cell_born() = lowest_coord().cell();

    // Initialize last cells from current cell
    for (int j = 0; j < n_coord(); ++j) {
      cell_last(j) = coord(j).cell();
    }
    n_coord_last() = n_coord();
  }

  // Write particle track.
  if (write_track())
    write_particle_track(*this);

  if (settings::check_overlaps)
    check_cell_overlap(*this);

  // Calculate microscopic and macroscopic cross sections
  if (material() != MATERIAL_VOID) {
    if (settings::run_CE) {
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
          density_mult() != density_mult_last()) {
        // If the material is the same as the last material and the
        // temperature hasn't changed, we don't need to lookup cross
        // sections again.
        model::materials[material()]->calculate_xs(*this);
      }
    } else {
      // Get the MG data; unlike the CE case above, we have to re-calculate
      // cross sections for every collision since the cross sections may
      // be angle-dependent
      data::mg.macro_xs_[material()].calculate_xs(*this);

      // Update the particle's group while we know we are multi-group
      g_last() = g();
    }
  } else {
    macro_xs().total = 0.0;
    macro_xs().absorption = 0.0;
    macro_xs().fission = 0.0;
    macro_xs().nu_fission = 0.0;
  }
}

void Particle::event_advance()
{
  // Find the distance to the nearest boundary
  boundary() = distance_to_boundary(*this);

  // Sample a distance to collision
  if (type() == ParticleType::electron() ||
      type() == ParticleType::positron()) {
    collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
  } else if (macro_xs().total == 0.0) {
    collision_distance() = INFINITY;
  } else {
    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
  }

  double speed = this->speed();
  double time_cutoff = settings::time_cutoff[type().transport_index()];
  double distance_cutoff =
    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;

  // Select smaller of the three distances
  double distance =
    std::min({boundary().distance(), collision_distance(), distance_cutoff});

  // Advance particle in space and time
  this->move_distance(distance);
  double dt = distance / speed;
  this->time() += dt;
  this->lifetime() += dt;

  // Score timed track-length tallies
  if (!model::active_timed_tracklength_tallies.empty()) {
    score_timed_tracklength_tally(*this, distance);
  }

  // Score track-length tallies
  if (!model::active_tracklength_tallies.empty()) {
    score_tracklength_tally(*this, distance);
  }

  // Score track-length estimate of k-eff
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
  }

  // Score flux derivative accumulators for differential tallies.
  if (!model::active_tallies.empty()) {
    score_track_derivative(*this, distance);
  }

  // Set particle weight to zero if it hit the time boundary
  if (distance == distance_cutoff) {
    wgt() = 0.0;
  }

  // Clear surface component if distance is long enough
  if (distance > TINY_BIT)
    surface() = SURFACE_NONE;
}

void Particle::event_cross_surface()
{
  // Saving previous cell data
  for (int j = 0; j < n_coord(); ++j) {
    cell_last(j) = coord(j).cell();
  }
  n_coord_last() = n_coord();

  // Set surface that particle is on and adjust coordinate levels
  surface() = boundary().surface();
  n_coord() = boundary().coord_level();

  if (boundary().lattice_translation()[0] != 0 ||
      boundary().lattice_translation()[1] != 0 ||
      boundary().lattice_translation()[2] != 0) {
    // Particle crosses lattice boundary

    int i_lattice = coord(boundary().coord_level() - 1).lattice();
    bool verbose = settings::verbosity >= 10 || trace();
    cross_lattice(*this, boundary(), verbose);
    event() = TallyEvent::LATTICE;

    // Score cell to cell partial currents
    if (!model::active_surface_tallies.empty()) {
      auto& lat {*model::lattices[i_lattice]};
      bool is_valid;
      Direction normal =
        lat.get_normal(boundary().lattice_translation(), is_valid);
      if (is_valid) {
        normal /= normal.norm();
        score_surface_tally(*this, model::active_surface_tallies, normal);
      }
    }

  } else {

    const auto& surf {*model::surfaces[surface_index()].get()};

    // Particle crosses surface
    // If BC, add particle to surface source before crossing surface
    if (surf.surf_source_ && surf.bc_) {
      add_surf_source_to_bank(*this, surf);
    }
    this->cross_surface(surf);
    // If no BC, add particle to surface source after crossing surface
    if (surf.surf_source_ && !surf.bc_) {
      add_surf_source_to_bank(*this, surf);
    }
    if (settings::weight_window_checkpoint_surface) {
      apply_weight_windows(*this);
    }
    event() = TallyEvent::SURFACE;

    // Score cell to cell partial currents
    if (!model::active_surface_tallies.empty()) {
      Direction normal = surf.normal(r());
      normal /= normal.norm();
      score_surface_tally(*this, model::active_surface_tallies, normal);
    }
  }
}

void Particle::event_collide()
{

  // Score collision estimate of keff
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
    keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
  }

  // Score surface current tallies -- this has to be done before the collision
  // since the direction of the particle will change and we need to use the
  // pre-collision direction to figure out what mesh surfaces were crossed

  if (!model::active_meshsurf_tallies.empty())
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);

  // Preserve whether the particle is still associated with a recently crossed
  // surface so that a direction change during a near-surface collision can be
  // reconciled afterward. The surface marker is no longer needed during the
  // collision itself.
  const bool near_surface = surface() != SURFACE_NONE;
  surface() = SURFACE_NONE;

  if (settings::run_CE) {
    collision(*this);
  } else {
    collision_mg(*this);
  }

  // Collision track feature to recording particle interaction
  if (settings::collision_track) {
    collision_track_record(*this);
  }

  // Score collision estimator tallies -- this is done after a collision
  // has occurred rather than before because we need information on the
  // outgoing energy for any tallies with an outgoing energy filter
  if (!model::active_collision_tallies.empty())
    score_collision_tally(*this);
  if (!model::active_analog_tallies.empty()) {
    if (settings::run_CE) {
      score_analog_tally_ce(*this);
    } else {
      score_analog_tally_mg(*this);
    }
  }

  if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {
    pht_collision_energy();
  }

  // Reset banked weight during collision
  n_bank() = 0;
  bank_second_E() = 0.0;
  wgt_bank() = 0.0;

  // Clear number of secondaries in this collision. This is
  // distinct from the number of created neutrons n_bank() above!
  n_secondaries() = 0;

  zero_delayed_bank();

  // Reset fission logical
  fission() = false;

  // Save coordinates for tallying purposes
  r_last_current() = r();

  // Set last material to none since cross sections will need to be
  // re-evaluated
  material_last() = C_NONE;

  // Set all directions to base level -- right now, after a collision, only
  // the base level directions are changed
  for (int j = 0; j < n_coord() - 1; ++j) {
    if (coord(j + 1).rotated()) {
      // If next level is rotated, apply rotation matrix
      const auto& m {model::cells[coord(j).cell()]->rotation_};
      const auto& u {coord(j).u()};
      coord(j + 1).u() = u.rotate(m);
    } else {
      // Otherwise, copy this level's direction
      coord(j + 1).u() = coord(j).u();
    }
  }

  // Score flux derivative accumulators for differential tallies.
  if (!model::active_tallies.empty())
    score_collision_derivative(*this);

#ifdef OPENMC_DAGMC_ENABLED
  history().reset();
#endif

  if (near_surface && alive())
    reconcile_cell_after_collision(*this);
}

void Particle::event_revive_from_secondary(const SourceSite& site)
{
  // Write final position for the previous track (skip if this is a freshly
  // constructed particle with no prior track, e.g., Phase 2 of shared
  // secondary transport)
  if (write_track() && n_event() > 0) {
    write_particle_track(*this);
  }

  from_source(&site);

  n_event() = 0;
  if (!settings::use_shared_secondary_bank) {
    n_tracks()++;
  }
  bank_second_E() = 0.0;

  // Subtract secondary particle energy from interim pulse-height results.
  // In shared secondary mode, this subtraction was already done on the parent
  // particle during create_secondary(), so skip it here.
  if (!settings::use_shared_secondary_bank &&
      !model::active_pulse_height_tallies.empty() && this->type().is_photon()) {
    // Since the birth cell of the particle has not been set we
    // have to determine it before the energy of the secondary particle can be
    // removed from the pulse-height of this cell.
    if (lowest_coord().cell() == C_NONE) {
      bool verbose = settings::verbosity >= 10 || trace();
      if (!exhaustive_find_cell(*this, verbose)) {
        mark_as_lost("Could not find the cell containing particle " +
                     std::to_string(id()));
        return;
      }
      // Set birth cell attribute
      if (cell_born() == C_NONE)
        cell_born() = lowest_coord().cell();

      // Initialize last cells from current cell
      for (int j = 0; j < n_coord(); ++j) {
        cell_last(j) = coord(j).cell();
      }
      n_coord_last() = n_coord();
    }
    pht_secondary_particles();
  }

  // Enter new particle in particle track file
  if (write_track())
    add_particle_track(*this);
}

void Particle::event_check_limit_and_revive()
{
  // If particle has too many events, display warning and kill it
  n_event()++;
  if (n_event() == settings::max_particle_events) {
    warning("Particle " + std::to_string(id()) +
            " underwent maximum number of events.");
    wgt() = 0.0;
  }

  // In non-shared-secondary mode, revive from local secondary bank
  if (!alive() && !settings::use_shared_secondary_bank &&
      !local_secondary_bank().empty()) {
    SourceSite& site = local_secondary_bank().back();
    event_revive_from_secondary(site);
    local_secondary_bank().pop_back();
  }
}

void Particle::event_death()
{
#ifdef OPENMC_DAGMC_ENABLED
  history().reset();
#endif

  // Finish particle track output.
  if (write_track()) {
    write_particle_track(*this);
    finalize_particle_track(*this);
  }

  // Contribute tally reduction variables to global accumulator
  const auto k_absorption = keff_tally_absorption();
  const auto k_collision = keff_tally_collision();
  const auto k_tracklength = keff_tally_tracklength();
  const auto leakage = keff_tally_leakage();

  if (settings::run_mode == RunMode::EIGENVALUE) {
    if (k_absorption != 0.0) {
#pragma omp atomic
      global_tally_absorption += k_absorption;
    }
    if (k_collision != 0.0) {
#pragma omp atomic
      global_tally_collision += k_collision;
    }
    if (k_tracklength != 0.0) {
#pragma omp atomic
      global_tally_tracklength += k_tracklength;
    }
  }
  if (leakage != 0.0) {
#pragma omp atomic
    global_tally_leakage += leakage;
  }

  // Reset particle tallies once accumulated
  keff_tally_absorption() = 0.0;
  keff_tally_collision() = 0.0;
  keff_tally_tracklength() = 0.0;
  keff_tally_leakage() = 0.0;

  if (!model::active_pulse_height_tallies.empty()) {
    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
  }

  // Accumulate track count for this particle history
  if (!settings::use_shared_secondary_bank) {
#pragma omp atomic
    simulation::simulation_tracks_completed += n_tracks();
  }

  // Record the number of progeny created by this particle.
  // This data will be used to efficiently sort the fission bank.
  if (settings::run_mode == RunMode::EIGENVALUE ||
      settings::use_shared_secondary_bank) {
    simulation::progeny_per_particle[current_work()] = n_progeny();
  }
}

void Particle::pht_collision_energy()
{
  // Adds the energy particles lose in a collision to the pulse-height

  // determine index of cell in pulse_height_cells
  auto it = std::find(model::pulse_height_cells.begin(),
    model::pulse_height_cells.end(), lowest_coord().cell());

  if (it != model::pulse_height_cells.end()) {
    int index = std::distance(model::pulse_height_cells.begin(), it);
    pht_storage()[index] += E_last() - E();

    // If the energy of the particle is below the cutoff, it will not be sampled
    // so its energy is added to the pulse-height in the cell
    int photon = ParticleType::photon().transport_index();
    if (E() < settings::energy_cutoff[photon]) {
      pht_storage()[index] += E();
    }
  }
}

void Particle::pht_secondary_particles()
{
  // Removes the energy of secondary produced particles from the pulse-height

  // determine index of cell in pulse_height_cells
  auto it = std::find(model::pulse_height_cells.begin(),
    model::pulse_height_cells.end(), cell_born());

  if (it != model::pulse_height_cells.end()) {
    int index = std::distance(model::pulse_height_cells.begin(), it);
    pht_storage()[index] -= E();
  }
}

void Particle::cross_surface(const Surface& surf)
{

  if (settings::verbosity >= 10 || trace()) {
    write_message(1, "    Crossing surface {}", surf.id_);
  }

  // Handle any applicable boundary conditions.
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
      settings::run_mode != RunMode::VOLUME) {
    surf.bc_->handle_particle(*this, surf);
    return;
  }

  // ==========================================================================
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL

#ifdef OPENMC_DAGMC_ENABLED
  // in DAGMC, we know what the next cell should be
  if (surf.geom_type() == GeometryType::DAG) {
    int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
                       lowest_coord().universe()) -
                     1;
    // save material, temperature, and density multiplier
    material_last() = material();
    sqrtkT_last() = sqrtkT();
    density_mult_last() = density_mult();
    // set new cell value
    lowest_coord().cell() = i_cell;
    auto& cell = model::cells[i_cell];

    cell_instance() = 0;
    if (cell->distribcell_index_ >= 0)
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);

    material() = cell->material(cell_instance());
    sqrtkT() = cell->sqrtkT(cell_instance());
    density_mult() = cell->density_mult(cell_instance());
    return;
  }
#endif

  bool verbose = settings::verbosity >= 10 || trace();
  if (neighbor_list_find_cell(*this, verbose)) {
    return;
  }

  // ==========================================================================
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS

  // Remove lower coordinate levels
  n_coord() = 1;
  bool found = exhaustive_find_cell(*this, verbose);

  if (settings::run_mode != RunMode::PLOTTING && (!found)) {
    // If a cell is still not found, there are two possible causes: 1) there is
    // a void in the model, and 2) the particle hit a surface at a tangent. If
    // the particle is really traveling tangent to a surface, if we move it
    // forward a tiny bit it should fix the problem.

    surface() = SURFACE_NONE;
    n_coord() = 1;
    r() += TINY_BIT * u();

    // Couldn't find next cell anywhere! This probably means there is an actual
    // undefined region in the geometry.

    if (!exhaustive_find_cell(*this, verbose)) {
      mark_as_lost("After particle " + std::to_string(id()) +
                   " crossed surface " + std::to_string(surf.id_) +
                   " it could not be located in any cell and it did not leak.");
      return;
    }
  }
}

void Particle::cross_vacuum_bc(const Surface& surf)
{
  // Score any surface current tallies -- note that the particle is moved
  // forward slightly so that if the mesh boundary is on the surface, it is
  // still processed

  if (!model::active_meshsurf_tallies.empty()) {
    // TODO: Find a better solution to score surface currents than
    // physically moving the particle forward slightly

    r() += TINY_BIT * u();
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
  }

  // Score to global leakage tally
  keff_tally_leakage() += wgt();

  // Kill the particle
  wgt() = 0.0;

  // Display message
  if (settings::verbosity >= 10 || trace()) {
    write_message(1, "    Leaked out of surface {}", surf.id_);
  }
}

void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
{
  // Do not handle reflective boundary conditions on lower universes
  if (n_coord() != 1) {
    mark_as_lost("Cannot reflect particle " + std::to_string(id()) +
                 " off surface in a lower universe.");
    return;
  }

  // Score surface currents since reflection causes the direction of the
  // particle to change. For surface filters, we need to score the tallies
  // twice, once before the particle's surface attribute has changed and
  // once after. For mesh surface filters, we need to artificially move
  // the particle slightly back in case the surface crossing is coincident
  // with a mesh boundary

  if (!model::active_surface_tallies.empty()) {
    Direction normal = surf.normal(r());
    normal /= normal.norm();
    score_surface_tally(*this, model::active_surface_tallies, normal);
  }

  if (!model::active_meshsurf_tallies.empty()) {
    Position r {this->r()};
    this->r() -= TINY_BIT * u();
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
    this->r() = r;
  }

  // Set the new particle direction
  u() = new_u;

  // Reassign particle's cell and surface
  coord(0).cell() = cell_last(0);
  surface() = -surface();

  // If a reflective surface is coincident with a lattice or universe
  // boundary, it is necessary to redetermine the particle's coordinates in
  // the lower universes.
  // (unless we're using a dagmc model, which has exactly one universe)
  n_coord() = 1;
  if (surf.geom_type() != GeometryType::DAG &&
      !neighbor_list_find_cell(*this)) {
    mark_as_lost("Couldn't find particle after reflecting from surface " +
                 std::to_string(surf.id_) + ".");
    return;
  }

  // Set previous coordinate going slightly past surface crossing
  r_last_current() = r() + TINY_BIT * u();

  // Diagnostic message
  if (settings::verbosity >= 10 || trace()) {
    write_message(1, "    Reflected from surface {}", surf.id_);
  }
}

void Particle::cross_periodic_bc(
  const Surface& surf, Position new_r, Direction new_u, int new_surface)
{
  // Do not handle periodic boundary conditions on lower universes
  if (n_coord() != 1) {
    mark_as_lost(
      "Cannot transfer particle " + std::to_string(id()) +
      " across surface in a lower universe. Boundary conditions must be "
      "applied to root universe.");
    return;
  }

  // Score surface currents since reflection causes the direction of the
  // particle to change -- artificially move the particle slightly back in
  // case the surface crossing is coincident with a mesh boundary
  if (!model::active_meshsurf_tallies.empty()) {
    Position r {this->r()};
    this->r() -= TINY_BIT * u();
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
    this->r() = r;
  }

  // Adjust the particle's location and direction.
  r() = new_r;
  u() = new_u;

  // Reassign particle's surface
  surface() = new_surface;

  // Figure out what cell particle is in now
  n_coord() = 1;

  if (!neighbor_list_find_cell(*this)) {
    mark_as_lost("Couldn't find particle after hitting periodic "
                 "boundary on surface " +
                 std::to_string(surf.id_) + ".");
    return;
  }

  // Set previous coordinate going slightly past surface crossing
  r_last_current() = r() + TINY_BIT * u();

  // Diagnostic message
  if (settings::verbosity >= 10 || trace()) {
    write_message(1, "    Hit periodic boundary on surface {}", surf.id_);
  }
}

void Particle::mark_as_lost(const char* message)
{
  // Print warning and write lost particle file
  warning(message);
  if (settings::max_write_lost_particles < 0 ||
      simulation::n_lost_particles < settings::max_write_lost_particles) {
    write_restart();
  }
  // Increment number of lost particles
  wgt() = 0.0;
#pragma omp atomic
  simulation::n_lost_particles += 1;

  // Count the total number of simulated particles (on this processor)
  auto n = simulation::current_batch * settings::gen_per_batch *
           simulation::work_per_rank;

  // Abort the simulation if the maximum number of lost particles has been
  // reached
  if (simulation::n_lost_particles >= settings::max_lost_particles &&
      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
    fatal_error("Maximum number of lost particles has been reached.");
  }
}

void Particle::write_restart() const
{
  // Dont write another restart file if in particle restart mode
  if (settings::run_mode == RunMode::PARTICLE)
    return;

  // Set up file name
  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
    simulation::current_batch, id());

#pragma omp critical(WriteParticleRestart)
  {
    // Create file
    hid_t file_id = file_open(filename, 'w');

    // Write filetype and version info
    write_attribute(file_id, "filetype", "particle restart");
    write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);
    write_attribute(file_id, "openmc_version", VERSION);
#ifdef GIT_SHA1
    write_attr_string(file_id, "git_sha1", GIT_SHA1);
#endif

    // Write data to file
    write_dataset(file_id, "current_batch", simulation::current_batch);
    write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
    write_dataset(file_id, "current_generation", simulation::current_gen);
    write_dataset(file_id, "n_particles", settings::n_particles);
    switch (settings::run_mode) {
    case RunMode::FIXED_SOURCE:
      write_dataset(file_id, "run_mode", "fixed source");
      break;
    case RunMode::EIGENVALUE:
      write_dataset(file_id, "run_mode", "eigenvalue");
      break;
    case RunMode::PARTICLE:
      write_dataset(file_id, "run_mode", "particle restart");
      break;
    default:
      break;
    }
    write_dataset(file_id, "id", id());
    write_dataset(file_id, "type", type().pdg_number());

    // Get source site data for the particle that got lost
    int64_t i = current_work();
    SourceSite site;
    if (settings::run_mode == RunMode::EIGENVALUE) {
      site = simulation::source_bank[i];
    } else if (settings::run_mode == RunMode::FIXED_SOURCE &&
               settings::use_shared_secondary_bank &&
               i < simulation::shared_secondary_bank_read.size()) {
      site = simulation::shared_secondary_bank_read[i];
    } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
      // Re-sample using the same seed used to generate the source particle.
      // current_work() is 0-indexed, compute_particle_id expects 1-indexed.
      int64_t id = compute_transport_seed(compute_particle_id(i + 1));
      uint64_t seed = init_seed(id, STREAM_SOURCE);
      site = sample_external_source(&seed);
    }
    write_dataset(file_id, "weight", site.wgt);
    write_dataset(file_id, "energy", site.E);
    write_dataset(file_id, "xyz", site.r);
    write_dataset(file_id, "uvw", site.u);
    write_dataset(file_id, "time", site.time);

    // Close file
    file_close(file_id);
  } // #pragma omp critical
}

void Particle::update_neutron_xs(
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
{
  // Get microscopic cross section cache
  auto& micro = this->neutron_xs(i_nuclide);

  // If the cache doesn't match, recalculate micro xs
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
      ncrystal_xs != micro.ncrystal_xs) {
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);

    // If NCrystal is being used, update micro cross section cache
    micro.ncrystal_xs = ncrystal_xs;
    if (ncrystal_xs >= 0.0) {
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
      ncrystal_update_micro(ncrystal_xs, micro);
    }
  }
}

//==============================================================================
// Non-method functions
//==============================================================================
void add_surf_source_to_bank(Particle& p, const Surface& surf)
{
  if (simulation::current_batch <= settings::n_inactive ||
      simulation::surf_source_bank.full()) {
    return;
  }

  // If a cell/cellfrom/cellto parameter is defined
  if (settings::ssw_cell_id != C_NONE) {

    // Retrieve cell index and storage type
    int cell_idx = model::cell_map[settings::ssw_cell_id];

    if (surf.bc_) {
      // Leave if cellto with vacuum boundary condition
      if (surf.bc_->type() == "vacuum" &&
          settings::ssw_cell_type == SSWCellType::To) {
        return;
      }

      // Leave if other boundary condition than vacuum
      if (surf.bc_->type() != "vacuum") {
        return;
      }
    }

    // Check if the cell of interest has been exited
    bool exited = false;
    for (int i = 0; i < p.n_coord_last(); ++i) {
      if (p.cell_last(i) == cell_idx) {
        exited = true;
      }
    }

    // Check if the cell of interest has been entered
    bool entered = false;
    for (int i = 0; i < p.n_coord(); ++i) {
      if (p.coord(i).cell() == cell_idx) {
        entered = true;
      }
    }

    // Vacuum boundary conditions: return if cell is not exited
    if (surf.bc_) {
      if (surf.bc_->type() == "vacuum" && !exited) {
        return;
      }
    } else {

      // If we both enter and exit the cell of interest
      if (entered && exited) {
        return;
      }

      // If we did not enter nor exit the cell of interest
      if (!entered && !exited) {
        return;
      }

      // If cellfrom and the cell before crossing is not the cell of
      // interest
      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
        return;
      }

      // If cellto and the cell after crossing is not the cell of interest
      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
        return;
      }
    }
  }

  SourceSite site;
  site.r = p.r();
  site.u = p.u();
  site.E = p.E();
  site.time = p.time();
  site.wgt = p.wgt();
  site.delayed_group = p.delayed_group();
  site.surf_id = surf.id_;
  site.particle = p.type();
  site.parent_id = p.id();
  site.progeny_id = p.n_progeny();
  simulation::surf_source_bank.thread_safe_append(site);
}

} // namespace openmc