double-bent monochromator component

[Peter Link]

Dear Kristian,
on the basis of the standard monochromator component I tried to develop
a realistic double-bent monochromator. I need this for my future triple
axis machine and I know of such monochromators for example at focus
(PSI). The monochromator is made of n rows and m columns of equally
sized slabs. Due to the mechanics there is a small spacing between
adjacent slabs. The focussing is done by tilting the individual slabs
around their horizontal axe and turning the columns around the column 
vertical axe.
In my component I determine first the slab hit by the neutron, then
calculate the tilt horizontal and vertical from row and column number
with the focussing radius ( horizontal and vertical focussing is
independent) In order to use the same asumptions as in the normal
monochromator component I then introduce a coordinate transformation
into the rotated coords. At the end I transform back to the ccords of
the hole component given by the central slab. I performed some tests
with my instrument definition and I think it works, but may be I got
some bug which I've overseen and therefore I send it to you. May be
there is someone else interested to test it further. I know the for
very  small focussing radii ( compared to the size ) this can never work
and also the small angle approximation may be violated, but i.g. radii
are of the dimension of meters compared to sizes of tenth of cm.
Hoping that I've produce no nonsense, best regards Peter 
-- 
******************************************
Dr. P. Link
IPC Uni Göttingen
Außenstelle Garching
Neue Forschungs- Neutronenquelle Garching
ZBE- FRM II- Bau
85747 Garching

Tel. 089 2891 4622
Fax. 089 2895 4622
mailto: plink at physik.tu-muenchen.de
******************************************

-------------- next part --------------
/*******************************************************************************
*
* McStas, version 1.0, released October 26, 1998
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Risoe National Laboratory, Roskilde, Denmark
*
* Component: Mon_2foc
*
* Written by: KL, HMR  June 16, 1997
* Rewritten by: KL  Oct. 16, 1997
* Added double bent feature by: Peter Link Feb. 12,1999
* 
* Double bent monochromator which uses a small-mosaicity approximation as well as
* the approximation vy^2 << vz^2 + vx^2.
* Second order scattering is neglected.
* For an unrotated monochromator component, the crystal plane lies in the y-z
* plane (ie. parallel to the beam).
*
* INPUT PARAMETERS:
*
* zwidth:  (horizontal) width of an individual slab
* ywidth:  (vertical) heigth of an individual slab
* gap   :   typical gap  between adjacent slabs
* NH    :   number of slabs horizontal ( columns )
* NV    :   number of slabs vertical   ( rows )
* mosaich: Horisontal mosaic (FWHM) (arc minutes)
* mosaicv: Vertical mosaic (FWHM) (arc minutes)
* R0:      Maximum reflectivity (1)
* Q:       Scattering vector (AA-1)
* RV :     radius of vertical focussing (m)
* RH :     radius of horizontal focussing (m)
*
*******************************************************************************/

DEFINE COMPONENT Mon_2foc
DEFINITION PARAMETERS (zwidth, ywidth, gap, NH, NV, mosaich, mosaicv, r0, Q, RV, RH)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
  %{
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
  %}

TRACE
  %{
    double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
    double zmin,zmax,ymin,ymax,zp,yp,row,col;
	double tilth,tiltv;         /* used to calculate tilt angle of slab */
	double sna,snb,csa,csb;
	double old_x = x, old_y = y, old_z = z, old_t = t;
    double dt;
    
	if(vx != 0.0 && (dt = -x/vx) >= 0.0)
    {
      zmax = ((NH*(zwidth+gap))-gap)/2;
	  zmin = -1*zmax;
	  ymax = ((NV*(ywidth+gap))-gap)/2;
	  ymin = -1*ymax;
      y += vy*dt; z += vz*dt; t += dt; x = 0.0;
	  zp = fmod ( (z-zmin),(zwidth+gap) );
	  yp = fmod ( (y-ymin),(ywidth+gap) );
	  
	/* hit a slab or a gap ? */
    
    if (z>zmin && z<zmax && y>ymin && y<ymax && zp<zwidth && yp< ywidth)
    {
      
      col = ceil ( (z-zmin)/(zwidth+gap));
	  row = ceil ( (y-ymin)/(ywidth+gap));
	  tilth = asin((col-(NH+1)/2)*(zwidth+gap)/RH);
	  tiltv = -asin((row-(NV+1)/2)*(ywidth+gap)/RV);
	  
	  /* rotate with tilth and tiltv */

      sna = sin(tilth);
	  snb = sin(tiltv);
	  csa = cos(tilth);
	  csb = cos(tiltv);
	  vx = vx*csa*csb+vy*snb-vz*sna*csb;
	  vy = -vx*csa*snb+vy*csb+vz*sna*snb;
	  vz = vx*sna+vz*csa;  

	  
	  /* First: scattering in plane */
      theta0 = atan2(vx,vz);           /* neutron angle to slab */
      v = sqrt(vx*vx+vy*vy+vz*vz);
      theta = asin(Q2V*Q/(2.0*v));               /* Bragg's law */
      if(theta0 < 0)
        theta = -theta;
      tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
      if(tmp3 > DIV_CUTOFF)
      {
        x = old_x; y = old_y; z = old_z; t = old_t;
      }
      else
      {
        p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
        tmp1 = 2*theta;
        cs = cos(tmp1);
        sn = sin(tmp1);
        tmp2 = cs*vx - sn*vz; 
        vy = vy;
        vz = cs*vz + sn*vx; 
        vx = tmp2;

        /* Second: scatering out of plane. 
           Approximation is that Debye-Scherrer cone is a plane */

        phi = atan2(vy,vz);  /* out-of plane angle */
        dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm();  /* MC choice: */
        /* Vertical angle of the crystallite */
        vy = vz*tan(phi+2*dphi*sin(theta));
        vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
        vz = vz*vratio;
        vy = vy*vratio;                             /* Renormalize v */
        vx = vx*vratio;
      }
    /* rotate v coords back */
    vx = vx*csb*csa-vy*snb*csa+vz*sna;
	vy = vx*snb+vy*csb;
	vz = -vx*csb*sna+vy*snb*sna+vz*csa;    

	}
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
  %}
END