/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: FermiChopper_ILL
*
* %Identification
*
* Written by: M. Poehlmann, C. Carbogno, H. Schober, E. Farhi
* Date:       May 2002
* Origin:     ILL Grenoble / TU Muenchen
* Modified by: K.Lieutenant, June 2005: added phase parameter. Comp validation.
* Modified by: E. Farhi, Jul 2008: uniformize parameter names ()
* Modified by: EF, Oct 2008: fix chopper orientation
* Modified by: EF, Mar 2009: fixed infinite recursion which may cause SEGV. Cleanup of code.
*
* Fermi Chopper with rotating frame.
*
* %D
* Models a Fermi chopper with optional supermirror coated blades
* supermirror facilities may be disabled by setting m = 0, R0=0
* Slit packages are straight. Chopper slits are separated by an infinitely
* thin absorbing material. The effective transmission (resulting from fraction
* of the transparent material and its transmission) may be specified.
* The chopper slit package width may be specified through the total width 'xwidth'
* of the full package or the width 'w' of each single slit. The other parameter
* is calculated by: xwidth = nslit*w.
*
* Example:
* FermiChopper_ILL(phase=-50.0, radius=0.04, nu=100,
*   yheight=0.08, w=0.00022475, nslit=200.0, R0=0.0,
*   Qc=0.02176, alpha=2.33, m=0.0, length=0.012, eff=0.95,
*   zero_time=0)
*
* Markus Poehlmann     <Markus.Poehlmann@ph.tum.de>
* Christian Carbogno   <carbogno@ph.tum.de>
* and Helmut Schober   <schober@ill.fr>
*
* %VALIDATION
* Apr 2005: extensive external test, most problems solved (cf. 'Bugs')
* Validated by: K. Lieutenant
*
* limitations:
* no blade width used
*
* %BUGS
* - overestimates peak width for long wavelengths
* - does not give the right pulse position, shape and width for slit widths below 0.1 mm
* - fails sometimes when using MPI
*
* %Parameters
* INPUT PARAMETERS:
*
* Geometrical chopper constants:
* radius: [m]   chopper cylinder radius
* yheight: [m]  Height of chopper
* nslit: [1]    number of chopper slits
* length: [m]   channel length of the Fermi chopper
* w: [m]        width of one chopper slit
* xwidth: [m]   optional total width of slit package
* nu:      [Hz]     chopper frequency
* eff:     [1]      efficiency = transmission x fraction of transparent material
* verbose: [1]      optional flag to display component statistics, use 1 or 3 (debuging)
*
* Supermirror constants:
* m:       [1]      m-value of material. Zero means completely absorbing.
* alpha:   [AA]     slope of reflectivity
* Qc:      [AA-1]   critical scattering vector
* W:      [AA-1]   width of supermirror cut-off
* R0:      [1]      low-angle reflectivity
*
* Constants to reset time of flight:
* zero_time: [1]    set time to zero: 0: no  1: once per half cycle
* phase:   [deg]    chopper phase at t=0
*
*
* %End
*****************************************************************************/

/* NOTE:
*   The initial component version (McStas version <= 1.12) was written
*   with an inverted coordinate frame orientation. This corresponds
*   to inverting the frequency and phase sign.
*/

DEFINE COMPONENT FermiChopper_ILL

SETTING PARAMETERS (phase=0, radius=0.04, nu=100,
yheight=0.08, w=0.00022475, nslit=200.0, R0=0.0,
Qc=0.02176, alpha=2.33, m=0.0, W=2e-3, length=0.012, eff=0.95,
zero_time=0, xwidth=0, verbose=0)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
#ifndef FCILL_TimeAccuracy
#define FCILL_TimeAccuracy 1e-8
#define FCILL_MAXITER      10
/* Definition of internal variable structure: all counters */
struct FermiChopper_ILL_struct {

/** other variables ********************************/
double omega, t0;  /* chopper rotation */
};


/**************** DECLARING FUNCTIONS ***************************************/

/*********** ORTHOGONAL TRANSFORMATION INTO ROTATING FRAME ******************/

/************ X - component ********************/
#pragma acc routine seq
double xstrich(double X, double Z, double T, double omega, double t0){
  return( X*cos(omega*(T-t0))+Z*sin(omega*(T-t0)) );
}

/************ Z - component ********************/
#pragma acc routine seq
double zstrich(double X, double Z, double T, double omega, double t0){
  return( Z*cos(omega*(T-t0))-X*sin(omega*(T-t0)) );
}

/*************************NUMERICAL METHODS *********************************/

/*************************** SECANT METHOD FOR... ***************************/

/****************************...X-component *********************************/
#pragma acc routine seq
double xsecant(double x, double z, double vx, double vz,
               double t, double dt, double d, double omega, double t0){

  double dt1     = 1;
  double counter = 0;
  double t1      = 0;
  double t2      = dt;
  double xr1     = xstrich(x,z,t, omega, t0)-d;
  double xr2     = xstrich(x+vx*t2,z+vz*t2,t+t2, omega, t0)-d;
  double sign;

  while ((fabs(dt1) > FCILL_TimeAccuracy) && (counter < FCILL_MAXITER) && (xr2-xr1)){
    counter++;
    dt1 = (t2-t1)*xr2/(xr2-xr1);
    t2  = t1;
    xr1 = xr2;
    t1 += dt1;
    xr2 = xstrich(x+vx*t1,z+vz*t1,t+t1, omega, t0)-d;
  }

  if(counter >= FCILL_MAXITER) t1 = -2;

  return(t1);
}


/****************************...Z-component *********************************/
#pragma acc routine seq
double zsecant(double x, double z, double vx, double vz,
               double t, double dt, double d, double omega, double t0) {

  double t1      = 0;
  double t2      = dt;
  double dt1     = 1;
  double counter = 0;
  double zr1     =  zstrich(x,z,t, omega, t0)-d;
  double zr2     =  zstrich(x+vx*t2,z+vz*t2,t+t2, omega, t0)-d;

  while ((fabs(dt1) > FCILL_TimeAccuracy) && (counter < FCILL_MAXITER) && (zr2-zr1)){
    counter++;
    dt1 = (t2-t1)*zr2/(zr2-zr1);
    t2  = t1;
    zr1 = zr2;
    t1 += dt1;
    zr2 = zstrich(x+vx*t1,z+vz*t1,t+t1, omega, t0)-d;
  }

  if(counter >= FCILL_MAXITER) t1=-1;

  return(t1);
}


/*************************** INTERPOLATION METHOD FOR... ********************/

/****************************...X-component *********************************/
#pragma acc routine seq
double xinterpolation(double x, double z, double vx, double vz,
                      double t, double dt, double d, double omega, double t0){

  double sign;
  double xr3=1, t3=0, t1=0, t2=dt, dt1=dt;
  double counter = 0;
  double xr1      =  xstrich(x,z,t, omega, t0)-d;
  double xr2      =  xstrich(x+vx*dt,z+vz*dt,t+dt, omega, t0)-d;

  while ((fabs(xr3) > FCILL_TimeAccuracy)&&(counter < FCILL_MAXITER)){
    counter++;
    t3 = (t1+t2)*0.5;
    xr3 = xstrich(x+(vx*(t3)),z+(vz*(t3)),t+t3, omega, t0)-d;
    xr2 = xstrich(x+(vx*(t2)),z+(vz*(t2)),t+t2, omega, t0)-d;
    if(xr2*xr3<0) t1=t3;
    else          t2=t3;
  }

  if(counter >= FCILL_MAXITER) t3=-1;

  return(t3);
}


/****************************...Z-component *********************************/
#pragma acc routine seq
double zinterpolation(double x, double z, double vx, double vz,
                      double t, double dt, double d, double omega, double t0){

  double counter = 0;
  double zr3=1,zr2=0,t3=0,t1=0,t2=dt;

  while ((fabs(zr3)>FCILL_TimeAccuracy)&&(counter<FCILL_MAXITER)) {
    counter++;
    t3 = (t1+t2)*0.5;
    zr3 = zstrich(x+(vx*(t3)),z+(vz*(t3)),t+t3, omega, t0)-d;
    zr2 = zstrich(x+(vx*(t2)),z+(vz*(t2)),t+t2, omega, t0)-d;
    if(zr2*zr3 < 0) t1=t3;
    else            t2=t3;
  }

  if(counter >= FCILL_MAXITER) t3=-1;

  return(t3);
}
#endif
%}

DECLARE
%{
  struct FermiChopper_ILL_struct FCVars;
%}

INITIALIZE
%{

/************************* INITIALIZE COUNTERS ******************************/

  int i;

/************************ CALCULATION CONSTANTS *****************************/
  FCVars.omega    = 2*PI*nu;
  if (nu && phase) FCVars.t0 = -phase/360.0/nu;

  /* check of input parameters */
  if (m < 0) m == 0;
  if (radius <= 0) {
    printf("FermiChopper_ILL: %s: FATAL: unrealistic cylinder radius radius=%g [m]\n", NAME_CURRENT_COMP, radius);
    exit(-1);
  }
  if (yheight <= 0) 
  	exit(printf("FermiChopper_ILL: %s: FATAL: unrealistic cylinder yheight =%g [m]\n", NAME_CURRENT_COMP, yheight));
  if (xwidth > 0 && xwidth < radius*2 && nslit > 0) {
    w = xwidth/nslit;
  }
  if (w <= 0) {
    printf("FermiChopper_ILL: %s: FATAL: Slits in the package have unrealistic width w=%g [m]\n", NAME_CURRENT_COMP, w);
    exit(-1);
  }
  if (nslit*w > radius*2) {
    nslit = floor(radius/w);
    printf("FermiChopper_ILL: %s: Too many slits to fit in the cylinder\n"
           "Adjusting nslit=%f\n", NAME_CURRENT_COMP, nslit);
  }
  if (length > radius*2) {
    length = sqrt(radius*radius - nslit*w*nslit*w/4);
    printf("FermiChopper_ILL: %s: Slit package is longer than the whole\n"
           "chopper cylinder. Adjusting length=%g [m]\n", NAME_CURRENT_COMP, length);
  }

  if (eff <= 0 || eff > 1) {
    eff = 0.95;
    printf("FermiChopper_ILL: %s: Efficiency is unrealistic\n"
           "Adjusting eff=%f\n", NAME_CURRENT_COMP, eff);
  }
  if (Qc <= 0) { Qc = 0.02176; m = 0; R0=0; }
  if (W <= 0) W=1e-6;
  
  if (verbose && nu)
    printf("FermiChopper_ILL: %s: frequency nu=%g [Hz] %g [rpm], time frame=%g [s] phase=%g [deg]\n"
      , NAME_CURRENT_COMP, nu, nu*60, 2/nu, -FCVars.t0*360*nu);
  
  /* fix for the wrong coordinate frame orientation to come back to McStas XYZ system */
  FCVars.omega *= -1;
%}

TRACE
%{

  /** local CALCULATION VARIABLES**************************************/

  /** Interaction with slitpacket ***************************/
  double slit_input;   /* length of the slits */
  double zr1,zr2;       /* distance to slitpacket entrance/exit in rotating frame */
  double xr1,xr2;       /* X entrance/exit position in rotating frame */

  /** Variables for calculating interaction with blades ***************/
  double m1,m2;    /* slope of the tangents */
  double b1,b2;    /* y-intersection of tangent */

  /**  Reflections ***********************************************/
  double t3a, t3b, distance_Wa, distance_Wb;
  double n1,n2,n3,n4;

  /** variables used for calculating new velocities after reflection **/
  double q;
  double vper, vpar;
  double arg;

  /**  Multiple Reflections  ******************************/
  int loopcounter=0;   /* How many reflections happen? */

  /** Time variables *********************************/
  double t3;      /* interaction time */
  double dt;      /* interaction intervals */
  double t1,t2;   /* cylinder intersection time */


/************** test, if the neutron interacts with the cylinder ***/
  if (cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight)) {

    if (t1 <= 0) {
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron started inside the cylinder, t1=%g (enter)\n", 
          NAME_CURRENT_COMP, t1);
      ABSORB;    /* Neutron started inside the cylinder */
    }

    dt=t2-t1;     /* total time of flight inside the cylinder  */
    PROP_DT(t1);  /* Propagates neutron to entrance of the cylinder */
    
    if (verbose > 2)
      printf("FermiChopper_ILL: %s: PROP_DT t1=%8.3g t2=%8.3g xyz=[%8.3g %8.3g %8.3g] v=[%8.3g %8.3g %8.3g] t=%8.3g (IN cyl).\n",
             NAME_CURRENT_COMP, t1, t2, x,y,z,vx,vy,vz,t);

    if(dt > fabs(0.5/FCVars.omega*2*PI) && verbose) {
      printf("FermiChopper_ILL: %s: Frequency too low. Method will fail.\n"
             "              Absorbing neutron\n", NAME_CURRENT_COMP);
      ABSORB;
    }

  /* Checks if neutron enters or leaves from top or bottom of cylinder. */
    if ( fabs(y) > yheight/2 ||
        fabs(y+vy*dt) > yheight/2 ) {
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (enter)\n", 
          NAME_CURRENT_COMP, y);
      ABSORB;
    }

  /*  checking wether the neutron can enter the chopper (slit channel) */
    xr1 = xstrich(x,z,t, FCVars.omega, FCVars.t0);
    if(fabs(xr1)>=nslit*w/2) {
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron X is outside cylinder aperture, xp1=%8.3g (enter)\n", 
          NAME_CURRENT_COMP, xr1);
      ABSORB;
    }

  /*********************** PROPAGATE TO SLIT PACKAGE **************************/


    /* Checking on which side of the Chopper the Neutron enters******/
    slit_input = 0.5*length;
    zr1 = zstrich(x,z,t, FCVars.omega, FCVars.t0);
    zr2 = zstrich(x+vx*dt,z+vz*dt,t+dt, FCVars.omega, FCVars.t0);
    if(zr1 < 0) slit_input *= -1;

    /*  Checking if the Neutron will hit the slits (Z) */
    zr1  -=  slit_input;
    zr2  -=  slit_input;

    if (zr2*zr1>0) {
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron Z does not change sign, zr1=%8.3g zr2=%8.3g (enter)\n", 
          NAME_CURRENT_COMP, zr1,zr2);
      ABSORB;
    }

    /* Calculating where/when Neutron hits the slits (Z) */
    t3 = zsecant(x,z,vx,vz,t,dt,slit_input, FCVars.omega, FCVars.t0);

    if((t3 < 0)||(t3 > dt)) {
      t3 = zinterpolation(x,z,vx,vz,t,dt,slit_input, FCVars.omega, FCVars.t0);
    }
    if((t3 < 0)||(t3 > dt)) {
      if (verbose) 
      	printf("FermiChopper_ILL: %s: Can not reach entrance of slits. dt=%g t3=%g\n", NAME_CURRENT_COMP, dt, t3);
      ABSORB;
    }

    /* Propagating whole system to that point */
    PROP_DT(t3);
    dt -= t3;
    SCATTER;
    
    if (verbose > 2)
      printf("FermiChopper_ILL: %s: PROP_DT t3=%8.3g dt=%8.3g xyz=[%8.3g %8.3g %8.3g] length=%g (slit enter).\n",
             NAME_CURRENT_COMP, t3, dt, x,y,z, slit_input);

    /* Checking if neutron hits the slits entrance window (X) */
    xr1 = xstrich(x,z,t, FCVars.omega, FCVars.t0);
    if(fabs(xr1) >= nslit*w/2) {
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron X is outside slit package, xp1=%8.3g (enter)\n", 
        NAME_CURRENT_COMP, xr1);
      ABSORB;
    }

    /* Calculating where/when Neutron leaves the slits (Z) */
    t3 = zsecant(x,z,vx,vz,t,dt,-slit_input, FCVars.omega, FCVars.t0);
    if((t3 < 0) || (t3 > dt)){
      t3 = zinterpolation(x,z,vx,vz,t,dt,-slit_input, FCVars.omega, FCVars.t0);
    }
    if((t3 <= 0) || (t3 > dt)){
      if (verbose) 
      	printf("FermiChopper_ILL: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g\n", NAME_CURRENT_COMP, dt, t3);
      ABSORB;
    } else dt=t3;

  /********************* PROPAGATION INSIDE THE SLIT PACKET *******************/

    /* Which slit was hit ? */
    n1 = floor(xr1/w);


  /******************* BEGIN LOOP FOR MULTIPLE REFLECTIONS ********************/

    for(loopcounter; loopcounter<=FCILL_MAXITER;loopcounter++){
    	double dt_to_tangent=0; /* time shift to go to tangent intersection */

  /* Calculate most probable time for interaction with blades by using tangents */
      m1 = xstrich(vx,vz,t, FCVars.omega, FCVars.t0)
         + FCVars.omega * zstrich(x,z,t, FCVars.omega, FCVars.t0);
      m2 = xstrich(vx,vz,t+dt, FCVars.omega, FCVars.t0)
         + FCVars.omega * zstrich(x+vx*dt,z+vz*dt,t+dt,FCVars.omega,FCVars.t0);

      b1 = xstrich(x,z,t, FCVars.omega, FCVars.t0) - m1*t;
      b2 = xstrich(x+vx*dt,z+vz*dt,t+dt, FCVars.omega, FCVars.t0) - m2*(t+dt);

      if (m1-m2) dt_to_tangent = ((b2-b1)/(m1-m2))-t;
      else       dt_to_tangent = -1;

      /* If method with tangents doesn't succeed, just take the middle of the interval */
      if((dt_to_tangent < 0)||(dt_to_tangent > dt)) dt_to_tangent=dt*0.5;

      /* Calculate different positions for the neutron to determine interaction. */

      /*...at the end of the slit: */
      n2 = floor(xstrich(x+(vx*dt),z+(vz*dt),t+dt, FCVars.omega, FCVars.t0)/w);

      /*...at the before calculated t3: tangent intersection point */
      n3 = floor(xstrich(x+(vx*dt_to_tangent),z+(vz*dt_to_tangent),t+dt_to_tangent, FCVars.omega, FCVars.t0)/w);

      if (verbose > 2)
        printf("FermiChopper_ILL: %s: t3=%8.3g n=[%g %g %g] (time at tangent intersection).\n",
             NAME_CURRENT_COMP, dt_to_tangent, n1, n2, n3);

      /* Does the neutron stay in the same slit ? */
      if((n2!=n1)||(n3!=n1)){

        /* Choosing the first time it isn't in the slit anymore */
        if(n3!=n1){
          n2 = n3;
        }

      /************ABSORB to save calculation time ******************/
        if (m == 0 || R0 == 0) {
          if (verbose > 2) 
          	printf("FermiChopper_ILL: %s: ABSORB Neutron hits absorbing coating (change slit).\n", 
              NAME_CURRENT_COMP);
          ABSORB;
        }


  /********************** WHEN DOES IT HIT THE BLADE? *************************/

        /*********** SECANT METHOD ****************************/
        /* get position of slit wall towards which neutron is propagating */
        if (n2 > n1) {     /* X' positive side of slit is in principle the first intersection to test*/
          distance_Wa = n1*w+w;
          distance_Wb = n1*w;
        } else {           /* X' negative side of slit */
        	distance_Wb = n1*w+w;
          distance_Wa = n1*w;
        }

        /* time shift to reach slit wall point in [0,dt_to_tangent]: X'=distance_W slit_wall */
        t3a = xsecant(x,z,vx,vz,t,dt,distance_Wa, FCVars.omega, FCVars.t0);
        t3b = xsecant(x,z,vx,vz,t,dt,distance_Wb, FCVars.omega, FCVars.t0);
        if (t3b < 0) t3 = t3a;
        else if (t3a < 0 && t3b > 0) t3 = t3b;
        else t3 = (t3a < t3b ? t3a : t3b);

        /***** INTERPOLATION USED WHEN SECANT METHOD FAILS ****/
        /* try second intersection method in case of failure */
        if ((t3 < 0) || (t3 > dt)) {
        	t3a = xinterpolation(x,z,vx,vz,t,dt,distance_Wa, FCVars.omega, FCVars.t0);
        	t3b = xinterpolation(x,z,vx,vz,t,dt,distance_Wb, FCVars.omega, FCVars.t0);
        	if (t3b < 0) t3 = t3a;
		      else if (t3a < 0 && t3b > 0) t3 = t3b;
		      else t3 = (t3a < t3b ? t3a : t3b);
        }

        /* Check for errors in calculation*******/
        if ((t3 < 0) || (t3 > dt)) {
          if (verbose) 
          	printf("FermiChopper_ILL: %s: Reflecting interpolation Problem. dt=%8.3g t3=%8.3g\n", 
          	NAME_CURRENT_COMP, dt, t3);
          ABSORB;
        }
        /* Propagate whole system to that point */
        PROP_DT(t3); dt -= t3;
        if (verbose > 2)
          printf("FermiChopper_ILL: %s: PROP_DT t3=%8.3g dt=%8.3g xyz=[%8.3g %8.3g %8.3g] (on wall).\n",
             NAME_CURRENT_COMP, t3, dt, x,y,z);

        /* Check if this point is inside the slit packet */
        if(fabs(zstrich(x,z,t, FCVars.omega, FCVars.t0)) > fabs(slit_input)){
          if (verbose > 2) 
          	printf("FermiChopper_ILL: %s: Neutron is outside slit pack (on slit wall).\n", 
                NAME_CURRENT_COMP);
          break;
        }

  /******************** REFLECTION ALGORITHM ********************************/
        vper    = xstrich(vx,vz,t, FCVars.omega, FCVars.t0); /* perpendicular velocity (to blade) */
        vpar    = zstrich(vx,vz,t, FCVars.omega, FCVars.t0);  /* parallel velocity (to blade) */
        q       = 2*MS2AA*(fabs(vper));

        if (q > Qc && W){
          arg = (q-m*Qc)/W;
          if (arg < 10.0) p *= 0.5*(1-tanh(arg))*(1-alpha*(q-Qc));
          else {
            if (verbose > 2) 
            	printf("FermiChopper_ILL: %s: ABSORB Neutron hits absorbing coating (on slit wall).\n", 
                NAME_CURRENT_COMP);
            ABSORB;
          }
        }

        if (R0 != 0.0){
          p *= R0;

          vper *= (-1);   /* Mirroring perpendicular velocity */

          /**************SET NEW VELOCITIES***********/
          vx =  vper*cos(FCVars.omega*(t-FCVars.t0))
             -  vpar*sin(FCVars.omega*(t-FCVars.t0));
          vz =  vper*sin(FCVars.omega*(t-FCVars.t0))
             +  vpar*cos(FCVars.omega*(t-FCVars.t0));
          SCATTER;
        } else {
          if (verbose > 2) 
          	printf("FermiChopper_ILL: %s: ABSORB Neutron hits absorbing coating (R0=0).\n", 
              NAME_CURRENT_COMP);
         ABSORB;
        }


        /* Recalculating when Neutron will leave the slitpacket */
        t3 = zsecant(x,z,vx,vz,t,dt,-slit_input,FCVars.omega,FCVars.t0);
        if((t3 < 0) || (t3 > dt)) {
          t3=zinterpolation(x,z,vx,vz,t,dt,-slit_input,
                             FCVars.omega,FCVars.t0);
        }
        /* Check for errors in calculation*******/
        if ((t3 < 0) || (t3 > dt)) {
          if (verbose) 
          	printf("FermiChopper_ILL: %s: Reflecting interpolation Problem. dt=%8.3g t3=%8.3g\n", 
          	NAME_CURRENT_COMP, dt, t3);
          ABSORB;
        } else  dt=t3;
      } /* end if n2 != n2 != n3 */
      else break;
    } /* end for */
  /********************* END OF THE FOR LOOP **********************************/

    /****New time of cylinder intersection will be calculated**********/
    if (!cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight)) {
    	if (verbose > 2) 
    		printf("FermiChopper_ILL: %s: ABSORB Neutron has unexpectidely exited cylinder ! (exiting)\n", 
    			NAME_CURRENT_COMP);
    	ABSORB; 
    }

    if (t1 > 0 && verbose) {
      printf("FermiChopper_ILL: %s: Neutrons are leaving chopper in the wrong direction! \n", NAME_CURRENT_COMP);
    }

    if (t2 <= 0 && verbose) {
      printf("FermiChopper_ILL: %s: Neutrons are leaving chopper without any control\n", NAME_CURRENT_COMP);
    }

  /*********** PROPAGATE TO CYLINDER SURFACE ***********************************/
    PROP_DT(t2);
    SCATTER;
    
    if (verbose > 2)
      printf("FermiChopper_ILL: %s: t1=%8.3g PROP_DT t2=%8.3g xyz=[%8.3g %8.3g %8.3g] (OUT cyl).\n",
             NAME_CURRENT_COMP, t1, t2, x,y,z);

    /*****Checking if the neutron left the cylinder by his top or bottom **/
    if  ( fabs(y) > yheight/2 ){
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (exiting)\n", 
          NAME_CURRENT_COMP, y);
      ABSORB;
    }


    /*****Checking if neutron hits chopper exit ***/
    if(fabs(xstrich(x,z,t,FCVars.omega,FCVars.t0))>=nslit*w/2){
      if (verbose > 2) 
      	printf("FermiChopper_ILL: %s: ABSORB Neutron X is outside slit package cylinder, xp1=%8.3g (exiting)\n", 
          NAME_CURRENT_COMP, xstrich(x,z,t,FCVars.omega,FCVars.t0));
      ABSORB;
    }

    /**** Transmission coefficent******/
    p = p*eff;          //finite cross section + transmission

  } /* end if cylinder_intersect */
  else {
    if (verbose > 2 && 0) 
    	printf("FermiChopper_ILL: %s: ABSORB Neutron has not interacted with FC\n", 
        NAME_CURRENT_COMP);
    ABSORB;
  }

/************************ TIME OF FLIGHT RESET ************************/
  if (zero_time && nu)
    t -= (((int)((t+1/(4*nu))/(1/(2*nu))))*(1/(2*nu)));
%}

MCDISPLAY
%{
  double index=0;
  double xpos, zpos;
  double ymax = yheight/2; 
  double ymin = -ymax;
  
  /* cylinder top/center/bottom  */
  circle("xz", 0,ymax,0,radius);
  circle("xz", 0,0   ,0,radius);
  circle("xz", 0,ymin,0,radius);
  /* vertical lines to make a kind of volume */
  line( 0  ,ymin,-radius, 0  ,ymax,-radius);
  line( 0  ,ymin, radius, 0  ,ymax, radius);
  line(-radius,ymin, 0  ,-radius,ymax, 0  );
  line( radius,ymin, 0  , radius,ymax, 0  );
  /* slit package */
  index = -nslit/2;
  zpos  = length/2;
  for (index = -nslit/2; index < nslit/2; index++) {
    xpos = index*w;
    multiline(5, xpos, ymin, -zpos,
                 xpos, ymax, -zpos,
                 xpos, ymax, +zpos,
                 xpos, ymin, +zpos,
                 xpos, ymin, -zpos);
  }
  /* cylinder inner sides containing slit package */
  xpos = nslit*w/2;
  zpos = sqrt(radius*radius - xpos*xpos);
  multiline(5,   xpos, ymin, -zpos,
                 xpos, ymax, -zpos,
                 xpos, ymax, +zpos,
                 xpos, ymin, +zpos,
                 xpos, ymin, -zpos);
  xpos *= -1;
  multiline(5,   xpos, ymin, -zpos,
                 xpos, ymax, -zpos,
                 xpos, ymax, +zpos,
                 xpos, ymin, +zpos,
                 xpos, ymin, -zpos);
%}
END