/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2003, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: SANS_Guinier * * %I * Written by: Henrich Frielinghaus * Date: Sept 2004 * Origin: FZ-Juelich/FRJ-2/IFF/KWS-2 * * Sample for Small Angle Neutron Scattering: Guinier model * * %D * Sample that scatters with a Guinier shape. This is just an example where analytically * an integral exists. The neutron paths are proportional to the intensity * (low intensity > few paths). * * Guinier function (Rg) * a = Rg*Rg/3 * propability_unscaled = q * exp(-a*q*q) * integral_prop_unscal = 1/(2*a) * (1 - exp(-a*q*q)) * propability_scaled = 2*a * q*exp(-a*q*q) / (1 - exp(-a*q*q)) * integral_prop_scaled = (1 - exp(-a*q*q)) / (1 - exp(-a*qmax*qmax)) * * In this simulation method many paths occur for high propability. * For simulation of low intensities see SANS_AnySamp. * * Sample components leave the units of flux for the probability * of the individual paths. That is more consitent than the * Sans_spheres routine. Furthermore one can simulate the * transmitted beam. This allows to determine the needed size of * the beam stop. Only absorption has not been included yet to * these sample-components. But thats really nothing. * * Example: SANS_Guinier(transm=0.5, Rg=100, qmax=0.03, xwidth=0.01, yheight=0.01, zdepth=0.001) * * %P * * INPUT PARAMETERS * * transm: [1] (coherent) transmission of sample for the optical path "zdepth" * Rg: [Angs] Radius of Gyration * qmax: [AA-1] Maximum scattering vector * xwidth: [m] horiz. dimension of sample, as a width * yheight: [m] vert.. dimension of sample, as a height * zdepth: [m] thickness of sample * * %Link * Sans_spheres component * * %E *******************************************************************************/ DEFINE COMPONENT SANS_Guinier DEFINITION PARAMETERS () SETTING PARAMETERS (transm=0.5, Rg=100, qmax=0.03, xwidth=0.01, yheight=0.01, zdepth=0.001) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ DECLARE %{ %} INITIALIZE %{ if (!xwidth || !yheight || !zdepth) { exit(fprintf(stderr,"SANS_Guinier: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP)); } %} TRACE %{ double a,qm,q,q_v; double transmr, t0, t1, v, l_full, l, dt, d_phi, theta; double axis_x, axis_y, axis_z; double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z; char intersect=0; transmr = transm; /* real transmission */ if (transmr<1e-10) transmr = 1e-10; if (transmr>1e0 ) transmr = 1e0; intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); 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 */ transmr = exp(log(transmr)*l_full/zdepth); /* real transmission */ dt = rand01()*(t1 - t0) + t0; /* Time of scattering */ PROP_DT(dt); /* Point of scattering */ l = v*dt; /* Penetration in sample */ a = Rg*Rg/3.0; qm= qmax; if (qm<1.0/Rg) qm = 1.0/Rg; if (qm>sqrt(log(1e6)/a)) qm = sqrt(log(1e6)/a); q = sqrt(-log(1.0-rand01()*(1.0-exp(-a*qm*qm)))/a); q_v = q*K2V; /* scattering possible ??? */ arg = q_v/(2.0*v); if(arg<1.0 && rand01()>transmr) { theta = asin(arg); /* Bragg scattering law */ d_phi = 2*PI*rand01(); vec_prod(axis_x, axis_y, axis_z, vx, vy, vz, 0, 1, 0); rotate(tmp_vx, tmp_vy, tmp_vz, vx, vy, vz, 2*theta, axis_x, axis_y, axis_z); rotate(vout_x, vout_y, vout_z, tmp_vx, tmp_vy, tmp_vz, d_phi, vx, vy, vz); vx = vout_x; vy = vout_y; vz = vout_z; if(!box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth)) fprintf(stderr, "SANS_Guinier: FATAL ERROR: Did not hit box from inside.\n"); } SCATTER; } %} MCDISPLAY %{ double radius = 0; double h = 0; { 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); } %} END