Fixed bugs in union-lib.c relating to calculation of normal vectors for box and cylinder.

Added possibility to specify scattering length density directly in Union_make_material, and simplified memory structure to primarily store SLD instead of more variables used to calculate it.

Improved debug prints in Union_master and adopted use of SLD.

Author: mads-bertelsen

mcstas-comps/share/union-lib.c

@@ -471,8 +471,7 @@ struct scattering_process_struct *p_scattering_array;
 
 // refraction related
 int has_refraction_info;
-double refraction_rho;
-double refraction_bc;
+double refraction_scattering_length_density; // [AA^-2]
 double refraction_Qc;
 };
 
@@ -2646,11 +2645,7 @@ int sample_box_intersect_advanced(double *t, double *nx, double *ny, double *nz,
 	Coords normal_vector_rotated;
 	Coords normal_vector;
 
-	
-	// Rotate back to master coordinate system	
-	normal_vector_rotated = coords_set(nx[0], nx[1], nx[2]);
-    normal_vector = rot_apply(geometry->rotation_matrix, normal_vector_rotated);
-    
+	// Sort solution according to intersection time and rotate normal vectors to master coordinate system 
     switch(*num_solutions) {
     case 2:
         if (t[0] > t[1]) {
@@ -2668,7 +2663,7 @@ int sample_box_intersect_advanced(double *t, double *nx, double *ny, double *nz,
 			ny[0] = temp;
 
 			temp = nz[1];
-			nz[1] = nx[0];
+			nz[1] = nz[0];
 			nz[0] = temp;
 
 			// Switch surface_index
@@ -4127,7 +4122,7 @@ int sample_cylinder_intersect(double *t, double *nx, double *ny, double *nz, int
     double radius = geometry->geometry_parameters.p_cylinder_storage->cyl_radius;
     double height = geometry->geometry_parameters.p_cylinder_storage->height;
 
-    // Declare variables for the function
+    // Declare position variables for the local coordinate system
     double x_new,y_new,z_new;
 
     // Coordinate transformation
@@ -4137,69 +4132,77 @@ int sample_cylinder_intersect(double *t, double *nx, double *ny, double *nz, int
 
     Coords coordinates = coords_set(x_new,y_new,z_new);
     Coords rotated_coordinates;
-    // printf("Cords coordinates = (%f,%f,%f)\n",coordinates.x,coordinates.y,coordinates.z);
-    
-    // debug
-    // Rotation rotation_matrix_debug[3][3];
-    // rot_set_rotation(rotation_matrix_debug,-1.0*geometry->rotation.x,-1.0*geometry->rotation.y,-1.0*geometry->rotation.z);
-    // rot_transpose(geometry->rotation_matrix,rotation_matrix_debug);
 
     // Rotate the position of the neutron around the center of the cylinder
     rotated_coordinates = rot_apply(geometry->transpose_rotation_matrix,coordinates);
-    // rotated_coordinates = rot_apply(rotation_matrix_debug,coordinates);
-    //     printf("Cords rotated_coordinates = (%f,%f,%f)\n",rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z);
 
     Coords velocity = coords_set(v[0],v[1],v[2]);
     Coords rotated_velocity;
-    //     printf("Cords velocity = (%f,%f,%f)\n",velocity.x,velocity.y,velocity.z);
 
     // Rotate the position of the neutron around the center of the cylinder
     rotated_velocity = rot_apply(geometry->transpose_rotation_matrix,velocity);
-    // rotated_velocity = rot_apply(rotation_matrix_debug,velocity);
-    //     printf("Cords rotated_velocity = (%f,%f,%f)\n",rotated_velocity.x,rotated_velocity.y,rotated_velocity.z);
-    
-    
-    
-    // Cases where the velocity is parallel with the cylinder axis have given problems, and is checked for explicitly
-    if (sqrt(rotated_velocity.x*rotated_velocity.x+rotated_velocity.z*rotated_velocity.z)/fabs(rotated_velocity.y) < 0.00001) {
-      // The velocity is parallel with the cylinder axis. Either there is two solutions
-      if (sqrt(rotated_coordinates.x*rotated_coordinates.x+rotated_coordinates.z*rotated_coordinates.z) > radius) {
-        *num_solutions = 0;
-        return 0;
-      } else {
-        *num_solutions = 2;
-        t[0] = (0.5*height - rotated_coordinates.y)/rotated_velocity.y;
-        t[1] = (-0.5*height - rotated_coordinates.y)/rotated_velocity.y;
-        return 1;
-      }
-    }
 
-    int output;
-    // Run McStas built in sphere intersect funtion (sphere centered around origin)
-    if ((output = cylinder_intersect(&t[0],&t[1],
-	                                 rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,
-	                                 rotated_velocity.x,rotated_velocity.y,rotated_velocity.z,radius,height)) == 0) {
-        *num_solutions = 0;t[0]=-1;t[1]=-1;
-	}
-    else if (t[1] != 0) *num_solutions = 2;
-    else {*num_solutions = 1; t[1]=-1;}
+	int output;    
+	// Cases where the velocity is parallel with the cylinder axis have given problems, and is checked for explicitly
+	if (sqrt(rotated_velocity.x*rotated_velocity.x+rotated_velocity.z*rotated_velocity.z)/fabs(rotated_velocity.y) < 0.00001) {
+	  // The velocity is parallel with the cylinder axis. Either there are no solutions or two solutions
+	  if (sqrt(rotated_coordinates.x*rotated_coordinates.x+rotated_coordinates.z*rotated_coordinates.z) > radius) {
+	    *num_solutions = 0;
+	    return 0;
+	  } else {
+	    *num_solutions = 2;
+	    t[0] = (0.5*height - rotated_coordinates.y)/rotated_velocity.y;
+		surface_index[0] = 1; // index indicating top
+		
+	    t[1] = (-0.5*height - rotated_coordinates.y)/rotated_velocity.y;
+		surface_index[1] = 2; // index indicating bottom
+		
+		// sort solutions
+		if (t[0] > t[1]) {
+		  double d_temp;
+		  
+		  d_temp = t[0];
+		  t[0] = t[1];
+		  t[1] = d_temp;
+		  
+		  int i_temp;
+		  
+		  i_temp = surface_index[0];
+		  surface_index[0] = surface_index[1];
+		  surface_index[1] = i_temp;
+		}
+	  }
+	} else {
+	  // velocity not parallel to cylinder axis, call standard mcstas cylinder intersect 
+
+	  // Run McStas built in sphere intersect funtion (sphere centered around origin)
+	  if ((output = cylinder_intersect(&t[0],&t[1],
+	                                   rotated_coordinates.x,rotated_coordinates.y,rotated_coordinates.z,
+	                                   rotated_velocity.x,rotated_velocity.y,rotated_velocity.z,radius,height)) == 0) {
+	      *num_solutions = 0;t[0]=-1;t[1]=-1;
+	  }
+	  else if (t[1] != 0) *num_solutions = 2;
+	  else {*num_solutions = 1; t[1]=-1;}
 
-	// decode output value
-    // Check the bitmask for entry and exit
-	if (*num_solutions > 0) {
-	  int entry_index = 0;
-      if (output & 2) entry_index = 1; // Entry intersects top cap
-      if (output & 4) entry_index = 2; // Entry intersects bottom cap
-	  surface_index[0] = entry_index;
-    }
+	  // decode output value
+	  // Check the bitmask for entry and exit
+	  if (*num_solutions > 0) {
+	    int entry_index = 0;
+	    if (output & 2) entry_index = 1; // Entry intersects top cap
+	    if (output & 4) entry_index = 2; // Entry intersects bottom cap
+	    surface_index[0] = entry_index;
+	  }
 
-	if (*num_solutions > 1) {
-	  int exit_index = 0;		
-      if (output & 8) exit_index = 1;  // Exit intersects top cap
-      if (output & 16) exit_index = 2; // Exit intersects bottom cap
-	  surface_index[1] = exit_index;
+	  if (*num_solutions > 1) {
+	    int exit_index = 0;		
+	    if (output & 8) exit_index = 1;  // Exit intersects top cap
+	    if (output & 16) exit_index = 2; // Exit intersects bottom cap
+	    surface_index[1] = exit_index;
+	  }
 	}
+    
 
+	// Calculate normal vectors from surface index and cylinder geometry
     int index;
     double x, y, z, dt;
     Coords normal_vector_rotated;
@@ -4238,7 +4241,7 @@ int sample_cylinder_intersect(double *t, double *nx, double *ny, double *nz, int
 };
 
 int r_within_cylinder(Coords pos,struct geometry_struct *geometry) {
-// Unpack parameters
+	// Unpack parameters
     double radius = geometry->geometry_parameters.p_cylinder_storage->cyl_radius;
     double height = geometry->geometry_parameters.p_cylinder_storage->height;
 

mcstas-comps/union/Single_crystal_process.comp

@@ -54,10 +54,8 @@
 * cy: [AA or AA^-1]                                       c on y axis
 * cz: [AA or AA^-1]                                       c on z axis
 * reflections: [string]                                   File name containing structure factors of reflections. Use empty ("") or NULL for incoherent scattering only
-* order: [1]                                              Limit multiple scattering up to given order (0: all, 1: first, 2: second, ...)
+* order: [1]                                              Limit multiple scattering up to given order (0: all, 1: first, 2: second, ...) (Not supported in Union)
 * p_transmit: [1]                                         Monte Carlo probability for neutrons to be transmitted without any scattering. Used to improve statistics from weak reflections
-* sigma_abs: [barns]                                      Absorption cross-section per unit cell at 2200 m/s
-* sigma_inc: [barns]                                      Incoherent scattering cross-section per unit cell Use -1 to inactivate
 * aa: [deg]                                               Unit cell angles alpha, beta and gamma. Then uses norms of vectors a,b and c as lattice parameters
 * bb: [deg]                                               Beta angle
 * cc: [deg]                                               Gamma angle

mcstas-comps/union/Union_make_material.comp

@@ -37,6 +37,7 @@
 * refraction_sigma_coh: [barn]     coherent cross section of refracting material. Use negative value to indicate a negative coherent scattering length
 * refraction_density: [g/cm3]      density of the refracting material. density < 0 means the material is outside/before the shape.
 * refraction_weight: [g/mol]       molar mass of the refracting material
+* refraction_SLD: [AA^-2]       scattering length density of material (overwrites sigma_coh, density and weight)
 * init:           [string] Name of Union_init component (typically "init", default)
 *
 * CALCULATED PARAMETERS:
@@ -53,7 +54,10 @@
 
 DEFINE COMPONENT Union_make_material
 
-SETTING PARAMETERS(string process_string="NULL", my_absorption,absorber=0, refraction_density=0, refraction_sigma_coh=0, refraction_weight=0, string init="init")
+SETTING PARAMETERS(string process_string="NULL", my_absorption,absorber=0, 
+refraction_density=0, refraction_sigma_coh=0, refraction_weight=0, 
+refraction_SLD=-1500,
+string init="init")
 
 
 /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
@@ -237,22 +241,32 @@ struct pointer_to_global_process_list *global_process_list = COMP_GETPAR3(Union_
 
   this_material.has_refraction_info = 0;
   if (refraction_density != 0 || refraction_weight != 0 || refraction_sigma_coh != 0) {
+	  
+	  double refraction_rho, refraction_bc;
       if (refraction_density==0 || refraction_weight<=0)
         exit(printf("Union_make_material: %s: FATAL: invalid material density or molar weight: density=%g weight=%g\n",
           NAME_CURRENT_COMP, refraction_density, refraction_weight));
-      this_material.refraction_rho = fabs(refraction_density)*6.02214179*1e23*1e-24/refraction_weight; // per at/Angs^3
+
+	  refraction_rho = fabs(refraction_density)*6.02214179*1e23*1e-24/refraction_weight; // per at/Angs^3  
 
       if (refraction_sigma_coh==0)
         exit(printf("Refractor: %s: FATAL: invalid material coherent cross section: sigma_coh=%g\n",
           NAME_CURRENT_COMP, refraction_sigma_coh));
-      this_material.refraction_bc = sqrt(fabs(refraction_sigma_coh)*100/4/PI)*1e-5;         // bound coherent scattering length
-      if (refraction_sigma_coh<0) this_material.refraction_bc *= -1;
-    
-      this_material.refraction_Qc = 4*sqrt(PI*this_material.refraction_rho*fabs(this_material.refraction_bc));
-    
+
+	  refraction_bc = sqrt(fabs(refraction_sigma_coh)*100/4/PI)*1e-5;         // bound coherent scattering length
+	  if (refraction_sigma_coh<0) refraction_bc *= -1.0;
+	  
+	  this_material.refraction_scattering_length_density = refraction_rho*refraction_bc;
+	  this_material.refraction_Qc = 4*sqrt(PI*refraction_rho*fabs(refraction_bc));
       this_material.has_refraction_info = 1;
   }
 
+  if (refraction_SLD > -1499.0) { // refraction_SLD can be negative, zero or positive, it can however not be this negative, so this is used as input check
+	  this_material.refraction_scattering_length_density = refraction_SLD;  	
+	  this_material.refraction_Qc = 4.0*sqrt(PI*fabs(refraction_SLD));  
+	  this_material.has_refraction_info = 1;
+  }
+  
   // packing the information into the global_material_element, which is then included in the global_material_list.
   sprintf(global_material_element.name,"%s",NAME_CURRENT_COMP);
   global_material_element.component_index = INDEX_CURRENT_COMP;