/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright(C) 2007 Risoe National Laboratory. * * %I * Written by: Mads Bertelsen * Date: 20.08.15 * Origin: University of Copenhagen * * Cylinder geometry component for Union components * * %D * 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 * 2) These are gathered into material definitions using Union_make_material * 3) Geometries are placed using Union_box/cylinder/sphere, 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 this master, but after the previous will be simulated here. * * There is a dedicated manual available for the Union_components * * The position of this component is the center of the cylinder, and it thus * extends yheight/2 up and down along y axis. * * It is allowed to overlap components, but it is not allowed to have two * parallel planes that coincides. This will crash the code on run time. * * * %P * INPUT PARAMETERS: * radius: [m] Outer radius volume in (x,z) plane * yheight: [m] Cylinder height in (y) direction * material_string: [string] Material name of this volume, defined using Union_make_material * priority: [1] Priotiry of the volume (can not be the same as another volume) A high priority is on top of low. * p_interact: [1] Probability to interact with this geometry [0-1] * visualize: [1] Set to 0 if you wish to hide this geometry in mcdisplay * number_of_activations: [1] Number of subsequent Union_master components that will simulate this geometry * mask_string: [string] Comma seperated list of geometry names which this geometry should mask * mask_setting: [string] "All" or "Any", should the masked volume be simulated when the ray is in just one mask, or all. * target_index: [1] Focuses on component a component this many steps further in the component sequence * target_x: [m] X Position of target to focus at * target_y: [m] Y Position of target to focus at * target_z: [m] Z Position of target to focus at * focus_aw: [deg] Horiz. angular dimension of a rectangular area * focus_ah: [deg] Vert. angular dimension of a rectangular area * focus_xw: [m] Horiz. dimension of a rectangular area * focus_xh: [m] Vert. dimension of a rectangular area * focus_r: [m] Focusing on circle with this radius * curved_surface: [string] Comma seperated string of Union surface definitions defined prior to this component, applied to curved wall of cylinder * top_surface: [string] Comma seperated string of Union surface definitions defined prior to this component, applied to +y top face of cylinder * bottom_surface: [string] Comma seperated string of Union surface definitions defined prior to this component, applied to -y bottom face of cylinder * all_face_surface: [string] Comma seperated string of Union surface definitions defined prior to this component, applied to all faces above * cut_surface: [string] Comma seperated string of Union surface definitions defined prior to this component, applied to all internal cuts * 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 * * %E ******************************************************************************/ DEFINE COMPONENT Union_cylinder SETTING PARAMETERS(string material_string=0, priority, radius, yheight, visualize=1, int target_index=0, target_x=0, target_y=0, target_z=0, focus_aw=0, focus_ah=0, focus_xw=0, focus_xh=0, focus_r=0, p_interact=0, string mask_string=0, string mask_setting=0, number_of_activations=1, string curved_surface=0, string top_surface=0, string bottom_surface=0, string all_face_surface=0, string cut_surface=0, string init="" ) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ %include "union-lib" void mcdisplay_cylinder_function (struct lines_to_draw* lines_to_draw_output, int index, struct geometry_struct** Geometries, int number_of_volumes) { // Function to call in mcdisplay section of the sample component for this volume // One can assume that Geometries[index] refers to a geometry as described in this file // The 4 lines describin the cylinders side are aligned to the local frame of the cylinder, // it would be nicer to have them alligned with the global frame so that they show up nicely in // pgplotters on mcdisplay. // One could get the current global rotation and use this to counteract this effect. double height = Geometries[index]->geometry_parameters.p_cylinder_storage->height; double radius = Geometries[index]->geometry_parameters.p_cylinder_storage->cyl_radius; Coords direction = Geometries[index]->geometry_parameters.p_cylinder_storage->direction_vector; Coords center = Geometries[index]->center; Coords bottom_point = coords_add (center, coords_scalar_mult (direction, 0.5 * height)); Coords top_point = coords_add (center, coords_scalar_mult (direction, -0.5 * height)); struct lines_to_draw lines_to_draw_temp; lines_to_draw_temp.number_of_lines = 0; lines_to_draw_temp = draw_circle_with_highest_priority (top_point, direction, radius, index, Geometries, number_of_volumes, 2); merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp); lines_to_draw_temp = draw_circle_with_highest_priority (bottom_point, direction, radius, index, Geometries, number_of_volumes, 2); merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp); Coords point1, point2; int iterate, number_of_points = 4; for (iterate = 0; iterate < number_of_points; iterate++) { point1 = point_on_circle (top_point, direction, radius, iterate, number_of_points); point2 = point_on_circle (bottom_point, direction, radius, iterate, number_of_points); lines_to_draw_temp = draw_line_with_highest_priority (point1, point2, index, Geometries, number_of_volumes, 2); merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp); } }; void initialize_cylinder_geometry_from_main_component (struct geometry_struct* cylinder) { // Function to be called in initialize of the main component // This is done as the rotation matrix needs to be relative to the main component instead of global // Everything done in initialize in this component file has the rotation matrix relative to global Coords simple_vector; Coords cyl_vector; // Start with vector that points along the cylinder in the local frame simple_vector = coords_set (0, 1, 0); // Rotate the direction vector of the cylinder to the master component frame of reference cyl_vector = rot_apply (cylinder->rotation_matrix, simple_vector); NORM (cyl_vector.x, cyl_vector.y, cyl_vector.z); cylinder->geometry_parameters.p_cylinder_storage->direction_vector.x = cyl_vector.x; cylinder->geometry_parameters.p_cylinder_storage->direction_vector.y = cyl_vector.y; cylinder->geometry_parameters.p_cylinder_storage->direction_vector.z = cyl_vector.z; // if (verbal == 1) printf("Cords vector1 = (%f,%f,%f)\n",cyl_vector.x,cyl_vector.y, } struct pointer_to_1d_coords_list cylinder_shell_points (struct geometry_struct* geometry, int max_number_of_points) { // Function that returns a number (less than max) of points on the geometry surface // If used, remember to free the space allocated. int points_per_circle = floor (max_number_of_points / 2.0); struct pointer_to_1d_coords_list cylinder_shell_array; cylinder_shell_array.elements = malloc (2 * points_per_circle * sizeof (Coords)); cylinder_shell_array.num_elements = 2 * points_per_circle; Coords cyl_direction = geometry->geometry_parameters.p_cylinder_storage->direction_vector; Coords center = geometry->center; double radius = geometry->geometry_parameters.p_cylinder_storage->cyl_radius; double height = geometry->geometry_parameters.p_cylinder_storage->height; Coords cyl_top_point = coords_add (center, coords_scalar_mult (cyl_direction, 0.5 * height)); Coords cyl_bottom_point = coords_add (center, coords_scalar_mult (cyl_direction, -0.5 * height)); points_on_circle (cylinder_shell_array.elements, cyl_top_point, cyl_direction, radius, points_per_circle); // Need to verify this pointer arithimatic works as intended points_on_circle (cylinder_shell_array.elements + points_per_circle, cyl_bottom_point, cyl_direction, radius, points_per_circle); return cylinder_shell_array; } #ifndef ANY_GEOMETRY_DETECTOR_DECLARE #define ANY_GEOMETRY_DETECTOR_DECLARE dummy // struct pointer_to_global_geometry_list global_geometry_list = {0,NULL}; #endif %} DECLARE %{ // Needed for transport to the main component struct global_geometry_element_struct global_geometry_element; int loop_index; int loop_2_index; int material_index; struct Volume_struct this_cylinder_volume; struct cylinder_storage this_cylinder_storage; struct surface_stack_struct curved_surface_stack; struct surface_stack_struct top_surface_stack; struct surface_stack_struct bottom_surface_stack; struct surface_stack_struct cut_surface_stack; %} INITIALIZE %{ struct union_state_struct* union_state_p = union_acquire (); geometry_struct_init (&(this_cylinder_volume.geometry)); // Initializes the focusing system for this volume including input sanitation. focus_initialize (&this_cylinder_volume.geometry, POS_A_COMP_INDEX (INDEX_CURRENT_COMP + target_index), POS_A_CURRENT_COMP, ROT_A_CURRENT_COMP, target_index, target_x, target_y, target_z, focus_aw, focus_ah, focus_xw, focus_xh, focus_r, NAME_CURRENT_COMP); // Input sanitation for this geometry if (radius <= 0) { printf ("\nERROR in Union_cylinder named %s, the radius is <= 0. \n", NAME_CURRENT_COMP); exit (1); } if (yheight <= 0) { printf ("\nERROR in Union_cylinder named %s, yheight is <= 0. \n", NAME_CURRENT_COMP); exit (1); } struct pointer_to_global_material_list* global_material_list = &union_state_p->u_material_list; // Use sanitation #ifdef MATERIAL_DETECTOR if (global_material_list->num_elements == 0) { // Here if the user have defined a material, but only after this material printf ("\nERROR: Need to define a material using Union_make_material before using a Union geometry component. \n"); printf (" %s was defined before first use of Union_make_material.\n", NAME_CURRENT_COMP); exit (1); } #endif #ifndef MATERIAL_DETECTOR printf ("\nERROR: Need to define a material using Union_make_material before using a Union geometry component. \n"); exit (1); #endif this_cylinder_volume.geometry.is_masked_volume = 0; this_cylinder_volume.geometry.is_exit_volume = 0; this_cylinder_volume.geometry.is_mask_volume = 0; struct pointer_to_global_geometry_list* global_geometry_list = &union_state_p->u_geometry_list; // Read the material input, or if it lacks, use automatic linking. if (mask_string && strlen (mask_string) && strcmp (mask_string, "NULL") && strcmp (mask_string, "0")) { // A mask volume is used to limit the extend of other volumes, called the masked volumes. These are specified in the mask_string. // In order for a ray to enter a masked volume, it needs to be both in the region covered by that volume AND the mask volume. // When more than this_cylinder_volume.geometry.mask_mode = 1; // Default is mask mode is ALL if (mask_setting && strlen (mask_setting) && strcmp (mask_setting, "NULL") && strcmp (mask_setting, "0")) { if (strcmp (mask_setting, "ALL") == 0 || strcmp (mask_setting, "All") == 0) this_cylinder_volume.geometry.mask_mode = 1; else if (strcmp (mask_setting, "ANY") == 0 || strcmp (mask_setting, "Any") == 0) this_cylinder_volume.geometry.mask_mode = 2; else { printf ("The mask_mode of component %s is set to %s, but must be either ALL or ANY.\n", NAME_CURRENT_COMP, mask_setting); exit (1); } } int found_geometries = 0; for (loop_index = 0; loop_index < global_geometry_list->num_elements; loop_index++) { // Add mask list if (1 == manual_linking_function (global_geometry_list->elements[loop_index].name, mask_string)) { add_element_to_int_list (&this_cylinder_volume.geometry.mask_list, global_geometry_list->elements[loop_index].component_index); add_element_to_int_list (&global_geometry_list->elements[loop_index].Volume->geometry.masked_by_list, INDEX_CURRENT_COMP); global_geometry_list->elements[loop_index].Volume->geometry.is_masked_volume = 1; if (this_cylinder_volume.geometry.mask_mode == 2) global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 2; if (this_cylinder_volume.geometry.mask_mode == 1) { if (global_geometry_list->elements[loop_index].Volume->geometry.is_masked_volume == 1 && global_geometry_list->elements[loop_index].Volume->geometry.mask_mode != 2) // If more than one mask is added to one volume, the ANY mode overwrites the (default) ALL mode. global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 1; } found_geometries = 1; } } if (found_geometries == 0) { printf ("The mask_string in geometry: %s did not find any of the specified volumes in the mask_string %s \n", NAME_CURRENT_COMP, mask_string); exit (1); } this_cylinder_volume.p_physics = malloc (sizeof (struct physics_struct)); this_cylinder_volume.p_physics->is_vacuum = 0; // Makes this volume a vacuum this_cylinder_volume.p_physics->number_of_processes = (int)0; // Should not be used. this_cylinder_volume.p_physics->my_a = 0; // Should not be used. sprintf (this_cylinder_volume.p_physics->name, "Mask"); this_cylinder_volume.geometry.is_mask_volume = 1; // Read the material input, or if it lacks, use automatic linking. } else if (material_string && strlen (material_string) && strcmp (material_string, "NULL") && strcmp (material_string, "0")) { // A geometry string was given, use it to determine which material if (0 == strcmp (material_string, "vacuum") || 0 == strcmp (material_string, "Vacuum")) { // One could have a global physics struct for vacuum instead of creating one for each this_cylinder_volume.p_physics = malloc (sizeof (struct physics_struct)); this_cylinder_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum this_cylinder_volume.p_physics->number_of_processes = (int)0; this_cylinder_volume.p_physics->my_a = 0; // Should not be used. sprintf (this_cylinder_volume.p_physics->name, "Vacuum"); } else if (0 == strcmp (material_string, "exit") || 0 == strcmp (material_string, "Exit")) { // One could have a global physics struct for exit instead of creating one for each this_cylinder_volume.p_physics = malloc (sizeof (struct physics_struct)); this_cylinder_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum this_cylinder_volume.p_physics->number_of_processes = (int)0; this_cylinder_volume.p_physics->my_a = 0; // Should not be used. this_cylinder_volume.geometry.is_exit_volume = 1; sprintf (this_cylinder_volume.p_physics->name, "Exit"); } else { for (loop_index = 0; loop_index < global_material_list->num_elements; loop_index++) { if (0 == strcmp (material_string, global_material_list->elements[loop_index].name)) { this_cylinder_volume.p_physics = global_material_list->elements[loop_index].physics; break; } if (loop_index == global_material_list->num_elements - 1) { printf ("\n"); printf ("ERROR: The material string \"%s\" in Union geometry \"%s\" did not match a specified material. \n", material_string, NAME_CURRENT_COMP); printf (" The materials available at this point (need to be defined before the geometry): \n"); for (loop_index = 0; loop_index < global_material_list->num_elements; loop_index++) printf (" %s\n", global_material_list->elements[loop_index].name); printf ("\n"); printf (" It is also possible to use one of the defualt materials avaiable: \n"); printf (" Vacuum (for a Volume without scattering or absorption)\n"); printf (" Exit (for a Volume where the ray exits the component if it enters)\n"); printf (" Mask (for a Volume that masks existing volumes specified in the mask_string\n"); exit (1); } } } } else { // Automatic linking, simply using the last defined material. #ifndef MATERIAL_DETECTOR printf ("Need to define a material before the geometry to use automatic linking %s.\n", NAME_CURRENT_COMP); exit (1); #endif this_cylinder_volume.p_physics = global_material_list->elements[global_material_list->num_elements - 1].physics; } // Handle surface input if (all_face_surface && strlen (all_face_surface) && strcmp (all_face_surface, "NULL") && strcmp (all_face_surface, "0")) { // Overwrite all other surfaces that have no input (if one wants a surface left without anything, None can be used) overwrite_if_empty (curved_surface, all_face_surface); overwrite_if_empty (top_surface, all_face_surface); overwrite_if_empty (bottom_surface, all_face_surface); } this_cylinder_volume.geometry.number_of_faces = 3; this_cylinder_volume.geometry.surface_stack_for_each_face = malloc (this_cylinder_volume.geometry.number_of_faces * sizeof (struct surface_stack_struct*)); // This could be inserted into the fill_surface_stack function this_cylinder_volume.geometry.surface_stack_for_each_face[0] = &curved_surface_stack; this_cylinder_volume.geometry.surface_stack_for_each_face[1] = &top_surface_stack; this_cylinder_volume.geometry.surface_stack_for_each_face[2] = &bottom_surface_stack; this_cylinder_volume.geometry.internal_cut_surface_stack = &cut_surface_stack; // This surface is used if the volume is cut by overlapping of higher priority volume struct pointer_to_global_surface_list* global_surface_list = &union_state_p->u_surface_list; fill_surface_stack (curved_surface, global_surface_list, NAME_CURRENT_COMP, &curved_surface_stack); fill_surface_stack (top_surface, global_surface_list, NAME_CURRENT_COMP, &top_surface_stack); fill_surface_stack (bottom_surface, global_surface_list, NAME_CURRENT_COMP, &bottom_surface_stack); fill_surface_stack (cut_surface, global_surface_list, NAME_CURRENT_COMP, &cut_surface_stack); sprintf (this_cylinder_volume.name, "%s", NAME_CURRENT_COMP); sprintf (this_cylinder_volume.geometry.shape, "cylinder"); this_cylinder_volume.geometry.eShape = cylinder; this_cylinder_volume.geometry.priority_value = priority; // Currently the coordinates will be in absolute space. this_cylinder_volume.geometry.center = POS_A_CURRENT_COMP; this_cylinder_volume.geometry.geometry_p_interact = p_interact; this_cylinder_storage.cyl_radius = radius; this_cylinder_storage.height = yheight; this_cylinder_volume.geometry.visualization_on = visualize; this_cylinder_volume.geometry.geometry_parameters.p_cylinder_storage = &this_cylinder_storage; this_cylinder_volume.geometry.within_function = &r_within_cylinder; this_cylinder_volume.geometry.intersect_function = &sample_cylinder_intersect; this_cylinder_volume.geometry.mcdisplay_function = &mcdisplay_cylinder_function; this_cylinder_volume.geometry.shell_points = &cylinder_shell_points; this_cylinder_volume.geometry.initialize_from_main_function = &initialize_cylinder_geometry_from_main_component; this_cylinder_volume.geometry.process_rot_allocated = 0; this_cylinder_volume.geometry.copy_geometry_parameters = &allocate_cylinder_storage_copy; rot_copy (this_cylinder_volume.geometry.rotation_matrix, ROT_A_CURRENT_COMP); rot_transpose (ROT_A_CURRENT_COMP, this_cylinder_volume.geometry.transpose_rotation_matrix); // Initialize loggers this_cylinder_volume.loggers.num_elements = 0; this_cylinder_volume.abs_loggers.num_elements = 0; // packing the information into the global_geometry_element, which is then included in the global_geometry_list. sprintf (global_geometry_element.name, "%s", NAME_CURRENT_COMP); global_geometry_element.activation_counter = number_of_activations; global_geometry_element.component_index = INDEX_CURRENT_COMP; global_geometry_element.Volume = &this_cylinder_volume; // Would be nicer if this m was a pointer, now we have the (small) data two places add_element_to_geometry_list (global_geometry_list, global_geometry_element); %} TRACE %{ %} FINALLY %{ union_release (NAME_CURRENT_COMP); %} END