/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Bender
*
* %Identification
* Written by: Philipp Bernhardt
* Date: Februar 7 1999
* Origin: Uni. Erlangen (Germany)
*
* Models a curved neutron guide.
*
* %Description
* Models a curved neutron guide with cylindrical walls.
*
* Bender radius, entrance width and height are necessary input data. To define
* the bender, you may either enter the deviation angle 'Win' or the length 'l'.
* The bender may consist of 'k' vertical channels, separated by partitioning walls
* of thickness 'd'. Three different reflectivity profiles can be given: for outer
* walls, for inner walls and for the top and bottom walls. The partitioning walls
* have the same coating as the exterior walls.
*
* The entrance lies in the X-Y plane, centered on the Z axis. The neutrons will
* also leave the bender in the X-Y plane at the z-value l=r*Win, i.e. they are
* centred on (0,0,r*Win); they have an (average) flight direction along the z-axis.
* Therefore, the following component is adjacent, if positioned AT (0,0,r*Win)
* without rotation.
* So, seen from outside, it behaves like a straight guide along the Z axis. As a
* consequence, it is shown straight in 'mcdisplay'.
* This behaviour results from a co-ordinate transformation inside the component.
* It is done to facilitate its use. Neither rotation nor shift along the x-axis
* need to be calculated; a new arm is not necessary. Internally, the bender is
* bent to the negative X axis;
*
* Example:
* Bender of 120 mm height, 50 mm width, 250 m radius and 0.04 rad (or 2.292 deg) curvature
* not channeled, with a standard m=2 coating on
*
* Bender(w=0.05,h=0.12,r=250,Win=0.04,
* R0a=0.99,Qca=0.021,alphaa=6.07,ma=2,Wa=0.003,
* R0i=0.99,Qci=0.021,alphai=6.07,mi=2,Wi=0.003,
* R0s=0.99,Qcs=0.021,alphas=6.07,ms=2,Ws=0.003)
*
* %BUGS
* Some users have reported potentially strange behaviours with this component.
* This component does not work with gravitation on.
*
* %Parameters
* INPUT PARAMETERS:
*
* w: [m] Width at the bender entry and exit
* h: [m] Height at the bender entry and exit
* r: [m] Radius of the bender
* Win: [rad] Angle of the deflection of the whole bender
* k: [1] Number of channels inside the bender
* d: [m] Thickness of one blade separating the channels
* R0a: [1] Low-angle reflectivity at the outer side of the bender
* Qca: [AA-1] Critical scattering vector
* alphaa: [AA] Slope of reflectivity
* ma: [1] m-value of material
* Wa: [AA-1] Width of supermirror cut-off
* R0i: [1] Low-angle reflectivity at the inner side of the bender
* Qci: [AA-1] Critical scattering vector
* alphai: [AA] Slope of reflectivity
* mi: [1] m-value of material
* Wi: [AA-1] Width of supermirror cut-off
* R0s: [1] Low-angle reflectivity at the top and bottom side of the bender
* Qcs: [AA-1] Critical scattering vector
* alphas: [AA] Slope of reflectivity
* ms: [1] m-value of material
* Ws: [AA-1] Width of supermirror cut-off
*
* Optional parameters:
* l: [m] length of bender l=r*Win
*
* OUTPUT PARAMETERS:
* bk: [m] Width of 1 channel + 1 separating blade
* mWin: [rad] Angle of the deflection of the whole bender
*
* %L
* Additional note from Philipp Bernhardt.
*
* %End
*******************************************************************************/
DEFINE COMPONENT Bender
DEFINITION PARAMETERS()
SETTING PARAMETERS (w,h,r,Win=0.04,k=1,d=0.001,l=0,
R0a=0.99,Qca=0.021,alphaa=6.07,ma=2,Wa=0.003,
R0i=0.99,Qci=0.021,alphai=6.07,mi=2,Wi=0.003,
R0s=0.99,Qcs=0.021,alphas=6.07,ms=2,Ws=0.003)
OUTPUT PARAMETERS ()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "ref-lib"
%}
DECLARE
%{
double bk;
double mWin;
%}
INITIALIZE
%{
if (r <0)
{ fprintf(stderr,"Bender: error: %s: to bend in the other direction\n", NAME_CURRENT_COMP);
fprintf(stderr," rotate comp on z-axis by 180 deg.\n"); exit(-1); }
if (k*d > w)
{ fprintf(stderr,"Bender: error: %s has (k*d > w).\n", NAME_CURRENT_COMP);
exit(-1); }
if (w*h*r*Win*k == 0)
{ fprintf(stderr,"Bender: error: %s has one of w,h,r,Win,k null.\n", NAME_CURRENT_COMP);
exit(-1); }
/* width of one channel + thickness d of partition */
mWin = Win;
if (l!= 0 && r != 0) mWin = (double)l/(double)r;
bk=(w+d)/k;
if (mcgravitation) fprintf(stderr,"WARNING: Bender: %s: "
"This component produces wrong results with gravitation !\n",
NAME_CURRENT_COMP);
%}
TRACE
%{
int i,num,numa,numi;
double dru,ab,dab,R,Q,Ta,vpl;
double einmWin,ausmWin,zykmWin,aeumWin,innmWin,ref,innref,aeuref;
double einzei,auszei,zykzei;
/* does the neutron hit the bender at the entrance? */
PROP_Z0;
if ((fabs(x)bk-d)
{
double inn[] = {R0i, Qci, alphai, mi, Wi};
/* reflection coefficient at the convex side */
innmWin=acos((R-dab)/(R-bk+d));
Q=2.0*V2K*vpl*sin(innmWin);
StdReflecFunc(Q, inn, &innref);
}
/* divergence of the neutron at the exit */
zykmWin=2.0*(aeumWin-innmWin);
ausmWin=fmod(mWin+einmWin+aeumWin-innmWin
*(1.0+SIGN(einmWin)),zykmWin)-zykmWin/2.0;
ausmWin+=innmWin*SIGN(ausmWin);
/* number of reflections at the concave side */
numa=(mWin+einmWin+aeumWin-innmWin*(1.0+SIGN(einmWin)))/zykmWin;
/* number of reflections at the convex side */
numi=numa;
if (ausmWin*einmWin<0)
{
if (ausmWin-einmWin>0)
numi++;
else
numi--;
}
/* is the reflection coefficient too small? */
if (((numa>0) && (aeuref<=0)) || ((numi>0) && (innref<=0)))
ABSORB;
/* calculation of the neutron probability weight p */
for (i=1;i<=numa;i++)
p*=aeuref;
for (i=1;i<=numi;i++)
p*=innref;
/* time to cross the bender */
Ta=(2*numa*(tan(aeumWin)-tan(innmWin))
+tan(ausmWin)-tan(einmWin)
-tan(innmWin)*(SIGN(ausmWin)-SIGN(einmWin)))
*(R-dab)/vpl;
t+=Ta;
/* distance between neutron and concave side of channel at the exit */
ab=R-(R-dab)/cos(ausmWin);
/* calculation of the exit coordinates in the XZ-plane */
x=w/2.0-ab-dru;
z=r*mWin;
vx=sin(ausmWin)*vpl;
vz=cos(ausmWin)*vpl;
/*** reflections at top and bottom side (Y axis) ***/
if (vy!=0.0)
{
double s[] = {R0s, Qcs, alphas, ms, Ws};
/* reflection coefficent at the top and bottom side */
Q=2.0*V2K*fabs(vy);
StdReflecFunc(Q, s, &ref);
/* number of reflections at top and bottom */
einzei=h/2.0/fabs(vy)+y/vy;
zykzei=h/fabs(vy);
num=(Ta+einzei)/zykzei;
/* time between the last reflection and the exit */
auszei=fmod(Ta+einzei,zykzei);
/* is the reflection coefficient too small? */
if ((num>0) && (ref<=0))
ABSORB;
/* calculation of the probability weight p */
for (i=1;i<=num;i++) {
p*=ref;
vy*=-1.0; }
/* calculation of the exit coordinate */
y=auszei*vy-vy*h/fabs(vy)/2.0;
} /* if (vy!=0.0) */
SCATTER;
} /* if (dab>bk-d) */
else
ABSORB; /* hit separating walls */
}
else /* if ((fabs(x)