1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright(C) 2007 Risoe National Laboratory.
*
* %I
* Written by: Mads Bertelsen
* Date: 20.08.15
* Version: $Revision: 0.1 $
* Origin: University of Copenhagen
*
* A sample component to separate geometry and phsysics
*
* %D
*
* 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
* sperates geometry and physics within McStas.
* 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:
* Described elsewhere
*
* %P
* INPUT PARAMETERS:
* sigma: [barns] Incoherent scattering cross section
* f_QE: [1] Fraction of quasielastic scattering (rest is elastic) [1]
* gamma: [meV] Lorentzian width of quasielastic broadening (HWHM) [1]
* 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] name of Union_init component (typically "init", default)
*
* 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="init")
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
#ifndef Union
#error "The Union_init component must be included before this Incoherent_process component"
#endif
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 v_i,v_f,E_i,dE,E_f;
if (rand01() < data_transfer.pointer_to_a_Incoherent_physics_storage_struct->QE_sampling_frequency) {
v_i = k_length * K2V;
E_i = VS2E*v_i*v_i;
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;
v_f = SE2V*sqrt(E_f);
k_length = v_f*V2K;
}
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
%}
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
%{
// 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;
// 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)
// The type of the process must be saved in the global enum process
This_process.eProcess = Incoherent;
// Packing the data into a structure that is transported to the main component
sprintf(This_process.name,"%s",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,"%s",NAME_CURRENT_COMP);
global_process_element.component_index = INDEX_CURRENT_COMP;
global_process_element.p_scattering_process = &This_process;
if (_getcomp_index(init) < 0) {
fprintf(stderr,"Incoherent_process:%s: Error identifying Union_init component, %s is not a known component name.\n",
NAME_CURRENT_COMP, init);
exit(-1);
}
struct pointer_to_global_process_list *global_process_list = COMP_GETPAR3(Union_init, init, global_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
%}
END
|