|
| 1 | +/******************************************************************************* |
| 2 | +* |
| 3 | +* McStas, neutron ray-tracing package |
| 4 | +* Copyright(C) 2007 Risoe National Laboratory. |
| 5 | +* |
| 6 | +* %I |
| 7 | +* Written by: Mads Bertelsen |
| 8 | +* Date: 20.08.15 |
| 9 | +* Version: $Revision: 0.1 $ |
| 10 | +* Origin: University of Copenhagen |
| 11 | +* |
| 12 | +* A Mirror surface process |
| 13 | +* |
| 14 | +* %D |
| 15 | +* |
| 16 | +* This is a Union surface process that describes a supermirror or other surface |
| 17 | +* that only have specular reflection. The reflectivity can be given as a file |
| 18 | +* or using the standard reflectivity inputs. To use this in a simulation an |
| 19 | +* instance of this component should be defined in the instrument file, then |
| 20 | +* attatched to one or more geometries in their surface stacks pertaining to |
| 21 | +* each face of the geometry. |
| 22 | +* |
| 23 | +* Part of the Union components, a set of components that work together and thus |
| 24 | +* sperates geometry and physics within McStas. |
| 25 | +* The use of this component requires other components to be used. |
| 26 | +* |
| 27 | +* 1) One specifies a number of processes using process components like this one |
| 28 | +* 2) These are gathered into material definitions using Union_make_material |
| 29 | +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material |
| 30 | +* 4) A Union_master component placed after all of the above |
| 31 | +* |
| 32 | +* Only in step 4 will any simulation happen, and per default all geometries |
| 33 | +* defined before the master, but after the previous will be simulated here. |
| 34 | +* |
| 35 | +* There is a dedicated manual available for the Union_components |
| 36 | +* |
| 37 | +* |
| 38 | +* Algorithm: |
| 39 | +* Described elsewhere |
| 40 | +* |
| 41 | +* %P |
| 42 | +* INPUT PARAMETERS: |
| 43 | +* R0: [1] Low-angle reflectivity |
| 44 | +* Qc: [AA-1] Critical scattering vector |
| 45 | +* alpha: [AA] Slope of reflectivity |
| 46 | +* m: [1] m-value of material. Zero means completely absorbing. |
| 47 | +* W: [AA-1] Width of supermirror cut-off |
| 48 | +* reflect: [str] Name of reflectivity file. Format q(Angs-1) R(0-1) |
| 49 | +* init: [string] Name of Union_init component (typically "init", default) |
| 50 | +* |
| 51 | +* %L |
| 52 | +* |
| 53 | +* %E |
| 54 | +******************************************************************************/ |
| 55 | + |
| 56 | +DEFINE COMPONENT Mirror_surface // Remember to change the name of process here |
| 57 | + |
| 58 | +SETTING PARAMETERS(string reflect=0, R0=0.99, Qc=0.0219, alpha=6.07, m=2, W=0.003, string init="init") |
| 59 | + |
| 60 | + |
| 61 | +SHARE |
| 62 | +%{ |
| 63 | +#ifndef Union |
| 64 | +#error "The Union_init component must be included before this Mirror_surface component" |
| 65 | +#endif |
| 66 | + |
| 67 | + |
| 68 | +%include "read_table-lib" |
| 69 | +%include "ref-lib" |
| 70 | + |
| 71 | +// Very important to add a pointer to this struct in the union-lib.c file |
| 72 | +struct Mirror_surface_storage_struct{ |
| 73 | + // Variables that needs to be transfered between any of the following places: |
| 74 | + // The initialize in this component |
| 75 | + // The function for calculating my |
| 76 | + // The function for calculating scattering |
| 77 | + |
| 78 | + // Avoid duplicates of setting parameters in naming |
| 79 | + double par[5]; |
| 80 | + t_Table pTable; |
| 81 | + int table_present; |
| 82 | +}; |
| 83 | + |
| 84 | +// Function for handling surface physics |
| 85 | +// The input for this function and its order may not be changed, but the names may be updated. |
| 86 | +int Mirror_surface_function(union surface_data_transfer_union data_transfer, // data in struct defined above, can have data stored in initialize and other calls to this function |
| 87 | + double *weight, double *wavevector, int *continues, // given weight and wavevector, to be updated, and pointer to continues which should be provided by function |
| 88 | + double *normal_vector, enum in_or_out inward_or_outward, // normal_vector of surface and information on whether the neutron is going into or out of this surface |
| 89 | + _class_particle *_particle) { |
| 90 | + |
| 91 | + // wavevector and normal_vector are a pointers to a simple array with 3 doubles, wavevector[0], wavevector[1], wavevector[2] which describes the vector |
| 92 | + |
| 93 | + // Retrieve data storage struct |
| 94 | + struct Mirror_surface_storage_struct *storage; |
| 95 | + storage = data_transfer.pointer_to_a_Mirror_surface_storage_struct; |
| 96 | + |
| 97 | + double k_n_dot = wavevector[0]*normal_vector[0] + wavevector[1]*normal_vector[1] + wavevector[2]*normal_vector[2]; |
| 98 | + |
| 99 | + // Calculate q normal |
| 100 | + double q_normal = 2.0*fabs(k_n_dot); |
| 101 | + |
| 102 | + // Calculate or read reflectivity |
| 103 | + double reflectivity; |
| 104 | + if (storage->table_present == 1) { |
| 105 | + TableReflecFunc(q_normal, &storage->pTable, &reflectivity); |
| 106 | + } else { |
| 107 | + StdReflecFunc(q_normal, storage->par, &reflectivity); |
| 108 | + } |
| 109 | + |
| 110 | + // Take monte carlo choice on whether to be reflected by the mirror or not |
| 111 | + if (rand01() > reflectivity) { |
| 112 | + *continues = 1; // Tells algorithm the ray goes through this surface |
| 113 | + } else { |
| 114 | + // mirror ray |
| 115 | + wavevector[0] = wavevector[0] - 2.0 * k_n_dot * normal_vector[0]; |
| 116 | + wavevector[1] = wavevector[1] - 2.0 * k_n_dot * normal_vector[1]; |
| 117 | + wavevector[2] = wavevector[2] - 2.0 * k_n_dot * normal_vector[2]; |
| 118 | + |
| 119 | + *continues = 0; // Tells algorithm the ray does not go through this surface |
| 120 | + } |
| 121 | + |
| 122 | + return 1; |
| 123 | + |
| 124 | + // There is access to the data_transfer from within this function: |
| 125 | + // int table_present = data_transfer.pointer_to_a_Mirror_physics_storage_struct->table_present; |
| 126 | + // Its even possible to write to this data structure and use results for a later ray if relevant |
| 127 | +}; |
| 128 | + |
| 129 | +// These lines help with future error correction, and tell other Union components |
| 130 | +// that at least one surface have been defined. |
| 131 | +#ifndef SURFACE_DETECTOR |
| 132 | + #define SURFACE_DETECTOR dummy |
| 133 | +#endif |
| 134 | +#ifndef SURFACE_PROCESS_MIRROR_DETECTOR |
| 135 | + #define SURFACE_PROCESS_MIRROR_DETECTOR dummy |
| 136 | +#endif |
| 137 | +%} |
| 138 | + |
| 139 | +DECLARE |
| 140 | +%{ |
| 141 | +// Declare for this component, to do calculations on the input / store in the transported data |
| 142 | +struct Mirror_surface_storage_struct Mirror_storage; // Replace Mirror with your own name here |
| 143 | + |
| 144 | +// Needed for transport to the main component, will be the same for all surface processes |
| 145 | +struct global_surface_element_struct global_surface_element; |
| 146 | +struct surface_process_struct This_surface; |
| 147 | +%} |
| 148 | + |
| 149 | +INITIALIZE |
| 150 | +%{ |
| 151 | + // Initialize done in the component |
| 152 | + if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) { |
| 153 | + if (Table_Read(&Mirror_storage.pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */ |
| 154 | + exit(fprintf(stderr,"Mirror_surface: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect)); |
| 155 | + Mirror_storage.table_present=1; |
| 156 | + } else { |
| 157 | + Mirror_storage.table_present=0; |
| 158 | + if (W < 0 || R0 < 0 || Qc < 0 || m < 0) { |
| 159 | + fprintf(stderr,"Mirror_surface: %s: W R0 Qc must be >0.\n", NAME_CURRENT_COMP); |
| 160 | + exit(-1); |
| 161 | + } |
| 162 | + } |
| 163 | + |
| 164 | + // Insert input parameter to par array as required for the standard reflectivity function |
| 165 | + Mirror_storage.par[0] = R0; |
| 166 | + Mirror_storage.par[1] = Qc; |
| 167 | + Mirror_storage.par[2] = alpha; |
| 168 | + Mirror_storage.par[3] = m; |
| 169 | + Mirror_storage.par[4] = W; |
| 170 | + |
| 171 | + // Will need to do similar system for surfaces |
| 172 | + // Need to specify if this process is isotropic |
| 173 | + //This_process.non_isotropic_rot_index = -1; // Yes (powder) |
| 174 | + //This_process.non_isotropic_rot_index = 1; // No (single crystal) |
| 175 | + //rot_copy(This_surface.rotation_matrix, ROT_A_CURRENT_COMP); |
| 176 | + |
| 177 | + // The type of the process must be saved in the global enum process located in union-lib.c |
| 178 | + This_surface.eSurface = Mirror; |
| 179 | + |
| 180 | + // Packing the data into a structure that is transported to the main component |
| 181 | + This_surface.data_transfer.pointer_to_a_Mirror_surface_storage_struct = &Mirror_storage; |
| 182 | + |
| 183 | + // This will be the same for all process's, and can thus be moved to an include. |
| 184 | + sprintf(This_surface.name,"%s",NAME_CURRENT_COMP); |
| 185 | + sprintf(global_surface_element.name,"%s",NAME_CURRENT_COMP); |
| 186 | + global_surface_element.component_index = INDEX_CURRENT_COMP; |
| 187 | + global_surface_element.p_surface_process = &This_surface; |
| 188 | + |
| 189 | + if (_getcomp_index(init) < 0) { |
| 190 | + fprintf(stderr,"Mirror_surface:%s: Error identifying Union_init component, %s is not a known component name.\n", |
| 191 | + NAME_CURRENT_COMP, init); |
| 192 | + exit(-1); |
| 193 | + } |
| 194 | + |
| 195 | + struct pointer_to_global_surface_list *global_surface_list = COMP_GETPAR3(Union_init, init, global_surface_list); |
| 196 | + add_element_to_surface_list(global_surface_list, global_surface_element); |
| 197 | + %} |
| 198 | + |
| 199 | +TRACE |
| 200 | +%{ |
| 201 | + // Trace should be empty, the simulation is perfomed in Union_master |
| 202 | +%} |
| 203 | + |
| 204 | +END |
| 205 | + |
0 commit comments