/*******************************************************************************
*
* McStas, neutron ray-tracing package
*
* Component: MCViNE_Phonon_CoherentInelastic_PolyXtal_process
*
* %I
* Written by: Fahima Islam (McStas port of the MCViNE phonon CoherentInelastic_PolyXtal kernel; MCViNE by J. Y. Y. Lin et al.)
* Date: 2026-09-29
* Origin: MCViNE, https://github.com/mcvine/mcvine (mccomponents/lib/kernels/sample/phonon/CoherentInelastic_PolyXtal.cc)
*
* Union process: Coherent one-phonon scattering from a powder, using a full phonon dispersion (energies + polarizations) on a grid.
*
* %D
* Powder average of the coherent one-phonon cross section, computed from phonon
* energies and polarization vectors tabulated on a grid over one reciprocal cell
* (MCViNE IDF format, e.g. from phonopy via MCViNE tools). A random branch and a
* random Q vector in a cube are drawn until the event is kinematically allowed;
* phonon creation/annihilation, Bose factor and Debye-Waller factor included.
* Scattering coefficient: total coherent cross section / Vc.
*
* <b>Union process.</b> Part of the Union components: define this process, collect it
* into a material with Union_make_material (absorption is set there with
* my_absorption, the absorption inverse penetration depth at 2200 m/s), assign
* the material to Union_box/Union_cylinder/Union_sphere/Union_mesh geometries,
* and add a Union_master after them. Geometry, attenuation and multiple
* scattering are handled by Union_master; this component provides the MCViNE
* kernel: the scattering coefficient and the final-state sampling (weight).
* Isotropic process (powder/liquid-like); rotation is irrelevant.
* The kernel code is shared with the standalone component MCViNE_Phonon_CoherentInelastic_PolyXtal (mcvine-lib.c).
* Uses share/mcvine-lib.h/.c and share/mcvine-union-lib.h/.c; the process type
* MCViNE is declared in share/union-lib.h and registered with Union_master in
* share/mcvine-union-lib.h.
*
* Example: MCViNE_Phonon_CoherentInelastic_PolyXtal_process(idf_dir="MCViNE/fcc_toy_phonons", atoms="MCViNE/fcc_toy_atoms.dat", T=300, max_omega=45)
*
* %P
* INPUT PARAMETERS:
* idf_dir: [str] Directory with MCViNE IDF phonon files: Qgridinfo, Omega2, Polarizations[, DOS]
* atoms: [str] Atoms file: one row per atom 'x y z mass b_coh sigma_inc sigma_abs' ([AA] cartesian, [amu], [fm], [barn], [barn]), same order as in the IDF files
* T: [K] Temperature
* dw_core: [AA^2] Debye-Waller core; <0: computed from the DOS (idf_dir/DOS or dos)
* dos: [str] Optional DOS file for the Debye-Waller factor (default: idf_dir/DOS)
* Vc: [AA^3] Unit cell volume; 0: (2pi)^3/|b1.(b2 x b3)| from Qgridinfo
* max_omega: [meV] Maximum phonon energy
* min_omega: [meV] Phonons below this energy are skipped
* unbiased: [1] 0: MCViNE sampling (rejection + empirical accessible reciprocal volume, a few % high in tests). 1: single Q sample in the cube with its exact volume, unbiased
* packing_factor: [1] Packing factor (scales the scattering coefficient)
* interact_fraction: [1] Union: fraction of interactions forced to this process (-1: by cross section)
* init: [string] Deprecated and unused, accepted like on the other Union components (see Union_init).
*
* %L
* MCViNE documentation: https://mcvine.github.io
*
* %E
*******************************************************************************/

DEFINE COMPONENT MCViNE_Phonon_CoherentInelastic_PolyXtal_process

SETTING PARAMETERS (string idf_dir=0, string atoms=0, T=300, dw_core=-1, string dos=0, Vc=0, max_omega=50, min_omega=0.01, int unbiased=0, packing_factor=1, interact_fraction=-1, string init="")

NOACC

SHARE
%{
  %include "union-lib"
  %include "read_table-lib"
  %include "mcvine-lib"
  %include "mcvine-union-lib"
  #ifndef PROCESS_DETECTOR
  #define PROCESS_DETECTOR dummy
  #endif
%}

DECLARE
%{
  mcvine_kernel_Phonon_CoherentInelastic_PolyXtal kernel;
  struct MCViNE_physics_storage_struct storage;
  struct scattering_process_struct This_process;
  struct global_process_element_struct global_process_element;
%}

INITIALIZE
%{
  struct union_state_struct* union_state_p = union_acquire ();
  double mu = 0, sig = 0;
  memset (&kernel, 0, sizeof (kernel));
  {
    int i;
    double mass = 0, xs_abs = 0, V;
    mcvine_dispersion* dp = (mcvine_dispersion*)calloc (1, sizeof (mcvine_dispersion));
    if (mcvine_dispersion_load_idf (dp, idf_dir, NAME_CURRENT_COMP))
      exit (-1);
    kernel.m_natoms = mcvine_atoms_load (&kernel.m_atoms, atoms, NAME_CURRENT_COMP);
    if (kernel.m_natoms < 0)
      exit (-1);
    if (kernel.m_natoms != dp->natoms) {
      fprintf (stderr, "%s: atoms file has %d atoms but the dispersion has %d\n", NAME_CURRENT_COMP, kernel.m_natoms, dp->natoms);
      exit (-1);
    }
    kernel.m_disp = dp;
    kernel.m_T = T;
    kernel.m_xs_coh_tot = 0;
    for (i = 0; i < kernel.m_natoms; i++) {
      mass += kernel.m_atoms[i].mass;
      kernel.m_xs_coh_tot += kernel.m_atoms[i].xs_coh;
      xs_abs += kernel.m_atoms[i].xs_abs;
    }
    mass /= kernel.m_natoms;
    if (dw_core >= 0)
      kernel.m_dw_core = dw_core;
    else {
      mcvine_dos d;
      if (dos && dos[0]) {
        if (mcvine_dos_load (&d, dos, 0, NAME_CURRENT_COMP))
          exit (-1);
      } else if (dp->has_dos)
        d = dp->m_dos;
      else {
        fprintf (stderr, "%s: no DOS for the Debye-Waller factor (give dw_core, dos, or idf_dir/DOS)\n", NAME_CURRENT_COMP);
        exit (-1);
      }
      kernel.m_dw_core = mcvine_dw_core_from_dos (&d, mass, T, 100);
    }
    V = Vc > 0 ? Vc : dp->ucvol;
    mu = mcvine_xs2coeff (xs_abs, V);
    sig = mcvine_xs2coeff (kernel.m_xs_coh_tot, V);
    printf ("%s: %d atoms, %d branches, grid %dx%dx%d, Vc=%g AA^3, sigma_coh=%g barn, DW core=%g AA^2\n", NAME_CURRENT_COMP, kernel.m_natoms, dp->nbranches,
            dp->n[0], dp->n[1], dp->n[2], V, kernel.m_xs_coh_tot, kernel.m_dw_core);
  }
  kernel.m_max_omega = max_omega;
  kernel.m_min_omega = min_omega;
  kernel.m_unbiased = unbiased;
  printf ("%s: absorption is not part of the process; set Union_make_material(my_absorption=%g) for this material\n", NAME_CURRENT_COMP, mu);
  storage.m_kernel = &kernel;
  storage.m_S = mcvine_S_Phonon_CoherentInelastic_PolyXtal;
  storage.m_my_scattering = sig * packing_factor;
  storage.m_kind = MCVINE_UNION_GENERIC;
  mcvine_union_register (&union_state_p->u_process_list, &This_process, &global_process_element, &storage, NAME_CURRENT_COMP, INDEX_CURRENT_COMP, interact_fraction, 0, ROT_A_CURRENT_COMP);
%}

TRACE
%{
  // the simulation is done in Union_master
%}

FINALLY
%{
  union_release (NAME_CURRENT_COMP);
%}

END