eimpact.cpp

#include "Matrix.hpp"
#include "PhysicsConstants.hpp"
#include "Scattering.hpp"
#include <fstream>
#include <iostream>
#include <string>

const std::string usage_text{R"(
Calculates electron-impact ionisation cross-section.

Usage:

Takes input directly from command line.
First argument is file containing Kion.
(File generated by ampsci in 'mat' format, assumed to be in atomic units).

e.g.,
  ./eimpact K_Xe_v_6_hp_orth_mat.txt > sigma_e_Xe.txt

sigma_e_Xe.txt will contain the cross-section, ready for plotting.

)"};

//==============================================================================
//==============================================================================
//==============================================================================
int main(int argc, char *argv[]) {

  // read input filename from command line:
  const std::string in_filename = argc > 1 ? argv[1] : "";
  if (in_filename.empty()) {
    std::cout << "Error: No file given (pass as command-line arg)\n";
    std::cout << usage_text;
    return 1;
  }

  // Read in Kion (and Energy/momentum grids)
  const auto [E_grid, q_grid, Kion] = Scattering::read_in_Kion(in_filename);

  // Check formats correct:
  assert(E_grid.size() == Kion.rows() && "Input file must be in valid format");
  assert(q_grid.size() == Kion.cols() && "Input file must be in valid format");

  // Calculate sig(E) for e-e impact ionisation, in atomic units
  const auto sig_impact = Scattering::sigtot(Kion, E_grid, q_grid);

  // Sanity check:
  assert(sig_impact.size() == E_grid.size());

  // Convert units: energy to keV, and sigma to cm^2
  using namespace VectorOverloads;
  const auto E_keV = E_grid * Physics::Conversions::E_to_keV;
  constexpr auto a02 = Physics::Constants::aB_cm * Physics::Constants::aB_cm;
  const auto sigma_cm2 = sig_impact * a02;

  std::cout
      << "# outputs sigma(E), where E is incident electron energy (in keV), "
         "and sigma is the total ionisation cross-section, in cm^2\n";
  std::cout << "E(keV) sigma(cm^2)";
  for (std::size_t i = 0; i < sig_impact.size(); ++i) {
    std::cout << E_keV[i] << " " << sigma_cm2[i] << "\n";
  }
}