/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2015, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: PerfectCrystal
*
* %I
* Written by: Markus Appel
* Date: 2015-12-21
* Origin: ILL / FAU Erlangen-Nuernberg
* Based on a perfect crystal component by: Miguel A. Gonzalez, A. Dianoux June 2013 (ILL)
*
* Changelog:
* Version 1.1
*  - BUGFIX: correct neutron energy shift in Doppler mode
*  - added option 'debyescherrer' to select analyzer geometry
*  - added option 'facette' to approximate analyzer sphere by small, flat crystals
*
* Version 1.0
*  - inital release
*
*
* %D
* Perfect crystal component, primarily for use as monochromator and analyzer in
* backscattering spectrometers. Reflection from a single Bragg reflex of a flat or
* spherical surface is simulated using a Darwin, Ewald or Gaussian profile. Doppler
* movement of the monochromator is supported as well.
*
* Orientation of the crystal surface is *different* from other monochromator components!
* Gravitational energy shift for tall analyzers should work in principle, but it not tested yet.
* See the parameter description on how to define the geometry and properties.
*
* [1] Website for Backscattering Spectroscopy: http://www.ill.eu/sites/BS-review/index.htm
*
* Examples:
* IN16B (ILL) Si111 large angle analyzers (approximate dimensions):
* PerfectCrystal(radius=2, lambda=6.2708, sigma=0.24e-3,
*                ttmin=40, ttmax=165, phimin=-45, phimax=+45, centerfocus=1)
*
* IN16B (ILL) Si111 Doppler monochromator:
* PerfectCrystal(radius=2.165, lambda=6.2708, sigma=0.24e-3,
*                width = 0.500, height = 0.250, centerfocus=0,
*                speed = 4.7, amplitude = 0.075, exclusive=1)
*
* %P
* (A) Size, shape and position:
* =============================
* radius: [m]                  Radius of curvature, set to 0 for a flat surface. Default: 0
* centerfocus: [0/1]           Component origin is the center of the analyzer sphere if 1, if set to 0 the origin is on the analyzer surface. Default: 0
*
* option 1
* ttmin: [deg]                 
* ttmax: [deg]                 analyzer coverage angle in horizontal xz-plane between -180 and 180
* phimin: [deg]                
* phimax: [deg]                vertical analyzer coverage between -90 and 90 (-180 and 180 if debyescherrer==1)
* option 2
* tt0: [deg]                   
* ttwidth: [deg]               horizontal coverage as center and full width
* phi0: [deg]                  
* phiwidth: [deg]              vertical coverage as center and full width
* option 3
* tt0: [deg]                   angular center position in the horizontal plane (only if centerfocus==1)
* width: [m]                   absolute width
* phi0: [deg]                  angular vertical position (only if centerfocus==1)
* height: [m]                  absolute height (only if debyescherrer==0)
*
* debyescherrer: [-180...180]  (0/1)   bend analyzer following a Debye-Scherrer ring along scattering angle tt (twotheta) (phi = 
* facette: [m]                 width of square crystal facettes arranged on the spherical surface (set 0 to disable). Default: 0 Warning: Facettation will fail at the poles along +-y axis.
* facette_xi: [deg]            random misalignemt of each facette. Default: 0
*
* (B) Neutron optics
* ===================
* tau: [A^-1]                  Scattering vector of the reflex (sometimes also called Q...)
* lambda: [A]                  Alternatively to tau: backscattered wavelength
* R0: []                       Peak reflectivity. Default: 1
*
* Ewald/Darwin mode
* dtauovertau: [1]             Plateau width of the Ewald/Darwin curve (see 
* dtauovertau_ext: []          Relative variation of tau (randomized for each trajectory, full width)
* ewald: [0/1]                 Use Ewald curve if 1, Darwin curve if 0. Default: 1 (Ewald)
*
* Gaussian mode
* sigma: [meV]                 Width of the Gaussian reflectivity curve in meV (corresponding to the energy resolution). (The width will be transformed and the Gaussian is actually calculated in k-space.)
*
* Mirror mode
* ismirror: [0/1]              Simply reflect all neutrons if 1. Good for debugging/visualization. Default: 0
*
* (C) Doppler movement for monochromators (sinusoidal speed profile):
* ===================================================================
* speed: [m/s]                 Maximum Doppler speed. The actual monochromator velocity is randomized between -speed and +speed. Default: 0
* amplitude: [m]               Amplitude of the Doppler movement. Default: 0
*
* smartphase: [0/1]            Optimize Doppler phase for better MC efficiency if set to 1. *WARNING:* Experimental option! Always compare to a simulation without smartphase. Better do not use smartwidth with Ewald/Darwin curves due to their endless tails. Default: 0
* smartwidth: [meV]            Half width of the possible energy reflection window for smartphase. Default: 5*sigma OR 10*2*dtauovertau*E0
*
* (D) Miscellaneous
* ==================
* exclusive: [0/1]             If set to 1, absorb all neutrons that missed the monochromator/analyzer surface. Default: 0
* transmit: [0...1]            Monte-Carlo probability of transmitting an event through the monochromator/analyzer surface. (Events with R=1.0 for DarwinE/Ewald curves are always reflected!). Default: 0
*
* Output parameters
* ==================
* tt: [deg]                    Position where the neutron hit (or missed) the analyzer sphere
* phi: [deg]                   Position where the neutron hit (or missed) the analyzer sphere
* xi: [deg]                    Reflection angle between neutron trajectory and analyzer surface normal
* phid: [rad]                  Actual doppler phase during reflection [0,2*PI]
* vd: [m/s]                    Actual doppler speed during reflection
* zd: [m]                      Actual doppler displacement during reflection (zd>0 => displaced towards inc. beam from -z)
* v0: [m/s]                    Backscattered neutron velocity
* E0: [meV]                    Backscattered neutron energy
* R: [0...1]                   Reflectivity value (on Gauss/Darwin/Ewald curve)
* eps: [see 1]                 'abs(y)' from the Darwin/Ewald reflectivity formula 
*
* %E
*******************************************************************************/

DEFINE COMPONENT PerfectCrystal

SETTING PARAMETERS ( ttmin=NAN, ttmax=NAN, tt0=NAN, ttwidth=NAN, width=NAN,
phimin=NAN, phimax=NAN, phi0=NAN, phiwidth=NAN, height=NAN, debyescherrer=0,
facette=0, facette_xi=0, centerfocus=0, radius=0, tau=NAN, lambda=NAN, R0=1.0,
dtauovertau=NAN, dtauovertau_ext=0,
ewald=1,sigma=NAN,ismirror=0,speed=0,amplitude=0,smartphase=0,
smartwidth=NAN, exclusive=0, transmit=0, verbose=0 )

SHARE
%{
// convert between spherical (r,tt,phi) and cartesian (x,y,z) coordinates (angles in deg)
   // debyescherrer swaps axes
   void sph2cart(double *x,double *y,double *z, double r, double tt, double phi, int debyescherrer)
   {
      tt *= DEG2RAD;
      phi *= DEG2RAD;
      if (debyescherrer) {
         double sintt = sin(tt);
         *x = - r * cos(phi) * sintt;
         *y = r * sin(phi) * sintt;
         *z = r * cos(tt);
      } else {
         double cosphi = cos(phi);
         *x = - r * cosphi * sin(tt);
         *y = r * sin(phi);
         *z = r * cosphi * cos(tt);
      }
   }

   /*******************************************************************************
   * grandvec_target_circle: Choose random direction towards target at (x,y,z)
   * with given radius and gaussian area distribution.
   * If radius is zero, choose random direction in full 4PI, no target.
   ******************************************************************************/
   void
   grandvec_target_circle(double *xo, double *yo, double *zo, double *solid_angle,
                  double xi, double yi, double zi, double radius)
   {
     double l2, phi, theta, nx, ny, nz, xt, yt, zt, xu, yu, zu;

     if(radius == 0.0)
     {
       /* No target, choose uniformly a direction in full 4PI solid angle. */
       theta = acos (1 - rand0max(2));
       phi = rand0max(2 * PI);
       if(solid_angle)
         *solid_angle = 4*PI;
       nx = 1;
       ny = 0;
       nz = 0;
       yi = sqrt(xi*xi+yi*yi+zi*zi);
       zi = 0;
       xi = 0;
     }
     else
     {
       double costheta0;
       l2 = xi*xi + yi*yi + zi*zi; /* sqr Distance to target. */
       costheta0 = sqrt(l2/(radius*radius+l2));
       if (radius < 0) costheta0 *= -1;
       if(solid_angle)
       {
         /* Compute solid angle of target as seen from origin. */
           *solid_angle = 2*PI*(1 - costheta0);
       }

       /* Now choose point uniformly on circle surface within angle theta0 */
       double costheta;
       costheta = (1 - fabs(randnorm()*(1 - costheta0)) );

       if (costheta < -1)
         costheta = -1;

       theta = acos(costheta); /* radius on circle */
       phi = rand0max(2 * PI); /* rotation on circle at given radius */
       /* Now, to obtain the desired vector rotate (xi,yi,zi) angle theta around a
          perpendicular axis u=i x n and then angle phi around i. */
       if(xi == 0 && zi == 0)
       {
         nx = 1;
         ny = 0;
         nz = 0;
       }
       else
       {
         nx = -zi;
         nz = xi;
         ny = 0;
       }
     }

     /* [xyz]u = [xyz]i x n[xyz] (usually vertical) */
     vec_prod(xu,  yu,  zu, xi, yi, zi,        nx, ny, nz);
     /* [xyz]t = [xyz]i rotated theta around [xyz]u */
     rotate  (xt,  yt,  zt, xi, yi, zi, theta, xu, yu, zu);
     /* [xyz]o = [xyz]t rotated phi around n[xyz] */
     rotate (*xo, *yo, *zo, xt, yt, zt, phi, xi, yi, zi);
   } /* randvec_target_circle */


%}
DECLARE
%{
   // official output variables
   double tt;
   double phi;
   double xi;
   double phid;
   double vd;
   double zd;
   double v0;
   double vmin;
   double vmax;
   double vperp;
   double E0;
   double R;
   double eps;

   // internal stuff
   double sin_facette_xi;
%}
INITIALIZE
%{
   // checks and balances ...
   if ( radius < 0 )
   {
      MPI_MASTER(fprintf(stderr,
         "PerfectCrystal: %s ERROR: Invalid parameters: negative radius\n",
         NAME_CURRENT_COMP););
      exit(-1);
   }

   //******************************************
   // position and size
   if ( !radius )
   {
      // flat analyzer surface
      if ( isnan(width) || isnan(height) )
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s ERROR: Invalid parameters: Need width and height for radius==0\n",
            NAME_CURRENT_COMP););
            exit(-1);
      }
   }
   else
   {
      // curved analyzer surface

      //******************************************
      // Determine analyzer width

      // ttmin / ttmax
      if ( !isnan(ttmin) && !isnan(ttmax) && isnan(tt0) && isnan(ttwidth) && isnan(width) )
      {
         // all right.
      }
      // tt0 / ttwidth
      else if ( isnan(ttmin) && isnan(ttmax) && !isnan(ttwidth) && isnan(width) )
      {
         if ( isnan(tt0) )
         {
             MPI_MASTER(fprintf(stderr,
               "PerfectCrystal: %s WARNING: Missing parameter: tt0 set to zero\n",
               NAME_CURRENT_COMP););
            tt0 = 0;
         }
         ttmin = tt0 - ttwidth/2.0;
         ttmax = tt0 + ttwidth/2.0;

      }
      // tt0 / width
      else if ( isnan(ttmin) && isnan(ttmax) && isnan(ttwidth) && !isnan(width) )
      {
         if ( isnan(tt0) )
         {
             MPI_MASTER(fprintf(stderr,
               "PerfectCrystal: %s WARNING: Missing parameter: tt0 set to zero\n",
               NAME_CURRENT_COMP););
            tt0 = 0;
         }
         ttwidth = 2.0 * asin( width / 2.0 / radius ) * RAD2DEG;
         ttmin = tt0 - ttwidth/2.0;
         ttmax = tt0 + ttwidth/2.0;
      }
      else
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s ERROR: Invalid parameters: Cannot determine analyzer width/tt\n",
            NAME_CURRENT_COMP););
         exit(-1);
      }

      //******************************************
      // Determine analyzer height

      // phimin / phimax
      if ( !isnan(phimin) && !isnan(phimax) && isnan(phi0) && isnan(phiwidth) && isnan(height) )
      {
         // all right, nothing to do.
      }
      // phi0 / phiwidth
      else if ( isnan(phimin) && isnan(phimax) && !isnan(phiwidth) && isnan(height) )
      {
         if ( isnan(phi0) )
         {
             MPI_MASTER(fprintf(stderr,
               "PerfectCrystal: %s WARNING: Missing parameter: phi0 set to zero\n",
               NAME_CURRENT_COMP););
            phi0 = 0;
         }
         phimin = phi0 - phiwidth/2.0;
         phimax = phi0 + phiwidth/2.0;
      }
      // phi0 / height
      else if ( isnan(phimin) && isnan(phimax) && isnan(phiwidth) && !isnan(height) )
      {
         if ( isnan(phi0) )
         {
             MPI_MASTER(fprintf(stderr,
               "PerfectCrystal: %s WARNING: Missing parameter: phi0 set to zero\n",
               NAME_CURRENT_COMP););
            phi0 = 0;
         }
         phiwidth = 2.0 * asin( height / 2.0 / radius ) * RAD2DEG;
         phimin = phi0 - phiwidth/2.0;
         phimax = phi0 + phiwidth/2.0;
      }
      else
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s ERROR: Invalid parameters: Cannot determine analyzer height/phi\n",
            NAME_CURRENT_COMP););
         exit(-1);
      }
   }

   if ( centerfocus && !radius )
   {
      MPI_MASTER(fprintf(stderr,
         "PerfectCrystal: %s ERROR: Invalid parameters: centerfocus doesn't make sense with radius==0\n",
         NAME_CURRENT_COMP););
      exit(-1);
   }

   //******************************************
   // neutron optics
   if ( !ismirror )
   {
      if ( isnan(lambda) == isnan(tau) )
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s ERROR: Invalid parameters: provide either tau or lambda\n",
            NAME_CURRENT_COMP););
         exit(-1);
      }

      if ( isnan(tau) )
         tau = 4*PI / lambda;

      v0 = tau / 2.0 * K2V;
      E0 = SQR(v0) * VS2E;

      if ( isnan(dtauovertau) == isnan(sigma) )
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s ERROR: Invalid parameters: provide either sigma or dtauovertau, or switch on ismirror\n",
            NAME_CURRENT_COMP););
         exit(-1);
      }
   }
   else
   {
      if ( !isnan(dtauovertau) || !isnan(sigma) )
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s WARNING: dtauovertau and/or sigma is ignored with ismirror==1\n",
            NAME_CURRENT_COMP););
      }
   }


   if ( !( R0>=0 && R0<=1.0 ) )
   {
      MPI_MASTER(fprintf(stderr,
         "PerfectCrystal: %s ERROR: Invalid parameter: R0 must be between 0 and 1\n",
         NAME_CURRENT_COMP););
      exit(-1);
   }

   if ( transmit < 0.0 || transmit > 1.0 )
   {
      MPI_MASTER(fprintf(stderr,
         "PerfectCrystal: %s ERROR: Invalid parameter: transmit must be between 0 and 1\n",
         NAME_CURRENT_COMP););
      exit(-1);
   }

   // transform Gaussian sigma from energy (meV) to neutron velocity (m/s)
   if ( !isnan(sigma) )
      sigma /= 3.29106e-3 * tau;   // constant is hbar/2 in appropriate units

   // transform smartwidth from energy (meV) to neutron velocity (m/s)
   if ( !isnan(smartwidth) )
      smartwidth /= 3.29106e-3 * tau;   // constant is hbar/2 in appropriate units

   // set standard widths for smartphase
   if ( isnan(smartwidth) && !isnan(sigma) )
      smartwidth = 5.0 * sigma;

   if ( isnan(smartwidth) && !isnan(dtauovertau) )
      smartwidth = 10.0 * dtauovertau * v0;

   if ( smartphase && !isnan(dtauovertau) )
   {
       MPI_MASTER(fprintf(stderr,
         "PerfectCrystal: %s WARNING: Using smartphase for ewald/darwin curves is probably a bad idea, because these curves have very long wings which will be cut off!!!\n",
         NAME_CURRENT_COMP););
   }

   if (facette_xi)
   {
      if (facette_xi < 0)
      {
         MPI_MASTER(fprintf(stderr,
            "PerfectCrystal: %s ERROR: facette_xi must be >= 0\n",
            NAME_CURRENT_COMP););
         exit(-1);
      }

      sin_facette_xi = sin(DEG2RAD*facette_xi);
   }

   // talk to the user ...
   if ( verbose )
   {
      #define PRINTVAR(str,val) fprintf(stderr,"%s --- %s : %g\n",NAME_CURRENT_COMP,str,val);
      MPI_MASTER(
         PRINTVAR("ttmin",ttmin);
         PRINTVAR("ttmax",ttmax);
         PRINTVAR("ttwidth",ttwidth);
         PRINTVAR("width",width);
         PRINTVAR("phimin",phimin);
         PRINTVAR("phimax",phimax);
         PRINTVAR("phi0",phi0);
         PRINTVAR("phiwidth",phiwidth);
         PRINTVAR("height",height);
         PRINTVAR("centerfocus",centerfocus);
         PRINTVAR("debyescherrer",debyescherrer);
         PRINTVAR("radius",radius);
         PRINTVAR("facette",facette);
         PRINTVAR("facette_xi",facette_xi);
         PRINTVAR("sin_facette_xi",sin_facette_xi);
         PRINTVAR("tau",tau);
         PRINTVAR("lambda",lambda);
         PRINTVAR("v0",v0);
         PRINTVAR("E0",E0);
         PRINTVAR("dtauovertau",dtauovertau);
         PRINTVAR("dtauovertau_ext",dtauovertau_ext);
         PRINTVAR("ewald",ewald);
         PRINTVAR("R0",R0);
         PRINTVAR("sigma[m/s]",sigma);
         PRINTVAR("ismirror",ismirror);
         PRINTVAR("speed",speed);
         PRINTVAR("amplitude",amplitude);
         PRINTVAR("smartphase",smartphase);
         PRINTVAR("smartwidth[m/s]",smartwidth);
         PRINTVAR("exclusive",exclusive);
      );
   }

%}
TRACE
%{
   double dt1, dt2, q0mod;
   double nx,ny,nz;
   int missed = 0;

   // determine phase of doppler movement
   if ( (speed!=0) && smartphase)
   {
      // do something smart
      vmin = v0 - sqrt(SQR(vx)+SQR(vy)+SQR(vz)) - smartwidth;
      vmax = vmin + 2.0*smartwidth;

      if ( (vmin>speed) || (vmax<-speed) )
         ABSORB;

      if ( vmin < -speed )
         vmin = -speed;

      if ( vmax > speed )
         vmax = speed;

//      fprintf(stderr,"%g  %g\n",vmin,vmax);

      vd = vmin + (vmax-vmin)*rand01();
      phid = acos(vd/speed) + (rand01()<0.5?0:PI);
      zd = amplitude * sin(phid);
      p *= (vmax-vmin) / 2.0 / speed;
   }
   else if ( speed!=0 )
   {
      // random selection of phase
//      phid = rand01() * 2.0 * PI;
//      zd = amplitude * sin(phid);
//      vd = speed * cos(phid);

      // random selection of speed
      phid = acos(randpm1()) + (rand01()<0.5?0:PI);
      vd = speed * cos(phid);
      zd = amplitude * sin(phid);
   }
   else
   {
      // no movement
      zd = 0;
      vd = 0;
   }

   // Propagate to analyzer and determine surface normal
   if ( !radius )
   {
      // flat analyzer, use height and width

      // propoagate to surface
      if (!vz)
         ABSORB;
      dt2 = (-z-zd) / vz;
      PROP_DT(dt2);

      // see if the covered area is hit
      missed = !inside_rectangle(x,y,width,height);

      // surface normal in this case is simple
      nx=0;ny=0;nz=-1;
   }
   else
   {
      // spherical analyzer, use spherical coordinates ...

      // compute neutron path intersection with analyzer sphere
      if ( centerfocus )
         missed = !sphere_intersect(&dt1,&dt2,x,y,z+zd,vx,vy,vz,radius);
      else
         missed = !sphere_intersect(&dt1,&dt2,x,y,z+radius+zd,vx,vy,vz,radius);

      if ( !missed )
      {
         // propoagate to surface
         PROP_DT(dt2);

         // tt (twotheta) is calculated as in IN16B, positive values downstream rightwards
         // select coordinate system depending on 'debyescherrer' switch
         double zprime = z+(centerfocus?0:radius)+zd;
         if (debyescherrer) {
            tt = RAD2DEG * atan2( sqrt(SQR(x)+SQR(y)) , zprime );
            phi = RAD2DEG * atan2(y,-x);
         } else {
            tt = - RAD2DEG * atan2(x,zprime);
            phi = RAD2DEG * asin( y/sqrt( SQR(x) + SQR(y) + SQR(zprime) ) );
         }

         missed = (tt<ttmin || tt>ttmax || phi<phimin || phi>phimax);

         if ( !missed )
         {
            // analyzer surface normal
            if ( facette ) {
               // calculate center of the facette hit
               // always use spherical coordinates with y main axis (as in debyescherrer=0),
               // otherwise the pole will be in the analyzer center!
               // for the sake of confusion, use radians in this part.
               double tt1, ttfacette, ttn, phi1, phifacette, phin;
               if (debyescherrer) {
                  tt1 = - atan2(x,zprime);
                  phi1 = asin( y/sqrt( SQR(x) + SQR(y) + SQR(zprime) ) );
               } else {
                  tt1 = DEG2RAD * tt;
                  phi1 = DEG2RAD * phi;
               }

               phifacette = facette / radius;
               phin = floor( phi1 / phifacette + 0.5 ) * phifacette;
               ttfacette = facette / (radius*cos(phin));
               ttn = floor( tt1 / ttfacette + 0.5 ) * ttfacette;

               nx = cos(phin)*sin(ttn);
               ny = -sin(phin);
               nz = -cos(phin)*cos(ttn);

               if (facette_xi)
               {
                  grandvec_target_circle(&nx,&ny,&nz,NULL,nx,ny,nz,sin_facette_xi);
               }
            } else {
               nx = -x;
               ny = -y;
               nz = -z-(centerfocus?0:radius)-zd;
               NORM(nx,ny,nz);
            }
         }
      }
   }

   // Do the reflection
   if ( !missed )
   {
      // velocity vector projected on surface normal
      // in moving doppler frame
      vperp = scalar_prod(vx,vy,vz+vd,nx,ny,nz);

      // angle between surface normal and velocity (only used as output parameter for monitoring)
      xi = RAD2DEG * acos( - vperp / sqrt( SQR(vx) + SQR(vy) + SQR(vz+vd) ) );

      // energy selection
      if (!ismirror)
      {
         if ( !isnan(dtauovertau) )
         {
            // eps is actually abs(y)
            // vperp is negative!
            double this_tau = tau;
            if ( dtauovertau_ext )
            {
               this_tau *= 1.0 + 0.5*dtauovertau_ext*randpm1();
            }

            eps = fabs( 4.0*vperp*V2K/this_tau + 2.0 ) / dtauovertau;

            // Darwin/Ewald curve
            if ( eps > 1 )
            {
               if ( ewald )
               {
                  // energy selection with Ewald curve
                  R = 1.0 - sqrt(SQR(eps)-1.0) / eps;
               }
               else
               {
                  // energy selection with Darwin curve
                  R = eps - sqrt(SQR(eps)-1.0);
                  R *= R;
               }

               R *= R0;
            }
            else
            {
               R = R0;
            }
         }
         else
         {
            // Gauss curve (vperp is negative!)
            eps = fabs(v0 + vperp) / sigma;
            R = exp( -SQR(eps) / 2.0 );

            R *= R0;
         }
      }
      else
      {
         R = 1;
         eps = NAN;
      }

      if ( transmit && (R!=1) && (rand01()<transmit) )
      {
         p *= (1.0 - R) / transmit;
         SCATTER;
      }
      else
      {
         // reflection: vector kf = ki - qmod0  where qmod0 = 2*n*(n dot ki'),
         // n is the surface normal vector and ki' the incident wavevector
         // in the moving frame
         q0mod = 2.0 * vperp;
         // q0mod is now directly in velocity units (and contains energy shift due to Doppler effect!)

         vx -= q0mod*nx;
         vy -= q0mod*ny;
         vz -= q0mod*nz;
         if (R!=1)
            p *= R / (1.0-transmit);

         if (R<1e-10)
            ABSORB;
         else
            SCATTER;
      }

   }
   else if ( !exclusive )
   {
      R=NAN;
      xi=NAN;
      RESTORE_NEUTRON (INDEX_CURRENT_COMP,x, y, z, vx, vy, vz, t, sx, sy, sz, p);
   }
   else
      ABSORB;

%}

MCDISPLAY
%{
   int ttsegments;
   int phisegments;
   int steps;
//   const double centerspheresize = 0.1*radius;

   int i, itt, iphi, istep;
   double phi,tt,dphi,dtt;
   double x1,y1,z1,x2,y2,z2;
   double z0;

   if (debyescherrer) {
      ttsegments = 4;
      phisegments = 20;
      steps = 50;
   } else {
      ttsegments = 10;
      phisegments = 10;
      steps = 30;
   }


   // mark the center with a sphere
//   circle("xy",0,0,0,centerspheresize);
//   circle("xz",0,0,0,centerspheresize);
//   circle("yz",0,0,z0,centerspheresize);

   // draw crystal surface
   if ( radius )
   {
      z0 = (centerfocus ? 0 : radius);

      // curved surface
      dtt = (ttmax - ttmin) / (double)steps;
      dphi = (phimax - phimin) / (double)steps;

      for (itt=0; itt<=ttsegments; itt++) {
         tt = ttmin + (ttmax - ttmin) * (double)itt / (double)ttsegments;
         for (istep=0; istep<steps; istep++) {
            phi = phimin + dphi*istep;
            sph2cart(&x1,&y1,&z1,radius,tt,phi,debyescherrer);
            sph2cart(&x2,&y2,&z2,radius,tt,phi+dphi,debyescherrer);
            line(x1,y1,z1-z0,x2,y2,z2-z0);
         }
      }

      for (iphi=0; iphi<=phisegments; iphi++) {
         phi = phimin + (phimax - phimin) * (double)iphi / (double)phisegments;
         for (istep=0; istep<steps; istep++) {
            tt = ttmin + dtt*istep;
            sph2cart(&x1,&y1,&z1,radius,tt,phi,debyescherrer);
            sph2cart(&x2,&y2,&z2,radius,tt+dtt,phi,debyescherrer);
            line(x1,y1,z1-z0,x2,y2,z2-z0);
         }
      }

      // additional dashed moving one
      if ( amplitude )
      {

         sph2cart(&x1,&y1,&z1,radius,ttmin,phimin,debyescherrer);
         dashed_line(x1,y1,z1-z0-amplitude,x1,y1,z1-z0+amplitude,6);
         sph2cart(&x1,&y1,&z1,radius,ttmax,phimin,debyescherrer);
         dashed_line(x1,y1,z1-z0-amplitude,x1,y1,z1-z0+amplitude,6);
         sph2cart(&x1,&y1,&z1,radius,ttmin,phimax,debyescherrer);
         dashed_line(x1,y1,z1-z0-amplitude,x1,y1,z1-z0+amplitude,6);
         sph2cart(&x1,&y1,&z1,radius,ttmax,phimax,debyescherrer);
         dashed_line(x1,y1,z1-z0-amplitude,x1,y1,z1-z0+amplitude,6);

         for (i=-1;i<=1;i+=2)
         {
            z0 = (centerfocus ? 0 : radius) + i * amplitude;

            for (itt=0; itt<=ttsegments; itt++) {
               tt = ttmin + (ttmax - ttmin) * (double)itt / (double)ttsegments;
               for (istep=0; istep<steps; istep++) {
                  phi = phimin + dphi*istep;
                  sph2cart(&x1,&y1,&z1,radius,tt,phi,debyescherrer);
                  sph2cart(&x2,&y2,&z2,radius,tt,phi+dphi,debyescherrer);
                  dashed_line(x1,y1,z1-z0,x2,y2,z2-z0,1);
               }
            }

            for (iphi=0; iphi<=phisegments; iphi++) {
               phi = phimin + (phimax - phimin) * (double)iphi / (double)phisegments;
               for (istep=0; istep<steps; istep++) {
                  tt = ttmin + dtt*istep;
                  sph2cart(&x1,&y1,&z1,radius,tt,phi,debyescherrer);
                  sph2cart(&x2,&y2,&z2,radius,tt+dtt,phi,debyescherrer);
                  dashed_line(x1,y1,z1-z0,x2,y2,z2-z0,1);
               }
            }
         }
      }

   }
   else
   {
      // flat surface
      rectangle("xy",0.0,0.0,0.0,width,height);
      // additional moving one
      if ( amplitude )
      {
         dashed_line(+width/2.0,+height/2.0,-amplitude,+width/2.0,+height/2.0,+amplitude,6);
         dashed_line(+width/2.0,-height/2.0,-amplitude,+width/2.0,-height/2.0,+amplitude,6);
         dashed_line(-width/2.0,+height/2.0,-amplitude,-width/2.0,+height/2.0,+amplitude,6);
         dashed_line(-width/2.0,-height/2.0,-amplitude,-width/2.0,-height/2.0,+amplitude,6);
         for (i=-1;i<=1;i+=2)
            rectangle("xy",0.0,0.0,i * amplitude,width,height);
      }
   }

%}

END