Bender

Philipp Bernhardt snphbern at jerry.rommel.stw.uni-erlangen.de
Wed Jan 20 15:57:26 CET 1999


Hello,

The following routine is simulating a curved neutron guide (bender). It is bent to the negative
X axis. It behaves like a parallel guide in the Y axis. To get a better transmission coefficient, you can 
split the bender into channels. These channels are separated by partitions with the thickness of d.

Because the angle of reflection doesn't change, the routine 
calculates the reflection coefficent for the concave and, if necessary, for the convex wall only onces, 
together
with the number of reflections.
Nevertheless the exact position, the time, and the divergence is calculated at the end of the bender, so
there aren't any approximations. 
 
I have compared the results with analytic results in the case of an ideal reflection and with the 
programme haupt and found that they fit well. 
  
I would appreciate hearing about any comments or other test results obtained from
this routine.
 
Thank you
Philipp Bernhardt

Universität Erlangen/Nürnberg
Lehrstuhl für Kristallographie und Strukturphysik
D-91054 Erlangen
E-mail: philipp.bernhardt at krist.uni-erlangen.de

-------------- next part --------------
/**********************************************************************************************
*
* Component: Bender
*
* Written by: Philipp Bernhardt, Januar 10 1999
*
* Models a curved neutron guide. 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 r*Win, so you do not have to calculate the real exit coordinates and you do not
* need a new arm. The bender is bent to the negative X axis; it behaves like a 
* parallel guide in the Y axis.
* There is no warning, if the parameters are wrong, so please make sure that
* - w,h,r,d,Win,k are positive numbers 
* - k*d is smaller than the width w
* - Win is big enough to prevent, that a neutron can pass the bender without a reflection
* otherwise you will get wrong results  
*
* 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
* d:       (m)    Thickness of the partition, which separates the channels
* Win:     (rad)  Angle of the deflection of the whole bender (r*Win is the length of the bender)
* k:       (AA)   Number of channels 
* R0a,Qca,alphaa,ma,Wa:  Parameters, which are describing the mirror at concave side of the bender
* R0i,Qci,alphai,mi,Wi:  Parameters, which are describing the mirror at convex side of the bender
* R0s,Qcs,alphas,ms,Ws:  Parameters, which are describing the mirror at the top and bottom side 
*                        of the bender 
*
* Example values: w=0.05 h=0.12 r=250 d=0.001 Win=0.04 k=3 ma=3 mi=2 ms=2
************************************************************************************************/
 
DEFINE COMPONENT Bender
DEFINITION PARAMETERS (w,h,r,d,Win,k,R0a,Qca,alphaa,ma,Wa,R0i,Qci,alphai,mi,Wi,R0s,Qcs,alphas,ms,Ws)
SETTING PARAMETERS ()
OUTPUT PARAMETERS ()
STATE PARAMETERS (x, y, z, vx, vy, vz, t, s1, s2, p)
DECLARE
 %{
        double bk; 
 %}

INITIALIZE
 %{
	bk=(w+d)/k;                            /* width of one channel + thickness d of partition */
 %}

TRACE
 %{
	double erszei,aktwin,aktzei,zykzei,zykwin,ab,dab,dru,einwin,Q,arg,ref,innref,aeuref,aeuwin,innwin,R,T,Ta;
	PROP_Z0;

        if ((fabs(x)>w/2) || (fabs(y)>h/2))    /* neutron is missing the bender */
	   ABSORB;

        dru=floor((w/2-x)/bk)*bk-w/2;          
        ab=-dru-x;                             /* distance between neutron and the concave side of neutron */ 
        R=r-dru;                               /* radius of the channel for the current neutron */ 
        if (ab>bk-d)                           /* neutron is hitting the partition */ 
           ABSORB; 

 /* reflections in the XZ-plane */
	einwin=atan(vx/vz);                    /* divergence of the current neutron at the entrance */	
        dab=R-cos(einwin)*(R-ab);              /* maximal distance between neutron and concave side */
        aeuwin=acos((R-dab)/R);                /* angle of reflection at the concave side */ 
        zykwin=2*aeuwin;                        
        T=2*sqrt(R*R-(R-dab)*(R-dab))/sqrt(vx*vx+vz*vz);  /* time period between two reflections at */ 
                                                          /* the concave side */
        innref=1.0;
        Q=2*V2K*sqrt(vx*vx+vz*vz)*sin(aeuwin); 
        if (Q<=Qca) 
           aeuref=R0a;
        else {
           arg=(Q-ma*Qca)/Wa;
           if (arg<10) 
              aeuref=.5*R0a*(1-tanh(arg))*(1-alphaa*(Q-Qca));   /* reflectivity at the concave side */ 
           else
              ABSORB; }  

        if (dab>bk-d) {                               /* there are also zick-zack-reflections */          
           innwin=acos(R*cos(aeuwin)/(R-bk+d));       /* angle of reflection at the convex side */ 
           Q=2*V2K*sqrt(vx*vx+vz*vz)*sin(innwin);
           if (Q<=Qci) 
              innref=R0i;
           else {
              arg=(Q-mi*Qci)/Wi;
              if (arg<10) 
                 innref=.5*R0i*(1-tanh(arg))*(1-alphai*(Q-Qci));     /* reflectivity at the convex side */
              else
                 ABSORB; }
           zykwin=2*(aeuwin-innwin); 
           T=2*sqrt(R*R+(R-bk+d)*(R-bk+d)-2*R*(R-bk+d)*cos(zykwin/2))/sqrt(vx*vx+vz*vz); }
        if (einwin<0) {                                               
           aktwin=Win-zykwin+aeuwin+einwin;
           p*=aeuref*innref; 
           Ta=sqrt((R-ab)*(R-ab)+(R-bk+d)*(R-bk+d)-2*(R-ab)*(R-bk+d)*cos(Win-aktwin-zykwin/2))/sqrt(vx*vx+vz*vz)+T/2; }
        else {
           aktwin=Win-aeuwin+einwin;
           p*=aeuref; 
           Ta=sqrt((R-ab)*(R-ab)+R*R-2*(R-ab)*R*cos(Win-aktwin))/sqrt(vx*vx+vz*vz); }
        while (aktwin>zykwin) {
           p*=aeuref*innref;
           aktwin-=zykwin; 
           Ta+=T; }
        if (aktwin>zykwin/2) {
           p*=innref;
           ab=(R*cos(aeuwin-zykwin+aktwin)-R+dab)/cos(aeuwin-zykwin+aktwin);
           Ta+=sqrt((R-ab)*(R-ab)+(R-bk+d)*(R-bk+d)-2*(R-ab)*(R-bk+d)*cos(aktwin-zykwin/2))/sqrt(vx*vx+vz*vz)+T/2;
           vx=sin(aeuwin-zykwin+aktwin)*sqrt(vx*vx+vz*vz);
           vz=vx/tan(aeuwin-zykwin+aktwin); }   
        else {
           ab=(R*cos(aeuwin-aktwin)-R+dab)/cos(aeuwin-aktwin);
           Ta+=sqrt((R-ab)*(R-ab)+R*R-2*(R-ab)*R*cos(aktwin))/sqrt(vx*vx+vz*vz);
           vx=-sin(aeuwin-aktwin)*sqrt(vx*vx+vz*vz);
           vz=-vx/tan(aeuwin-aktwin); }
        x=-ab-dru;
        z=r*Win;
        t+=Ta;

 /* reflections at top an bottom (Y axis) */ 
        if (vy!=0.0) {
           zykzei=h/fabs(vy);                           /* time period between two reflections */ 
           Q=2*V2K*fabs(vy);
           if (Q<=Qcs) 
              ref=R0s;
           else {
              arg=(Q-ms*Qcs)/Ws;
              if (arg<10) 
                 ref=.5*R0s*(1-tanh(arg))*(1-alphas*(Q-Qcs));        /* reflectivity at top and bottom */    
              else 
                 ABSORB; }
           erszei=(vy/fabs(vy)*h/2-y)/vy;                /* time till neutron hit the top or bottom side 
                                                            for the first time */  
           if (erszei<Ta) {
              p*=ref;
              vy*=-1;
              aktzei=Ta-erszei;
              while (aktzei>zykzei) {
                 p*=ref;
                 aktzei-=zykzei;
                 vy*=-1; }
              y=-vy/fabs(vy)*h/2+aktzei*vy; }
           else 
              y+=Ta*vy; }
%}

FINALLY
 %{
 %}

END



More information about the mcstas-users mailing list