/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2003, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Sans_spheres * * %I * Written by: P. Willendrup, K. Lefmann, L. Arleth * Date: 19.12.2003 * Origin: Risoe * Modified by: KL, 7 June 2005 * * Sample for Small Angle Neutron Scattering - hard spheres in thin solution, mono disperse. * * %D * Sample for use in a SANS instrument, models hard, mono disperse spheres in thin solution. * The shape of the sample may be a filled box with dimensions * xwidth, yheight, zdepth, a cylinder with dimensions radius and yheight, * a filled sphere with radius. * * Example: Sans_spheres(R = 100, Phi = 1e-3, Delta_rho = 0.6, sigma_abs = 50, xwidth=0.01, yheight=0.01, zdepth=0.005) * * %P * * INPUT PARAMETERS * * R: [AA] Radius of scattering hard spheres * Phi: [1] Particle volume fraction * Delta_rho: [fm/AA^3] Excess scattering length density * sigma_abs: [m^-1] Absorption cross section density at 2200 m/s * radius: [m] Outer radius of sample in (x,z) plane for cylinder/sphere * xwidth: [m] horiz. dimension of sample, as a width * yheight: [m] vert . dimension of sample, as a height for cylinder/box * zdepth: [m] depth of sample * target_index: [1] Relative index of component to focus at, e.g. next is +1 * focus_xw: [m] horiz. dimension of a rectangular area * focus_yh: [m] vert. dimension of a rectangular area * focus_aw: [deg] horiz. angular dimension of a rectangular area * focus_ah: [deg] vert. angular dimension of a rectangular area * focus_r: [m] Detector (disk-shaped) radius * * Optional parameters: * target_x: [m] * target_y: [m] position of target to focus at * target_z: [m] * * Variables calculated in the component * * my_s: Attenuation factor due to scattering [m^-1] * my_a: Attenuation factor due to absorbtion [m^-1] * * %Link * The test/example instrument SANS.instr. %L * Some alternative implementations exist as contributed components. * %E *******************************************************************************/ DEFINE COMPONENT Sans_spheres DEFINITION PARAMETERS () SETTING PARAMETERS (R=100, Phi=1e-3, Delta_rho=0.6, sigma_abs=0.05, xwidth=0, yheight=0, zdepth=0, radius=0, target_x = 0, target_y = 0, target_z = 6, int target_index=0, focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, focus_r=0) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ DECLARE %{ double my_s_pre; double my_a_v; double shape; %} INITIALIZE %{ shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere */ if (xwidth && yheight && zdepth) shape=1; /* box */ else if (radius > 0 && yheight) shape=0; /* cylinder */ else if (radius > 0 && !yheight) shape=2; /* sphere */ if (shape < 0) exit(fprintf(stderr,"Sans_spheres: %s: sample has invalid dimensions.\n" "ERROR Please check parameter values.\n", NAME_CURRENT_COMP)); /* now compute target coords if a component index is supplied */ if (!target_index && !target_x && !target_y && !target_z) target_index=1; if (target_index) { Coords ToTarget; ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP); ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget); coords_get(ToTarget, &target_x, &target_y, &target_z); } if (!(target_x || target_y || target_z)) { printf("Sans_spheres: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP); target_z=1; } my_a_v = sigma_abs*2200*100; /* Is not yet divided by v. 100: Convert barns -> fm^2 */ my_s_pre = Phi * 4*PI*R*R*R/3 * Delta_rho*Delta_rho; %} TRACE %{ double t0, t1, v, l_full, l, l_1, dt, d_phi, theta, my_s; double aim_x=0, aim_y=0, aim_z=1, axis_x, axis_y, axis_z; double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z; double f, solid_angle, vx_i, vy_i, vz_i, qx, qy, qz,q; char intersect=0; /* Intersection neutron trajectory / sample (sample surface) */ if (shape == 0) intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight); else if (shape == 1) intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); else if (shape == 2) intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius); if(intersect) { if(t0 < 0) ABSORB; /* Neutron enters at t=t0. */ v = sqrt(vx*vx + vy*vy + vz*vz); l_full = v * (t1 - t0); /* Length of full path through sample */ dt = rand01()*(t1 - t0) + t0; /* Time of scattering */ PROP_DT(dt); /* Point of scattering */ l = v*(dt-t0); /* Penetration in sample */ vx_i=vx; vy_i=vy; vz_i=vz; if ((target_x || target_y || target_z)) { aim_x = target_x-x; /* Vector pointing at target (anal./det.) */ aim_y = target_y-y; aim_z = target_z-z; } if(focus_aw && focus_ah) { randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP); } else if(focus_xw && focus_yh) { randvec_target_rect(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); } else { randvec_target_circle(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r); } NORM(vx, vy, vz); vx *= v; vy *= v; vz *= v; qx = V2K*(vx_i-vx); qy = V2K*(vy_i-vy); qz = V2K*(vz_i-vz); q = sqrt(qx*qx+qy*qy+qz*qz); f = 3 * (sin(q*R) - q*R*cos(q*R))/(q*R*q*R*q*R); l_1 = v*t1; double pmul=l_full*solid_angle/(4*PI)*my_s_pre*f*f*exp(-my_a_v*(l+l_1)/v); p = p * pmul; SCATTER; } %} MCDISPLAY %{ if (shape == 0) { /* cylinder */ circle("xz", 0, yheight/2.0, 0, radius); circle("xz", 0, -yheight/2.0, 0, radius); line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0); line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0); line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius); line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius); } else if (shape == 1) { /* box */ double xmin = -0.5*xwidth; double xmax = 0.5*xwidth; double ymin = -0.5*yheight; double ymax = 0.5*yheight; double zmin = -0.5*zdepth; double zmax = 0.5*zdepth; multiline(5, xmin, ymin, zmin, xmax, ymin, zmin, xmax, ymax, zmin, xmin, ymax, zmin, xmin, ymin, zmin); multiline(5, xmin, ymin, zmax, xmax, ymin, zmax, xmax, ymax, zmax, xmin, ymax, zmax, xmin, ymin, zmax); line(xmin, ymin, zmin, xmin, ymin, zmax); line(xmax, ymin, zmin, xmax, ymin, zmax); line(xmin, ymax, zmin, xmin, ymax, zmax); line(xmax, ymax, zmin, xmax, ymax, zmax); } else if (shape == 2) { /* sphere */ circle("xy", 0, 0.0, 0, radius); circle("xz", 0, 0.0, 0, radius); circle("yz", 0, 0.0, 0, radius); } %} END