/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2003, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: SANS_AnySamp * * %I * Written by: Henrich Frielinghaus * Date: Sept 2004 * Origin: FZ-Juelich/FRJ-2/IFF/KWS-2 * * Sample for Small Angle Neutron Scattering. To be customized. * * %D * Sample for Small Angle Neutron Scattering. * May be customized for Any Sample model. * Copy this component and modify it to your needs. * Just give shape of scattering function. * Normalization of scattering will be done in INITIALIZE. * * Some remarks: * Scattering function should be a defined the same way in INITIALIZE and TRACE sections * There exist maybe better library functions for the integral. * * Here for comparison: Guinier. * Transmitted paths set to 3% of all paths. In this simulation method paths are * well distributed among transmission and scattering (equally in Q-space). * * 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_AnySamp(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] horizontal dimension of sample, as a width * yheight: [m] vertical dimension of sample, as a height * zdepth: [m] thickness of sample * * %Link * Sans_spheres component * * %E *******************************************************************************/ DEFINE COMPONENT SANS_AnySamp 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 %{ double isq; %} INITIALIZE %{ if (!xwidth || !yheight || !zdepth) { exit(fprintf(stderr,"SANS_AnySamp: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP)); } int iqmax=30000,iq=iqmax; // number of intervals double q,sq; double a=Rg*Rg/3.0; // internal for function sq isq = 0.0; // integral = 0 while (iq > 1) // start integrating with low intensities at large q { // MODIFIY HERE q = (iq-0.5)*qmax/iqmax; // q always slightly larger than 0 sq = exp(-a*q*q); // define this function in INITIALIZE and TRACE part sq*= q; isq+= sq; --iq; } isq*=qmax/iqmax; %} TRACE %{ double a,q,qm,sq,q_v; double transsim=0.03; // portion of paths for transmission 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>0.99 ) transmr = 0.99; 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 */ qm = qmax; // adjust maximal q if (qm > 2.0*v/K2V) qm = 2.0*v/K2V; // should not be totally wrong q = sqrt(rand01())*qm; // otherwise normalization with isq is wrong q_v = q*K2V; /* scattering possible ??? */ arg = q_v/(2.0*v); if(rand01()>transsim) { a = Rg*Rg/3.0; // MODIFIY HERE sq= exp(-a*q*q); // identical to INITIALIZE p*= sq*(qmax*qmax*0.5)*(1.0-transmr)/(1.0-transsim)/isq; 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_AnySamp: FATAL ERROR: Did not hit box from inside.\n"); } else { p*= transmr / transsim; } 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