/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Tunneling_sample * * %I * Written by: Kim Lefmann * Date: 10.05.07 * Origin: Risoe * * A Double-cylinder shaped all-incoherent scatterer * with elastic, quasielastic (Lorentzian), and tunneling (sharp) * components. * * %D * A Double-cylinder shaped all-incoherent scatterer * with both elastic, quasielastic (Lorentzian), and tunneling (sharp) * components. No multiple scattering. Absorbtion included. * The shape of the sample may be a box with dimensions xwidth, yheight, zdepth. * The area to scatter to is a disk of radius 'focus_r' situated at the target. * This target area may also be rectangular if specified focus_xw and focus_yh * or focus_aw and focus_ah, respectively in meters and degrees. * The target itself is either situated according to given coordinates (x,y,z), * or defined with the relative target_index of the component to focus * to (next is +1). * This target position will be set to its AT position. When targeting to * centered components, such as spheres or cylinders, define an Arm component * where to focus to. * * The outgoing polarization is calculated as for nuclear spin incoherence: * P' = 1/3*P-2/3P = -1/3P * As above multiple scattering is ignored . * * Example: Tunneling_sample(thickness=0.001,radius=0.01,yheight=0.02,focus_r=0.035, * target_index=1) * * %P * INPUT PARAMETERS: * radius: [m] Outer radius of sample in (x,z) plane * yheight: [m] vert. dimension of sample, as a height * thickness: [m] Thickness of cylindrical sample in (x,z) plane * focus_r: [m] Radius of disk containing target. Use 0 for full space * target_index: [1] relative index of component to focus at, e.g. next is +1 * xwidth: horiz. dimension of sample, as a width [m] * zdepth: depth of sample [m] * focus_xw: horiz. dimension of a rectangular area [m] * focus_yh: vert. dimension of a rectangular area [m] * focus_aw: horiz. angular dimension of a rectangular area [deg] * focus_ah: vert. angular dimension of a rectangular area [deg] * sigma_abs:Absorbtion cross section pr. unit cell [barns] * sigma_inc:Total incoherent scattering cross section pr. unit cell [barns] * Vc: Unit cell volume [AA^3] * p_interact: MC Probability for scattering the ray; otherwise transmit [1] * f_QE: Fraction of quasielastic scattering [1] * f_tun: Fraction of tunneling scattering (f_QE+f_tun < 1) [1] * gamma: Lorentzian width of quasielastic broadening (HWHM) [meV] * E_tun: Tunneling energy [meV] * target_x: X-position of target to focus at [m] * target_y: Y-position of target to focus at [m] * target_z: Z-position of target to focus at [m] * * Variables calculated in the component * * V_my_s: Attenuation factor due to scattering [m^-1] * V_my_a: Attenuation factor due to absorbtion [m^-1] * * %L * Test * results (not up-to-date). * * %E *******************************************************************************/ DEFINE COMPONENT Tunneling_sample DEFINITION PARAMETERS () SETTING PARAMETERS (thickness=0, radius=0.01, focus_r = 0, p_interact = 1, f_QE=0, f_tun=0, gamma=0, E_tun=0, target_x = 0, target_y = 0, target_z = 0.235, focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, xwidth=0, yheight=0.05, zdepth=0, sigma_abs=5.08, sigma_inc=4.935, Vc=13.827, int target_index=0) OUTPUT PARAMETERS() /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ struct StructVarsV { double sigma_a; /* Absorption cross section per atom (barns) */ double sigma_i; /* Incoherent scattering cross section per atom (barns) */ double rho; /* Density of atoms (AA-3) */ double my_s; double my_a_v; char isrect; /* true when sample is a box */ double distance; /* when non zero, gives rect target distance */ double aw,ah; /* rectangular angular dimensions */ double xw,yh; /* rectangular metrical dimensions */ double tx,ty,tz; /* target coords */ }; %} DECLARE %{ struct StructVarsV VarsV; double ftun; double fQE; %} INITIALIZE %{ if (!xwidth || !yheight || !zdepth) /* Cannot define a rectangle */ if (!radius || !yheight) /* Cannot define a cylinder either */ exit(fprintf(stderr,"V_sample: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP)); else /* It is a cylinder */ VarsV.isrect=0; else /* It is a rectangle */ VarsV.isrect=1; VarsV.sigma_a=sigma_abs; VarsV.sigma_i=sigma_inc; VarsV.rho = (1/Vc); VarsV.my_s=(VarsV.rho * 100 * VarsV.sigma_i); VarsV.my_a_v=(VarsV.rho * 100 * VarsV.sigma_a); /* now compute target coords if a component index is supplied */ VarsV.tx= VarsV.ty=VarsV.tz=0; 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, &VarsV.tx, &VarsV.ty, &VarsV.tz); } else { VarsV.tx = target_x; VarsV.ty = target_y; VarsV.tz = target_z; } if (!(VarsV.tx || VarsV.ty || VarsV.tz)) printf("Tunneling_sample: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP); VarsV.distance=sqrt(VarsV.tx*VarsV.tx+VarsV.ty*VarsV.ty+VarsV.tz*VarsV.tz); /* different ways of setting rectangular area */ VarsV.aw = VarsV.ah = 0; if (focus_xw) { VarsV.xw = focus_xw; } if (focus_yh) { VarsV.yh = focus_yh; } if (focus_aw) { VarsV.aw = DEG2RAD*focus_aw; } if (focus_ah) { VarsV.ah = DEG2RAD*focus_ah; } /* Check that probabilities are positive and do not exceed unity */ if (f_tun<0) ftun=0; else ftun=f_tun; if(f_QE<0) fQE=0; else fQE=f_QE; if ((ftun+fQE)>1) { ftun=0; printf("Tunneling_sample: Sum of inelastic probabilities > 1. Setting f_tun=0"); if (fQE>1) { fQE=0; printf("Tunneling_sample: Probability fQE > 1. Setting fQE=0."); } } %} TRACE %{ double t0, t3; /* Entry/exit time for outer cylinder */ double t1, t2; /* Entry/exit time for inner cylinder */ double v; /* Neutron velocity */ double dt0, dt1, dt2, dt; /* Flight times through sample */ double l_full; /* Flight path length for non-scattered neutron */ double l_i, l_o=0; /* Flight path lenght in/out for scattered neutron */ double my_a=0; /* Velocity-dependent attenuation factor */ double solid_angle=0; /* Solid angle of target as seen from scattering point */ double aim_x=0, aim_y=0, aim_z=1; /* Position of target relative to scattering point */ double v_i, v_f, E_i, E_f; /* initial and final energies and velocities */ double dE; /* Energy transfer */ double scatt_choice; /* Representing random choice of scattering type */ int intersect=0; if (VarsV.isrect) intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth); else intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight); if(intersect) { if(t0 < 0) ABSORB; /* we already passed the sample; this is illegal */ /* Neutron enters at t=t0. */ if(VarsV.isrect) t1 = t2 = t3; else if(!thickness || !cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius-thickness, yheight)) t1 = t2 = t3; dt0 = t1-t0; /* Time in sample, ingoing */ dt1 = t2-t1; /* Time in hole */ dt2 = t3-t2; /* Time in sample, outgoing */ v = sqrt(vx*vx + vy*vy + vz*vz); l_full = v * (dt0 + dt2); /* Length of full path through sample */ if (v) my_a = VarsV.my_a_v*(2200/v); if (p_interact >= 1 || rand01() dt0) dt += dt1; /* jump to 2nd side of cylinder */ PROP_DT(dt+t0); /* Point of scattering */ if ((VarsV.tx || VarsV.ty || VarsV.tz)) { aim_x = VarsV.tx-x; /* Vector pointing at target (anal./det.) */ aim_y = VarsV.ty-y; aim_z = VarsV.tz-z; } if(VarsV.aw && VarsV.ah) { randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, VarsV.aw, VarsV.ah, ROT_A_CURRENT_COMP); } else if(VarsV.xw && VarsV.yh) { randvec_target_rect(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, VarsV.xw, VarsV.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); scatt_choice = rand01(); /* chooses type of scattering */ v_i = v; /* Store initial velocity in case of inel. */ E_i = VS2E*v_i*v_i; if (scatt_choice<(fQE+ftun)) /* Inelastic choices */ { if (scatt_choice0) dE = E_tun; else dE = -E_tun; } E_f = E_i + dE; if (E_f <= 0) ABSORB; v_f = SE2V*sqrt(E_f); v = v_f; } vx *= v; vy *= v; vz *= v; if(!VarsV.isrect) { if(!cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight)) { /* ??? did not hit cylinder */ printf("FATAL ERROR: Did not hit cylinder from inside.\n"); exit(1); } dt = t3; /* outgoing point */ if(thickness && cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius-thickness, yheight) && t2 > 0) dt -= (t2-t1); /* Subtract hollow part */ } else { if(!box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth)) { /* ??? did not hit box */ printf("FATAL ERROR: Did not hit box from inside.\n"); exit(1); } dt = t3; } l_o = v*dt; /* trajectory after scattering point: absorption only */ p *= v/v_i*l_full*VarsV.my_s*exp(-my_a*(l_i+v_i/v*l_o)-VarsV.my_s*l_i); /* We do not consider scattering from 2nd part (outgoing) */ p /= 4*PI/solid_angle; p /= p_interact; /* Polarisation part (1/3 NSF, 2/3 SF) */ sx *= -1.0/3.0; sy *= -1.0/3.0; sz *= -1.0/3.0; SCATTER; } else /* Transmitting; always elastic */ { p *= exp(-(my_a+VarsV.my_s)*l_full); p /= (1-p_interact); } } %} MCDISPLAY %{ if (!VarsV.isrect) { 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); if (thickness) { double radius_i=radius-thickness; circle("xz", 0, yheight/2.0, 0, radius_i); circle("xz", 0, -yheight/2.0, 0, radius_i); line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0); line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0); line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i); line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i); } } else { 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