/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2006, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Exact_radial_coll * * %I * Written by: Roland Schedler * Modified by: using Collimator_radial Component * Date: October 2006 * Origin: HMI * Modified by: E. Farhi, uniformize parameter names (Jul 2008) * * An exact radial Soller collimator. * * %D * Radial Soller collimator with rectangular opening, specified length and * specified foil thickness. * The collimator is made of many trapezium shaped nslit stacked radially. * The nslit are separated by absorbing foils, the whole stuff is inside * an absorbing housing. * The component should be positioned at the radius center. The model is exact. * The neutron beam outside the collimator area is transmitted unaffected. * * Example: Exact_radial_coll(theta_min=-5, theta_max=5, nslit=100, * radius=1.0, length=.3, h_in=.2, h_out=.3, d=0.0001) * * * %P * INPUT PARAMETERS: * * theta_min: [deg] Minimum Theta angle for the radial setting * theta_max: [deg] Maximum Theta angle for the radial setting * nslit: [1] Number of channels in the theta range * radius: [m] Inner radius (focus point to foil start point). * length: [m] Length of the foils / collimator * h_in: [m] Input window height * h_out: [m] Output window height * d: [m] Thickness of the absorbing foils * verbose: [0/1] Gives additional information * * %E *******************************************************************************/ DEFINE COMPONENT Exact_radial_coll DEFINITION PARAMETERS () SETTING PARAMETERS (theta_min=-5, theta_max=5, nslit=100, radius=1.0, length=.5, h_in=.3, h_out=.4, d=0.0001, verbose=0) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ DECLARE %{ double alpha_in; double alpha_out; double beta_in; double beta_out; double theta; double out_radius; double iw; double ow; double divergence; %} INITIALIZE %{ /* check for input parameters */ if (radius <= 0) exit(printf("Exact_radial_coll: %s: radius must be positive\n", NAME_CURRENT_COMP)); if (h_in <= 0) exit(printf("Exact_radial_coll: %s: h_in must be positive\n", NAME_CURRENT_COMP)); if (h_out <= 0) exit(printf("Exact_radial_coll: %s: h_out must be positive\n", NAME_CURRENT_COMP)); if (d <= 0) exit(printf("Exact_radial_coll: %s: d must be positive\n", NAME_CURRENT_COMP)); if (nslit <= 0) exit(printf("Exact_radial_coll: %s: number of channels must be positive\n", NAME_CURRENT_COMP)); if ((nslit - floor (nslit)) > 0) exit(printf("Exact_radial_coll: %s: number of channels must be an integer\n", NAME_CURRENT_COMP)); if (length <= 0) exit(printf("Exact_radial_coll: %s: collimator length must be positive\n", NAME_CURRENT_COMP)); if (theta_max <= theta_min) exit(printf("Exact_radial_coll: %s: theta_max must be greater than theta_min\n", NAME_CURRENT_COMP)); theta_max *= DEG2RAD; theta_min *= DEG2RAD; theta = theta_max - theta_min; out_radius = radius + length; beta_in = 2*asin(d / (2 * radius)); beta_out= 2*asin(d / (2 * out_radius)); if (theta < nslit*beta_in) exit(printf("Exact_radial_coll: %s: the %6.0f foils of %g [meter]\n" "do not fit within the angular range theta = %4.2f [deg]\n", NAME_CURRENT_COMP, nslit, d, theta*RAD2DEG)); alpha_in = (theta - nslit*beta_in)/nslit; alpha_out = (theta - nslit*beta_out)/nslit; iw = 2*radius*sin((alpha_in/2)); ow = 2*out_radius*sin((alpha_out/2)); divergence=(iw+ow)/(sqrt(4*length*length-(ow-iw)*(ow-iw))); if (verbose) { printf("Exact_radial_coll: %s: foil thickness is %.2g [millimeter]\n", NAME_CURRENT_COMP, d*1000); printf(" opening each input slit [%.3g:%.0f] [millimeter]\n", iw*1000, h_in*1000); printf(" opening each output slit [%.3g:%.0f] [millimeter]\n", ow*1000, h_out*1000); printf(" divergence per channel is %g [min] \n", divergence*RAD2MIN); } %} TRACE %{ double phi, t0, t1, t2, t3; int intersect; long input_chan, output_chan; double input_theta, output_theta; double input_center,output_center; double window_theta; char ok=0; /* first compute intersection time with input cylinder */ intersect=cylinder_intersect(&t0,&t3,x,y,z,vx,vy,vz,radius,h_in); if (!intersect) ABSORB; else if (t3 > t0) t0 = t3; intersect=cylinder_intersect(&t1,&t2,x,y,z,vx,vy,vz,out_radius,h_out); if (!intersect) ABSORB; else if (t2 > t1) t1 = t2; /* get index of input slit */ if (t0 > 0 && t1 > t0) { PROP_DT(t0); input_theta = atan2(x, z); /* channel number (start at 0) */ window_theta = (theta_max - theta_min)/nslit; input_chan = floor((input_theta - theta_min)/window_theta); if (input_chan >= 0 && input_chan < nslit && fabs(y) < h_in/2) ok=1; if (ok) { input_center= theta_min + input_chan*window_theta + (window_theta)/2; /* are we outside the soller or in the foil? */ phi = input_theta - input_center; if (fabs(phi) > alpha_in/2) ABSORB; /* inside the foil*/ SCATTER; /* propagate to output radius */ PROP_DT(t1-t0); SCATTER; output_theta = atan2(x, z); /* channel number (start at 0) */ output_chan = floor((output_theta - theta_min)/window_theta); /* did we change channel ? */ if (output_chan != input_chan) ABSORB; /* changed slit */ output_center= theta_min + output_chan*window_theta + (window_theta)/2; /* are we outside the soller */ phi = output_theta -output_center; if (fabs(phi) > alpha_out/2 || fabs(y) > h_out/2) ABSORB; /* outside output slit */ } /* else we pass aside the entrance window of radial collimator */ else { /* propagate to output radius */ PROP_DT(t1-t0); SCATTER; output_theta = atan2(x, z); /* channel number (start at 0) */ output_chan = floor((output_theta - theta_min)/window_theta); /* are we come from outside into the soller or in the foil?*/ if (output_chan >= 0 || output_chan < nslit) ABSORB; } /* else we pass aside the exit window of radial collimator */ } /* else did not encounter collimator */ %} MCDISPLAY %{ int i; double theta1, theta2, theta3, theta4; double x_in_l, z_in_l, x_in_r, z_in_r; double x_out_l, z_out_l, x_out_r, z_out_r; double window_theta, y1, y2; window_theta = alpha_in + beta_in; y1 = h_in/2; y2 = h_out/2; theta1 = theta_min; theta3 = theta1+beta_in/2; theta4 = theta1+beta_out/2; z_in_l = radius*cos(theta1); x_in_l = radius*sin(theta1); z_in_r = radius*cos(theta3); x_in_r = radius*sin(theta3); z_out_l = out_radius*cos(theta1); x_out_l = out_radius*sin(theta1); z_out_r = out_radius*cos(theta4); x_out_r = out_radius*sin(theta4); multiline(5, x_in_l, -y1, z_in_l, x_in_l, y1, z_in_l, x_out_l, y2, z_out_l, x_out_l,-y2, z_out_l, x_in_l, -y1, z_in_l); line(x_in_l, y1, z_in_l, x_in_r, y1, z_in_r); line(x_in_l, -y1, z_in_l, x_in_r, -y1, z_in_r); line(x_out_l, y2, z_out_l, x_out_r, y2, z_out_r); line(x_out_l, -y2, z_out_l, x_out_r,-y2, z_out_r); multiline(5, x_in_r, -y1, z_in_r, x_in_r, y1, z_in_r, x_out_r, y2, z_out_r, x_out_r,-y2, z_out_r, x_in_r, -y1, z_in_r); for (i = 1; i < nslit; i++) { theta1 = i*window_theta+theta_min-beta_in/2; theta2 = i*window_theta+theta_min+beta_in/2; theta3 = i*window_theta+theta_min-beta_out/2; theta4 = i*window_theta+theta_min+beta_out/2; z_in_l = radius*cos(theta1); x_in_l = radius*sin(theta1); z_in_r = radius*cos(theta2); x_in_r = radius*sin(theta2); z_out_l = out_radius*cos(theta3); x_out_l = out_radius*sin(theta3); z_out_r = out_radius*cos(theta4); x_out_r = out_radius*sin(theta4); /* left side */ multiline(5, x_in_l, -y1, z_in_l, x_in_l, y1, z_in_l, x_out_l, y2, z_out_l, x_out_l,-y2, z_out_l, x_in_l, -y1, z_in_l); /* left -> right lines */ line(x_in_l, y1, z_in_l, x_in_r, y1, z_in_r); line(x_in_l, -y1, z_in_l, x_in_r, -y1, z_in_r); line(x_out_l, y2, z_out_l, x_out_r, y2, z_out_r); line(x_out_l, -y2, z_out_l, x_out_r,-y2, z_out_r); /* right side */ multiline(5, x_in_r, -y1, z_in_r, x_in_r, y1, z_in_r, x_out_r, y2, z_out_r, x_out_r,-y2, z_out_r, x_in_r, -y1, z_in_r); } /* remaining bits */ theta1 = theta_max; theta3 = theta1-beta_in/2; theta4 = theta1-beta_out/2; z_in_l = radius*cos(theta1); x_in_l = radius*sin(theta1); z_in_r = radius*cos(theta3); x_in_r = radius*sin(theta3); z_out_l = out_radius*cos(theta1); x_out_l = out_radius*sin(theta1); z_out_r = out_radius*cos(theta4); x_out_r = out_radius*sin(theta4); multiline(5, x_in_l, -y1, z_in_l, x_in_l, y1, z_in_l, x_out_l, y2, z_out_l, x_out_l,-y2, z_out_l, x_in_l, -y1, z_in_l); line(x_in_l, y1, z_in_l, x_in_r, y1, z_in_r); line(x_in_l, -y1, z_in_l, x_in_r, -y1, z_in_r); line(x_out_l, y2, z_out_l, x_out_r, y2, z_out_r); line(x_out_l, -y2, z_out_l, x_out_r,-y2, z_out_r); multiline(5, x_in_r, -y1, z_in_r, x_in_r, y1, z_in_r, x_out_r, y2, z_out_r, x_out_r,-y2, z_out_r, x_in_r, -y1, z_in_r); %} END