C.1  Run-time calls and functions (mcstas-r)

Here we list a number of preprogrammed macros which may ease the task of writing component and instrument definitions.

C.1.1  Neutron propagation

Propagation routines perform all necessary operations to transport neutron rays from one point to an other. Except when using the special ALLOW_BACKPROP; call prior to executing any PROP_* propagation, the neutron rays which have negative propagation times are removed automatically.

• ABSORB. This macro issues an order to the overall McStas simulator to interrupt the simulation of the current neutron history and to start a new one.

• PROP_Z0. Propagates the neutron to the \(z=0\) plane, by adjusting \((x,y,z)\) and \(t\) accordingly from knowledge of the neutron velocity \((vx,vy,vz)\). If the propagation time is negative, the neutron ray is absorbed, except if a ALLOW_BACKPROP; preceeds it.

  For components that are centered along the \(z\)-axis, use the _intersect functions to determine intersection time(s), and then a PROP_DT call.

• PROP_DT\((dt)\). Propagates the neutron through the time interval \(dt\), adjusting \((x,y,z)\) and \(t\) accordingly from knowledge of the neutron velocity. This macro automatically calls PROP_GRAV_DT when the --gravitation option has been set for the whole simulation.

• PROP_GRAV_DT\((dt,Ax,Ay,Az)\). Like PROP_DT, but it also includes gravity using the acceleration \((Ax,Ay,Az)\). In addition to adjusting \((x,y,z)\) and \(t\), also \((vx,vy,vz)\) is modified.

• PROP_NO_G_DT\((dt)\). Propagates the neutron through the time interval \(dt\) using straight-line motion, regardless of the global gravitation setting. The negative-time and back-propagation rules are the same as for PROP_DT.

• PROP_NO_G_X0, PROP_NO_G_Y0, and PROP_NO_G_Z0. Propagate the neutron to the corresponding local coordinate plane using straight-line motion, regardless of the global gravitation setting.

• ALLOW_BACKPROP. Indicates that the next propagation routine will not remove the neutron ray, even if negative propagation times are found. Subsequent propagations are not affected.

• SCATTER. This macro is used to denote a scattering event inside a component. It should be used e.g to indicate that a component has interacted with the neutron ray (e.g. scattered or detected). This does not affect the simulation (see, however, Beamstop), and it is mainly used by the MCDISPLAY section and the GROUP modifier See also the SCATTERED variable (below).

C.1.2  Coordinate and component variable retrieval

• MC_GETPAR\((comp, outpar)\). This may be used in e.g. the FINALLY section of an instrument definition to reference the output parameters of a component.

• NAME_CURRENT_COMP gives the name of the current component as a string.

• POS_A_CURRENT_COMP gives the absolute position of the current component. A component of the vector is referred to as POS_A_CURRENT_COMP.\(i\) where \(i\) is \(x\), \(y\) or \(z\).

• ROT_A_CURRENT_COMP and ROT_R_CURRENT_COMP give the orientation of the current component as rotation matrices (absolute orientation and the orientation relative to the previous component, respectively). A component of a rotation matrix is referred to as ROT_A_CURRENT_COMP\([m][n]\), where \(m\) and \(n\) are 0, 1, or 2 standing for \(x,y\) and \(z\) coordinates respectively.

• POS_A_COMP\((comp)\) gives the absolute position of the component with the name comp. Note that comp is not given as a string. A component of the vector is referred to as POS_A_COMP\((comp).i\) where \(i\) is \(x\), \(y\) or \(z\).

• ROT_A_COMP\((comp)\) and ROT_R_COMP\((comp)\) give the orientation of the component comp as rotation matrices (absolute orientation and the orientation relative to its previous component, respectively). Note that comp is not given as a string. A component of a rotation matrice is referred to as ROT_A_COMP\((comp)[m][n]\), where \(m\) and \(n\) are 0, 1, or 2.

• INDEX_CURRENT_COMP is the number (index) of the current component (starting from 1).

• POS_A_COMP_INDEX\((index)\) is the absolute position of component \(index\).
  POS_A_COMP_INDEX (INDEX_CURRENT_COMP) is the same as
  POS_A_CURRENT_COMP. You may use
  POS_A_COMP_INDEX (INDEX_CURRENT_COMP+1)
  to make, for instance, your component access the position of the next component (this is usefull for automatic targeting). A component of the vector is referred to as POS_A_COMP_INDEX\((index).i\) where \(i\) is \(x\), \(y\) or \(z\).

• POS_R_COMP_INDEX works the same as above, but with relative coordinates.

• STORE_NEUTRON\((index, x, y, z, vx, vy, vz, t, sx, sy, sz, p)\) stores the current neutron state in the trace-history table, in local coordinate system. \(index\) is usually INDEX_CURRENT_COMP. This is automatically done when entering each component of an instrument.

• RESTORE_NEUTRON\((index, x, y, z, vx, vy, vz, t, sx, sy, sz, p)\) restores the neutron state to the one at the input of the component \(index\). To ignore a component effect, use RESTORE_NEUTRON (INDEX_CURRENT_COMP,
  \(x, y, z, vx, vy, vz, t, sx, sy, sz, p\)) at the end of its TRACE section, or in its EXTEND section. These neutron states are in the local component coordinate systems.

• SCATTERED is a variable set to 0 when entering a component, which is incremented each time a SCATTER event occurs. This may be used in the EXTEND sections to determine whether the component interacted with the current neutron ray.

• extend_list(\(n\), &arr, &len, elemsize). Given an array arr with len elements each of size elemsize, make sure that the array is big enough to hold at least \(n\) elements, by extending arr and len if necessary. Typically used when reading a list of numbers from a data file when the length of the file is not known in advance.

• mcset_ncount\((n)\). Sets the number of neutron histories to simulate to \(n\).

• mcget_ncount(). Returns the number of neutron histories to simulate (usually set by option -n).

• mcget_run_num(). Returns the number of neutron histories that have been simulated until now.

C.1.3  Coordinate transformations

• coords_set\((x,y,z)\) returns a Coord structure (like POS_A_CURRENT_COMP) with \(x\), \(y\) and \(z\) members.

• coords_get\((P,\) &\(x\), &\(y\), &\(z)\) copies the \(x\), \(y\) and \(z\) members of the Coord structure \(P\) into \(x,y,z\) variables.

• coords_add\((a,b)\), coords_sub\((a,b)\), coords_neg\((a)\) enable to operate on coordinates, and return the resulting Coord structure.

• rot_set_rotation(Rotation t, \(\phi _x, \phi _y, \phi _z\)) Get transformation matrix for rotation first \(\phi _x\) around x axis, then \(\phi _y\) around y, and last \(\phi _z\) around z. \(t\) should be a ’Rotation’ ([3][3] ’double’ matrix).

• rot_mul(Rotation t1, Rotation t2, Rotation t3) performs \(t3 = t1 . t2\).

• rot_copy(Rotation dest, Rotation src) performs \(dest = src\) for Rotation arrays.

• rot_transpose(Rotation src, Rotation dest) performs \(dest = src^t\).

• rot_apply(Rotation t, Coords a) returns a Coord structure which is \(t.a\)

C.1.4  Mathematical routines

• NORM\((x,y,z)\). Normalizes the vector \((x,y,z)\) to have length 1 \(^*\).

• scalar_prod\((a_x,a_y,a_z, b_x,b_y,b_z)\). Returns the scalar product of the two vectors \((a_x,a_y,a_z)\) and \((b_x,b_y,b_z)\).

• vec_prod(\(a_x\),\(a_y\),\(a_z\), \(b_x\),\(b_y\),\(b_z\), \(c_x\),\(c_y\),\(c_z\)). Sets \((a_x,a_y,a_z)\) equal to the vector product \((b_x,b_y,b_z) \times (c_x,c_y,c_z)\) \(^*\).

• rotate(\(x\),\(y\),\(z\), \(v_x\),\(v_y\),\(v_z\), \(\varphi \), \(a_x\),\(a_y\),\(a_z\)). Set \((x,y,z)\) to the result of rotating the vector \((v_x,v_y,v_z)\) the angle \(\varphi \) (in radians) around the vector \((a_x,a_y,a_z)\) \(^*\).

• normal_vec(\(n_x\), \(n_y\), \(n_z\), \(x\), \(y\), \(z\)). Computes a unit vector \((n_x, n_y, n_z)\) normal to the vector \((x,y,z)\).\(^*\)

• solve_2nd_order(&\(t_1\), &\(t_2\) \(A\), \(B\), \(C\)). Solves the 2\(^{nd}\) order equation \(At^2 + Bt + C = 0\) and returns the solutions into pointers *\(t_1\) and *\(t_2\). if \(t_2\) is specified as NULL, only the smallest positive solution is returned in \(t_1\).

(\(^*\) The experienced c-programmer may be puzzled that these routines can return information without the use of pass by reference, the reason is that these calls are implemented as macros / #define wrapped functions.)

C.1.5  Output from detectors

Details about using these functions are given in the McStas User Manual.

• DETECTOR_OUT_0D\((...)\). Used to output the results from a single detector. The name of the detector is output together with the simulated intensity and estimated statistical error. The output is produced in a format that can be read by McStas front-end programs.

• DETECTOR_OUT_1D\((...)\). Used to output the results from a one-dimensional detector. Integrated intensities error etc. is also reported as for DETECTOR_OUT_0D.

• DETECTOR_OUT_2D\((...)\). Used to output the results from a two-dimentional detector. Integrated intensities error etc. is also reported as for DETECTOR_OUT_0D.

• DETECTOR_OUT_3D\((...)\). Used to output the results from a three-dimentional detector. Arguments are the same as in DETECTOR_OUT_2D, but with an additional \(z\) axis. Resulting data files are treated as 2D data, but the 3rd dimension is specified in the \(type\) field. Integrated intensities error etc. is also reported as for DETECTOR_OUT_0D.

• mcinfo_simulation(FILE *f, mcformat, char *pre, char *name) is used to append the simulation parameters into file \(f\) (see for instance Res_monitor). Internal variable \(mcformat\) should be used as specified. Please contact the authors for further information.

C.1.6  Ray-geometry intersections

• inside_rectangle(\(x\), \(y\), \(xw\), \(yh\)). Return 1 if \(-xw/2 \leq x \leq xw/2\) AND \(-yh/2 \leq y \leq yh/2\). Else return 0.

• box_intersect(&\(t_1\), &\(t_2\), \(x\), \(y\), \(z\), \(v_x\), \(v_y\), \(v_z\), \(d_x\), \(d_y\), \(d_z\)). Calculates the (0, 1, or 2) intersections between the neutron path and a box of dimensions \(d_x\), \(d_y\), and \(d_z\), centered at the origin for a neutron with the parameters \((x,y,z,v_x,v_y,v_z)\). The times of intersection are returned in the variables \(t_1\) and \(t_2\), with \(t_1 < t_2\). In the case of less than two intersections, \(t_1\) (and possibly \(t_2\)) are set to zero. The function returns true if the neutron intersects the box, false otherwise.

• cylinder_intersect(&\(t_1\), &\(t_2\), \(x\), \(y\), \(z\), \(v_x\), \(v_y\), \(v_z\), \(r\), \(h\)). Similar to box_intersect, but using a cylinder of height \(h\) and radius \(r\), centered at the origin.

• sphere_intersect(&\(t_1\), &\(t_2\), \(x\), \(y\), \(z\), \(v_x\), \(v_y\), \(v_z\), \(r\)). Similar to box_intersect, but using a sphere of radius \(r\).

C.1.7  Random numbers

• rand01(). Returns a random number distributed uniformly between 0 and 1.

• randnorm(). Returns a random number from a normal distribution centered around 0 and with \(\sigma =1\). The algorithm used to sample the normal distribution is explained in Ref. [Pre+86, ch.7].

• randpm1(). Returns a random number distributed uniformly between -1 and 1.

• randtriangle(). Returns a random number from a triangular distribution between -1 and 1.

• randvec_target_circle(&\(v_x\), &\(v_y\), &\(v_z\), &\(d\Omega \), aim\(_x\), aim\(_y\), aim\(_z\), \(r_f\)). Generates a random vector \((v_x, v_y, v_z)\), of the same length as (aim\(_x\), aim\(_y\), aim\(_z\)), which is targeted at a disk centered at (aim\(_x\), aim\(_y\), aim\(_z\)) with radius \(r_f\) (in meters), and perpendicular to the aim vector.. All directions that intersect the circle are chosen with equal probability. The solid angle of the circle as seen from the position of the neutron is returned in \(d\Omega \). This routine was previously called randvec_target_sphere (which still works).

• randvec_target_rect_angular(&\(v_x\), &\(v_y\), &\(v_z\), &\(d\Omega \), aim\(_x\), aim\(_y\), aim\(_z\),\(h, w, Rot\)) does the same as randvec_target_circle but targetting at a rectangle with angular dimensions \(h\) and \(w\) (in radians, not in degrees as other angles). The rotation matrix \(Rot\) is the coordinate system orientation in the absolute frame, usually ROT_A_CURRENT_COMP.

• randvec_target_rect(&\(v_x\), &\(v_y\), &\(v_z\), &\(d\Omega \), aim\(_x\), aim\(_y\), aim\(_z\),\(height, width, Rot\)) is the same as randvec_target_rect_angular but \(height\) and \(width\) dimensions are given in meters. This function is useful to e.g. target at a guide entry window or analyzer blade.