/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* %Identification
* Written by: Mads Bertelsen and Erik B Knudsen
* Date: 20.08.15
* Version: $Revision: 0.1 $
* Origin: ESS DMSC & DTU Physics
*
* Component implementing a true incoherent scattering process 
*
* %Description
*
* This Union_process is based on the Incoherent.comp component originally written
*  by Kim Lefmann and Kristian Nielsen.
*
* Part of the Union components, a set of components that work together and thus
*  separates geometry and physics within McXtrace.
* The use of this component requires other components to be used.
*
* 1) One specifies a number of processes using process components like this one
* 2) These are gathered into material definitions using Union_make_material
* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material
* 4) A Union_master component placed after all of the above
*
* Only in step 4 will any simulation happen, and per default all geometries
*  defined before the master, but after the previous will be simulated here.
*
* There is a dedicated manual available for the Union_components
*
* Algorithm:
* This component is an approximation. It scatters completely isotropically.
*
* %Parameters
* INPUT PARAMETERS:
* sigma: [barns]          Incoherent scattering cross section
* f_QE: [1]               Fraction of quasielastic scattering (rest is elastic)
* gamma: [1]              Lorentzian width of quasielastic broadening (HWHM)
* packing_factor: [1]      How dense is the material compared to optimal 0-1
* unit_cell_volume: [AA^3] Unit_cell_volume
* interact_fraction: [1]   How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1)
* init:              [string]    Deprecated and unused. Accepted so that instruments written for McStas/McXtrace 3.8.7 and earlier, which name the Union_init component here, still compile.
*
* CALCULATED PARAMETERS:
*
* %L
* The test/example instrument <a href="../examples/Test_Phonon.instr">Test_Phonon.instr</a>.
*
* %E
******************************************************************************/

DEFINE COMPONENT Incoherent_process
SETTING PARAMETERS(sigma=5.08, f_QE=0, gamma=0, packing_factor=1, unit_cell_volume=13.8, interact_fraction=-1, string init="")

SHARE
%{
  %include "union-lib"

  struct Incoherent_physics_storage_struct {
    // Variables that needs to be transfered between any of the following places:
    // The initialize in this component
    // The function for calculating my
    // The function for calculating scattering

    double my_scattering;
    double QE_sampling_frequency;
    double lorentzian_width;
  };

  // Function for calculating my in Incoherent case
  int
  Incoherent_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data,
                         _class_particle* _particle) {
    *my = data_transfer.pointer_to_a_Incoherent_physics_storage_struct->my_scattering;
    return 1;
  };

  // Function for basic incoherent scattering event
  int
  Incoherent_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer,
                                 struct focus_data_struct* focus_data, _class_particle* _particle) {

    // New version of incoherent scattering
    double k_length = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]);

    Coords k_out;
    // Here is the focusing system in action, get a vector
    double solid_angle;
    focus_data->focusing_function (&k_out, &solid_angle, focus_data);
    NORM (k_out.x, k_out.y, k_out.z);
    *weight *= solid_angle * 0.25 / PI;

    double E_i, dE, E_f;

    if (rand01 () < data_transfer.pointer_to_a_Incoherent_physics_storage_struct->QE_sampling_frequency) {
      E_i = K2E * k_length;
      dE = data_transfer.pointer_to_a_Incoherent_physics_storage_struct->lorentzian_width * tan (PI / 2 * randpm1 ());
      E_f = E_i + dE;
      if (E_f <= 0)
        return 0;
      k_length = E_f * E2K;
    }

    k_final[0] = k_out.x * k_length;
    k_final[1] = k_out.y * k_length;
    k_final[2] = k_out.z * k_length;
    return 1;
  };

  #ifndef PROCESS_DETECTOR
  #define PROCESS_DETECTOR dummy
  #endif

  #ifndef PROCESS_INCOHERENT_DETECTOR
  #define PROCESS_INCOHERENT_DETECTOR dummy
  #endif

  // Register this process with the Union_master dispatch (see union-lib.h)
  #undef UNION_CASE_PHYSICS_MY_INCOHERENT
  #define UNION_CASE_PHYSICS_MY_INCOHERENT(out, ...) case Incoherent: out = Incoherent_physics_my(__VA_ARGS__); break;
  #undef UNION_CASE_PHYSICS_SCATTERING_INCOHERENT
  #define UNION_CASE_PHYSICS_SCATTERING_INCOHERENT(out, ...) case Incoherent: out = Incoherent_physics_scattering(__VA_ARGS__); break;
%}

DECLARE
%{
  // Needed for transport to the main component
  struct global_process_element_struct global_process_element;
  struct scattering_process_struct This_process;

  // Declare for this component, to do calculations on the input / store in the transported data
  struct Incoherent_physics_storage_struct Incoherent_storage;
  double effective_my_scattering;
%}

INITIALIZE
%{
  struct union_state_struct* union_state_p = union_acquire ();
  // Initialize done in the component
  effective_my_scattering = ((packing_factor / unit_cell_volume) * 100 * sigma);
  Incoherent_storage.my_scattering = effective_my_scattering;

  Incoherent_storage.QE_sampling_frequency = f_QE;
  Incoherent_storage.lorentzian_width = gamma;

  // The type of the process must be saved in the global enum process
  This_process.eProcess = Incoherent;

  // Need to specify if this process is isotropic
  This_process.non_isotropic_rot_index = -1; // Yes (powder)
  // This_process.non_isotropic_rot_index =  1;  // No (single crystal)

  // Packing the data into a structure that is transported to the main component
  sprintf (This_process.name, NAME_CURRENT_COMP);
  This_process.process_p_interact = interact_fraction;
  This_process.data_transfer.pointer_to_a_Incoherent_physics_storage_struct = &Incoherent_storage;
  // This_process.data_transfer.pointer_to_a_Incoherent_physics_storage_struct->my_scattering = effective_my_scattering;
  This_process.probability_for_scattering_function = &Incoherent_physics_my;
  This_process.scattering_function = &Incoherent_physics_scattering;

  // This will be the same for all process's, and can thus be moved to an include.
  sprintf (global_process_element.name, NAME_CURRENT_COMP);
  global_process_element.component_index = INDEX_CURRENT_COMP;
  global_process_element.p_scattering_process = &This_process;

  struct pointer_to_global_process_list* global_process_list = &union_state_p->u_process_list;
  add_element_to_process_list (global_process_list, global_process_element);
 %}

TRACE
%{
%}

FINALLY
%{
  // Since the process and it's storage is a static allocation, there is nothing to deallocate
  union_release (NAME_CURRENT_COMP);
%}

END