/*******************************************************************************
*
* McStas, version 1.0, released October 26, 1998
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Risoe National Laboratory, Roskilde, Denmark
*
* Component: Mon_2foc
*
* %Identification
* Written by: PL
* Date: Feb. 12,1999
* Origin: Uni. Gottingen (Germany)
* version: 1.0
*
* Double bent monochromator with multiple slabs
* 
* %Description
* 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).
*
* Added double bent feature by: Peter Link Feb. 12,1999
* corrected bug in rotation of v-coords: Peter Link Sep. 24. 1999 
*
* %Parameters
* INPUT PARAMETERS:
*
* zwidth:  horizontal width of an individual slab (m)
* ywidth:  vertical heigth of an individual slab (m)
* gap:     typical gap  between adjacent slabs (m)
* 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)
*
* %Link 
* <a href="http://neutron.risoe.dk/mcstas/support/misc/Mon_2foc.txt">Additional note</a> from <a href="mailto:plink@physik.tu-muenchen.de">Peter Link</a>.
*
* %End
*******************************************************************************/

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,tmp4,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,vxp,vyp,vzp;
    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));
        if (RH != 0) tilth = asin((col-(NH+1)/2)*(zwidth+gap)/RH);
	else tilth=0;
        if (RV != 0) tiltv = -asin((row-(NV+1)/2)*(ywidth+gap)/RV);
	else tiltv=0;

        /* rotate with tilth and tiltv */

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

          
        /* First: scattering in plane */
        /* theta0 = atan2(vx,vz);  neutron angle to slab Risoe version */
      
	v = sqrt(vxp*vxp+vyp*vyp+vzp*vzp);
	theta0 = asin(vxp/v);                /* correct neutron angle to slab */

        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*vxp - sn*vzp; 
          vyp = vyp;
          /* vz = cs*vz + sn*vx; diese Zeile wurde durch die folgende ersetzt */
          tmp4 = vyp/vzp;  /* korrigiert den schrägen Einfall aufs Plättchen  */
          vzp = cs*(-vyp*sin(tmp4)+vzp*cos(tmp4)) + sn*vxp;  
          vxp = tmp2;

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

          phi = atan2(vyp,vzp);  /* out-of plane angle */
          dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm();  /* MC choice: */
          /* Vertical angle of the crystallite */
          vyp = vzp*tan(phi+2*dphi*sin(theta));
          vratio = v/sqrt(vxp*vxp+vyp*vyp+vzp*vzp);
          vzp = vzp*vratio;
          vyp = vyp*vratio;                             /* Renormalize v */
          vxp = vxp*vratio;
	}
	/* rotate v coords back */
	vx = vxp*csb*csa-vyp*snb*csa+vzp*sna;
	vy = vxp*snb+vyp*csb;
	vz = -vxp*csb*sna+vyp*snb*sna+vzp*csa;    
	v=sqrt(vx*vx+vy*vy+vz*vz); 
      }
      else
      {
        x = old_x; y = old_y; z = old_z; t = old_t;
      }
    }
  %}

MCDISPLAY
%{
  double zmin,zmax,ymin,ymax;
  int ih,iv;

  magnify("xy");
  for(ih = 0; ih < NH; ih++)
  {
    for(iv = 0; iv < NV; iv++)
    {
      zmin = (zwidth+gap)*(ih-NH/2.0)+gap/2;
      zmax = zmin+zwidth;
      ymin = (ywidth+gap)*(iv-NV/2.0)+gap/2;
      ymax = ymin+ywidth;
      multiline(5, 0.0, (double)ymin, (double)zmin,
               0.0, (double)ymax, (double)zmin,
               0.0, (double)ymax, (double)zmax,
               0.0, (double)ymin, (double)zmax,
               0.0, (double)ymin, (double)zmin);
     }
   }
%}

END