/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright 1997-2002, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: Source_div * * %I * Written by: KL * Date: November 20, 1998 * Modified by: KL, 8 October 2001 * Origin: Risoe * * Neutron source with Gaussian or uniform divergence * * %D * The routine is a rectangular neutron source, which has a gaussian or uniform * divergent output in the forward direction. * The neutron energy is distributed between lambda0-dlambda and * lambda0+dlambda or between E0-dE and E0+dE. The flux unit is specified * in n/cm2/s/st/energy unit (meV or Angs). * In the case of uniform distribution (gauss=0), angles are uniformly distributed * between -focus_aw and +focus_aw as well as -focus_ah and +focus_ah. * For Gaussian distribution (gauss=1), 'focus_aw' and 'focus_ah' define the * FWHM of a Gaussian distribution. Energy/wavelength distribution is also * Gaussian. * * Example: Source_div(xwidth=0.1, yheight=0.1, focus_aw=2, focus_ah=2, E0=14, dE=2, gauss=0) * * %VALIDATION * Feb 2005: tested by Kim Lefmann (o.k.) * Apr 2005: energy distribution used in external tests of Fermi choppers (o.k.) * Jun 2005: wavelength distribution used in external tests of velocity selectors (o.k.) * Validated by: K. Lieutenant * * %BUGS * distribution is uniform in (hor. and vert.) angle (relative to moderator normal), * therefore not suited for large angles * * %P * xwidth: [m] Width of source * yheight: [m] Height of source * focus_aw: [deg] FWHM (Gaussian) or maximal (uniform) horz. width divergence * focus_ah: [deg] FWHM (Gaussian) or maximal (uniform) vert. height divergence * E0: [meV] Mean energy of neutrons. * dE: [meV] Energy half spread of neutrons. * lambda0: [Ang] Mean wavelength of neutrons (only relevant for E0=0) * dlambda: [Ang] Wavelength half spread of neutrons. * gauss: [0|1] Criterion: 0: uniform, 1: Gaussian distributions * flux: [1/(s cm 2 st energy_unit)] flux per energy unit, Angs or meV * * OUTPUT PARAMETERS: * sigmah: [rad] parameter 'sigma' of the Gaussian distribution for horizontal divergence * sigmav: [rad] parameter 'sigma' of the Gaussian distribution for vertical divergence * p_init: [1] normalisation factor 1/'neutron_count' * * %E *******************************************************************************/ DEFINE COMPONENT Source_div DEFINITION PARAMETERS () SETTING PARAMETERS (xwidth, yheight, focus_aw, focus_ah, E0=0.0, dE=0.0, lambda0=0.0, dlambda=0.0, gauss=0, flux=1) OUTPUT PARAMETERS () /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ DECLARE %{ double sigmah; double sigmav; double p_init; double dist; double focus_xw; double focus_yh; %} INITIALIZE %{ sigmah = DEG2RAD*focus_aw/(sqrt(8.0*log(2.0))); sigmav = DEG2RAD*focus_ah/(sqrt(8.0*log(2.0))); if (xwidth < 0 || yheight < 0 || focus_aw < 0 || focus_ah < 0) { printf("Source_div: %s: Error in input parameter values!\n" "ERROR Exiting\n", NAME_CURRENT_COMP); exit(0); } if ((!lambda0 && !E0 && !dE && !dlambda)) { printf("Source_div: %s: You must specify either a wavelength or energy range!\n ERROR - Exiting\n", NAME_CURRENT_COMP); exit(0); } if ((!lambda0 && !dlambda && (E0 <= 0 || dE < 0 || E0-dE <= 0)) || (!E0 && !dE && (lambda0 <= 0 || dlambda < 0 || lambda0-dlambda <= 0))) { printf("Source_div: %s: Unmeaningful definition of wavelength or energy range!\n ERROR - Exiting\n", NAME_CURRENT_COMP); exit(0); } /* compute distance to next component */ Coords ToTarget; double tx,ty,tz; ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+1),POS_A_CURRENT_COMP); ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget); coords_get(ToTarget, &tx, &ty, &tz); dist=sqrt(tx*tx+ty*ty+tz*tz); /* compute target area */ if (dist) { focus_xw=dist*tan(focus_aw*DEG2RAD); focus_yh=dist*tan(focus_ah*DEG2RAD); } p_init = flux*1e4*xwidth*yheight/mcget_ncount(); if (!focus_aw || !focus_ah) exit(printf("Source_div: %s: Zero divergence defined. \n" "ERROR Use non zero values for focus_aw and focus_ah.\n", NAME_CURRENT_COMP)); p_init *= 2*fabs(DEG2RAD*focus_aw*sin(DEG2RAD*focus_ah/2)); /* solid angle */ if (dlambda) p_init *= 2*dlambda; else if (dE) p_init *= 2*dE; %} TRACE %{ double E,lambda,v; double tan_h; double tan_v; double thetah; double thetav; p=p_init; z=0; t=0; x=randpm1()*xwidth/2.0; y=randpm1()*yheight/2.0; if(lambda0==0) { if (!gauss) { E=E0+dE*randpm1(); /* Choose from uniform distribution */ } else { E=E0+randnorm()*dE; } v=sqrt(E)*SE2V; } else { if (!gauss) { lambda=lambda0+dlambda*randpm1(); } else { lambda=lambda0+randnorm()*dlambda; } v = K2V*(2*PI/lambda); } if (gauss==1) { thetah = randnorm()*sigmah; thetav = randnorm()*sigmav; } else { thetah = randpm1()*focus_aw*DEG2RAD/2; thetav = randpm1()*focus_ah*DEG2RAD/2; } tan_h = tan(thetah); tan_v = tan(thetav); /* Perform the correct treatment - no small angle approx. here! */ vz = v / sqrt(1 + tan_v*tan_v + tan_h*tan_h); vy = tan_v * vz; vx = tan_h * vz; %} MCDISPLAY %{ multiline(5, -xwidth/2.0, -yheight/2.0, 0.0, xwidth/2.0, -yheight/2.0, 0.0, xwidth/2.0, yheight/2.0, 0.0, -xwidth/2.0, yheight/2.0, 0.0, -xwidth/2.0, -yheight/2.0, 0.0); if (dist) { dashed_line(0,0,0, -focus_xw/2,-focus_yh/2,dist, 4); dashed_line(0,0,0, focus_xw/2,-focus_yh/2,dist, 4); dashed_line(0,0,0, focus_xw/2, focus_yh/2,dist, 4); dashed_line(0,0,0, -focus_xw/2, focus_yh/2,dist, 4); } %} END