/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: FermiChopper.
*
* %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: EF, Nov 2005: completely rewrote comp.
* Modified by: EF, Oct 2008: fix chopper orientation
* Modified by: EF, Mar 2009: improved intersection algorithm
*
* 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. The slit package may be made curved and use
* super-mirror coating.
*
* Example:
* FermiChopper(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)
*
* %VALIDATION
* Apr 2005: extensive external test, most problems solved (cf. 'Bugs')
* Validated by: K. Lieutenant, E. Farhi
*
* limitations:
* no absorbing blade width used
*
* %Parameters
* INPUT PARAMETERS:
*
* Geometrical chopper constants:
* radius: [m]       chopper cylinder radius
* 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
* yheight: [m]      height of slit package
* nu: [Hz]          chopper frequency. Omega=2*PI*nu in rad/s, nu*60 in rpm. Positive value corresponds to a counter-clockwise rotation around y.
* eff: [1]          efficiency = transmission x fraction of transparent material
* verbose: [1]      set to 1,2 or 3 gives debugging information
* curvature: [m-1]  Curvature of slits (1/radius of curvature).
*
* 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, 2=auto adjust phase
* phase: [deg]      chopper phase at t=0
* delay: [s]        sets phase so that transmision is centered on 'delay'
*
* CALCULATED PARAMETERS:
* FCVars :  []     structure
*
* %L
* <a href="Vitess_ChopperFermi.html">Vitess_ChopperFermi</a> component by
* G. Zsigmond, imported from Vitess by K. Lieutenant.
*
* %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



SETTING PARAMETERS (phase=0, radius=0.04, nu=100,
w=0.00022475, nslit=200, 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, yheight=0.08,
curvature=0,delay=0)



SHARE
%{
%include "ref-lib"
#ifndef FermiChopper_TimeAccuracy
#define FermiChopper_TimeAccuracy 1e-9
#define FermiChopper_MAXITER      100

/* Definition of internal variable structure: all counters */
struct FermiChopper_struct
{
double omega;  /* chopper rotation */
double ph0;    /* chopper rotation */
double t0;     /* chopper rotation */
double C_slit;          /* slit curvature radius in [m] */
double L_slit;          /* slit package length [m] */
double sum_t;
double sum_v;
double sum_N;
double sum_N_pass;
/* events */
long absorb_alreadyinside;
long absorb_topbottom;
long absorb_cylentrance;
long absorb_sideentrance;
long absorb_notreachentrance;
long absorb_packentrance;
long absorb_slitcoating;
long warn_notreachslitwall;
long absorb_exitslitpack;
long absorb_maxiterations;
long absorb_wrongdirection;
long absorb_nocontrol;
long absorb_cylexit;
long warn_notreachslitoutput;
char compcurname[256];
};

/*****************************************************************************
* FC_zrot: returns Z' in rotating frame, from X,Z and t,omega,ph0
****************************************************************************/
#pragma acc routine seq
double FC_zrot(double X, double Z, double T, struct FermiChopper_struct FCs)
{
  double omega =FCs.omega;
  double ph0   =FCs.ph0;

  return( Z*cos(omega*T+ph0)-X*sin(omega*T+ph0) );
}

/*****************************************************************************
 * FC_xrot: returns X' in rotating frame, from X,Z and omega,t,ph0
 *          additional coordinate shift in case of curved slits
 ****************************************************************************/
#pragma acc routine seq
double FC_xrot(double X, double Z, double T, struct FermiChopper_struct FCs)
{
  double omega =FCs.omega;
  double ph0   =FCs.ph0;
  double C_slit=FCs.C_slit;
  double ret, tmp;

  ret = X*cos(omega*T+ph0)+Z*sin(omega*T+ph0);

  if (C_slit) {
    tmp  = fabs(FC_zrot(X, Z, T, FCs));
    if (tmp < FCs.L_slit/2) {
      tmp  = (FCs.L_slit/2 - tmp)*C_slit;
      ret += (1-sqrt(1-tmp*tmp))/C_slit;
    }
  }
  return( ret );
}

/*****************************************************************************
 * FC_xzrot_dt(x,z,vx,vz, t,dt, type='x' or 'z', FCs)
 *   returns X' or Z' in rotating frame, from X,Z and t,omega,ph0
 *              taking into account propagation with velocity during time dt
 ****************************************************************************/
#pragma acc routine seq
double FC_xzrot_dt(double x, double z, double vx, double vz,
                   double t, double dt, char type, struct FermiChopper_struct FCs)
{
  if (dt) /* with propagation */
    return( (type == 'x' ? FC_xrot(x+vx*dt, z+vz*dt, t+dt, FCs)
                         : FC_zrot(x+vx*dt, z+vz*dt, t+dt, FCs)) );
  else    /* without propagation */
    return( (type == 'x' ? FC_xrot(x,z,t,FCs)
                         : FC_zrot(x,z,t,FCs)) );
}

/*****************************************************************************
 * FC_xzbrent(x,z,vx,vz, t,dt, type='x' or 'z', d, FCs)
 *   solves X'=d and Z'=d with Brent algorithm in time interval [0, dt].
 *           Returns time within [0,dt], from NumRecip in C, chap 9, p360 (zbrent)
 *           ERRORS: return -1 not used
 *                          -2 if exceed MAX iteration
 *                          -3 no sign change in range
 ****************************************************************************/
#pragma acc routine seq
double FC_xzbrent(double x, double z, double vx, double vz,
                  double t, double dt,
                  char type, double d, struct FermiChopper_struct FCs)
{
  int iter;
  double a=0,b=dt;
  double c,dd,e,min1,min2;
  double tol=FermiChopper_TimeAccuracy;
  double EPS=FermiChopper_TimeAccuracy;
  double fa=FC_xzrot_dt(x,z,vx,vz, t,a, type, FCs) - d;
  double fb=FC_xzrot_dt(x,z,vx,vz, t,b, type, FCs) - d;
  double fc,p,q,r,s,tol1,xm;

  if (fb*fa > 0.0) return -3;
  fc=fb;
  for (iter=1;iter<=FermiChopper_MAXITER;iter++) {
    if (fb*fc > 0.0) {
      c=a;
      fc=fa;
      e=dd=b-a;
    }
    if (fabs(fc) < fabs(fb)) {
      a=b;
      b=c;
      c=a;
      fa=fb;
      fb=fc;
      fc=fa;
    }
    tol1=2.0*EPS*fabs(b)+0.5*tol;
    xm=0.5*(c-b);
    if (fabs(xm) <= tol1 || fb == 0.0) return b;
    if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
      s=fb/fa;
      if (a == c) {
        p=2.0*xm*s;
        q=1.0-s;
      } else {
        q=fa/fc;
        r=fb/fc;
        p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
        q=(q-1.0)*(r-1.0)*(s-1.0);
      }
      if (p > 0.0)  q = -q;
      p=fabs(p);
      min1=3.0*xm*q-fabs(tol1*q);
      min2=fabs(e*q);
      if (2.0*p < (min1 < min2 ? min1 : min2)) {
        e=dd;
        dd=p/q;
      } else {
        dd=xm;
        e=dd;
      }
    } else {
      dd=xm;
      e=dd;
    }
    a=b;
    fa=fb;
    if (fabs(dd) > tol1)
      b += dd;
    else
      b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
    fb=FC_xzrot_dt(x,z,vx,vz, t,b, type, FCs) - d;
  }
  return -2;
} /* FC_xzbrent */

/*****************************************************************************
 * Wrappers to intersection algorithms
 ****************************************************************************/
#pragma acc routine seq
double FC_xintersect(double x, double z, double vx, double vz,
                   double t, double dt,
                   double d, struct FermiChopper_struct FCs)
{
  return(FC_xzbrent(x, z, vx, vz, t, dt, 'x', d, FCs));
}
#pragma acc routine seq
double FC_zintersect(double x, double z, double vx, double vz,
                   double t, double dt,
                   double d, struct FermiChopper_struct FCs)
{
  return(FC_xzbrent(x, z, vx, vz, t, dt, 'z', d, FCs));
}

#endif
%}

DECLARE
%{
  struct FermiChopper_struct FCVars;
%}

INITIALIZE
%{

/************************ CALCULATION CONSTANTS *****************************/
  strcpy(FCVars.compcurname, NAME_CURRENT_COMP);

  FCVars.omega    = 2*PI*nu;
  if (!phase && delay) {
     FCVars.ph0= fmod(-delay*nu*360,360)*DEG2RAD;
  } else FCVars.ph0      = phase*DEG2RAD;
  FCVars.sum_t=FCVars.sum_v=FCVars.sum_N=FCVars.sum_N_pass=0;

  /* check of input parameters */
  if (nslit < 1) nslit=1;
  if (yheight <= 0) exit(printf("FermiChopper: %s: FATAL: unrealistic cylinder yheight =%g [m]\n", NAME_CURRENT_COMP, yheight));

  if (m <= 0) { m=0; R0=0; }
  if (radius <= 0) {
    printf("FermiChopper: %s: FATAL: Unrealistic cylinder radius radius=%g [m]\n", NAME_CURRENT_COMP, radius);
    exit(-1);
  }
  if (xwidth > 0 && xwidth < radius*2 && nslit > 0) {
    w = xwidth/nslit;
  }
  if (w <= 0) {
    printf("FermiChopper: %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: %s: Too many slits to fit in the cylinder\n"
           "    Adjusting nslit=%f\n", NAME_CURRENT_COMP, nslit);
  }
  if (length > radius*2) {
    length = 2*sqrt(radius*radius - nslit*w*nslit*w/4);
    printf("FermiChopper: %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: %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 (curvature) {
    FCVars.C_slit = curvature;
    if (1 < fabs(radius*curvature))
      exit(printf("FermiChopper: %s: Slit curvature is unrealistic\n",
           NAME_CURRENT_COMP));
  }
  FCVars.L_slit = length;
  if (verbose && nu)
    printf("FermiChopper: %s: Frequency nu=%g [Hz] %g [rpm], time frame=%g [s] phase=%g [deg]\n"
      , NAME_CURRENT_COMP, nu, nu*60, 2/nu, FCVars.ph0*RAD2DEG);

  FCVars.absorb_alreadyinside    = 0;
  FCVars.absorb_topbottom        = 0;
  FCVars.absorb_cylentrance      = 0;
  FCVars.absorb_sideentrance     = 0;
  FCVars.absorb_notreachentrance = 0;
  FCVars.absorb_packentrance     = 0;
  FCVars.absorb_slitcoating      = 0;
  FCVars.warn_notreachslitwall   = 0;
  FCVars.absorb_exitslitpack     = 0;
  FCVars.absorb_maxiterations    = 0;
  FCVars.absorb_wrongdirection   = 0;
  FCVars.absorb_nocontrol        = 0;
  FCVars.absorb_cylexit          = 0;
  FCVars.warn_notreachslitoutput = 0;

  /* fix for the wrong coordinate frame orientation to come back to McStas XYZ system */
  FCVars.omega *= -1;
  FCVars.ph0   *= -1;
  FCVars.t0     = -FCVars.ph0/FCVars.omega;
%}

TRACE
%{

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

  /** Interaction with slit package ***************************/
  double slit_input; /* length of the slits */

  /** Variables for calculating interaction with blades ***************/
  double xp1, zp1, vxp1;

  /**  Reflections ***********************************************/
  double n1;

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

  /** Time variables *********************************/
  double t3=0;      /* interaction time at n3 position (slit wall) */
  double t1=0,t2=0; /* cylinder intersection time (at entry and exit of slit pack - n1 n2 - or cylinder) */
  double dt;        /* interaction intervals (e.g. time for exiting slit pack) */

  double X[5],Z[5]; /* position of interaction locations in the Fermi chopper rotating frame
		       [0]=cylinder input
		       [1]=slit input
		       [2]=slit wall reflection
		       [3]=slit exit
		       [4]=cylinder exit
		    */
  double ref = 0;
  double par[5] = {R0, Qc, alpha, m, W};
  
  /************************ TIME OF FLIGHT RESET ************************/
  if (zero_time == 1 && nu)
    t -= floor( (t+1/(4*nu))*(2*nu) )/(2*nu);

  /* zero arrays used to store positions */
  for (loopcounter=0; loopcounter < 5; loopcounter++)
    X[loopcounter] = Z[loopcounter]=0;

  /************** test, if the neutron interacts with the cylinder ***/
  if (cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight))
  {
    if (t1 <= 0) { /* Neutron started inside the cylinder */
      if (verbose > 0 && FCVars.absorb_alreadyinside<FermiChopper_MAXITER) printf("FermiChopper: %s: ABSORB Neutron started inside the cylinder, t1=%8.3g (enter).\n",
        NAME_CURRENT_COMP, t1);
      #pragma acc atomic
      FCVars.absorb_alreadyinside = FCVars.absorb_alreadyinside + 1;
      ABSORB;
    }
    if (verbose > 2)
      printf("FermiChopper: %s:         t1=%8.3g t2=%8.3g xyz=[%8.3g %8.3g %8.3g] v=[%8.3g %8.3g %8.3g] t=%8.3g (init).\n",
           NAME_CURRENT_COMP, t1, t2, x,y,z,vx,vy,vz,t);

    dt=t2-t1;     /* total time of flight inside the cylinder  */
    PROP_DT(t1);  /* Propagates neutron to entrance of the cylinder */
    SCATTER;

    if (verbose > 2)
      printf("FermiChopper: %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);

    /* neutron must not enter or leave from top or bottom of cylinder. */
    if (fabs(y) >= yheight/2.0 || fabs(y+vy*dt) >= yheight/2.0) {
      if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (enter).\n",
        NAME_CURRENT_COMP, y);
      #pragma acc atomic
      FCVars.absorb_topbottom =  FCVars.absorb_topbottom + 1;
      ABSORB;
    }

    vxp1 = sqrt(vx*vx+vy*vy+vz*vz);
    double tmp=p*vxp1;
    #pragma acc atomic
    FCVars.sum_v = FCVars.sum_v + tmp;
    tmp = p*t;
    #pragma acc atomic
    FCVars.sum_t = FCVars.sum_t + tmp;
    #pragma acc atomic
    FCVars.sum_N = FCVars.sum_N + p;

    if (zero_time > 1 && FCVars.sum_N) { /* automatic phase adjustment */
      double mean_t, mean_phase;
      mean_t     = FCVars.sum_t/FCVars.sum_N;
      mean_phase = fmod(mean_t*nu*2*PI, 2*PI);
      /* now we shift the phase so that the neutron pulse is centered on the slit pack */
      mean_phase+= radius/vxp1*2*PI*nu;
      #pragma acc atomic write
      FCVars.ph0 = mean_phase;
    }


    /* neutron must enter the cylinder opening: |X'| < full package width*/
    xp1 = FC_xrot(x,z, t,FCVars); /* X'(t) */
    if (fabs(xp1) >= nslit*w/2) {
      if (verbose > 2)
        printf("FermiChopper: %s: ABSORB Neutron X is outside cylinder aperture, x'=%8.3g > %g (enter).\n",
        NAME_CURRENT_COMP, xp1, nslit*w/2);
      #pragma acc atomic
      FCVars.absorb_cylentrance = FCVars.absorb_cylentrance + 1;
      ABSORB;
    }

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

    /* zp1 = Z' at entrance of cylinder Z'(t) */
    zp1  = FC_zrot(x,z, t, FCVars);

    X[0] = xp1; Z[0] = zp1;

    /* here we should have sqrt(x*x+z*z) == sqrt(xp1*xp1+zp1*zp1) */
    t3 = sqrt(x*x+z*z) / sqrt(xp1*xp1+zp1*zp1);
    if (t3 < 0.99 || 1.01 < t3) {
      if (verbose > 1 && FCVars.absorb_cylentrance < FermiChopper_MAXITER)
        printf("FermiChopper: %s: ABSORB Neutron radius on cylinder in static frame r=%g does not match that of rotating frame r'=%g.\n",
        NAME_CURRENT_COMP, sqrt(x*x+z*z), sqrt(xp1*xp1+zp1*zp1));
      #pragma acc atomic
      FCVars.absorb_cylentrance = FCVars.absorb_cylentrance +1 ;
      ABSORB;
    }

    /* Checking on which side of the Chopper the Neutron enters: sign(Z') */
    slit_input = (zp1 > 0 ? length/2 : -length/2);

    /* time shift to reach slit package in [0,time to exit cylinder]: Z'=slit_input */
    /* t3 is used here as a tmp variable, will be redefined in for loop  */
    t3 = FC_zintersect(x,z,vx,vz, t,dt, slit_input, FCVars);

    if( (t3 < 0)||(t3 > dt) ) {
      if (verbose > 2 && FCVars.absorb_notreachentrance < FermiChopper_MAXITER) {
        printf("FermiChopper: %s: Can not reach entrance of slits. dt=%8.3g t3=%8.3g (intersection:1).\n",
        NAME_CURRENT_COMP, dt, t3);
        if (t3 == -3)
            printf("          No sign change to determine intersection\n");
        else if (t3 == -2)
            printf("          Max iterations reached\n");
        else if (t3 < 0)
            printf("          Error when solving intersection\n");
      }
      #pragma acc atomic 
      FCVars.absorb_notreachentrance = FCVars.absorb_notreachentrance + 1;
      ABSORB; /* neutron can not reach slit entrance */
    }

    /* Propagating to the slit package entrance */
    PROP_DT(t3); /* dt = t2-t1: time in cylinder */
    dt -= t3; /* remaining time from slit pack entry to exit of cylinder */
    xp1 = FC_xrot(x,z, t, FCVars); /* should be slit_input */
    zp1 = FC_zrot(x,z, t, FCVars);
    X[1] = xp1; Z[1] = zp1;

#ifndef OPENACC
    if (mcdotrace) {
      /* indicate position of neutron in mcdisplay */
      double xp2 = x; double zp2 = z; x = xp1; z=zp1; SCATTER; x=xp2; z=zp2;
    } else 
#endif
SCATTER;

    if (verbose > 2)
      printf("FermiChopper: %s: PROP_DT t=%8.3g dt=%8.3g x'=%8.3g z'=%8.3g length=%g (slit enter).\n",
           NAME_CURRENT_COMP, t, dt, xp1, zp1, slit_input);

    /* must have X'< slit package width at package Z'=slit_input */
    if (fabs(xp1) >= nslit*w/2) {
      if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package, x'=%8.3g > %g (enter).\n",
        NAME_CURRENT_COMP, xp1, nslit*w/2);
      #pragma acc atomic
      FCVars.absorb_packentrance =  FCVars.absorb_packentrance + 1;
      ABSORB;
    }

    /* solve Z'=-slit_input for time of exit of slit package */
    /* t3 is used here as a tmp variable, will be redefined in for loop  */
    t3 = FC_zintersect(x,z,vx,vz, t,dt*1.1, -slit_input, FCVars);

    if((t3 < FermiChopper_TimeAccuracy)||(t3 > dt)) {
      if (verbose > 1 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) {
        printf("FermiChopper: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g (intersection:2).\n",
          NAME_CURRENT_COMP, dt, t3);
        if (t3 == -3)
            printf("          No sign change to determine intersection\n");
        else if (t3 == -2)
            printf("          Max iterations reached\n");
        else if (t3 < 0)
            printf("          Error when solving intersection\n");
      }
      #pragma acc atomic
      FCVars.warn_notreachslitoutput =  FCVars.warn_notreachslitoutput + 1;
      /* we estimate analytically the time to the slit exit */
      t3 = length/vxp1;
    }
    dt = t3; /* reduce time interval to [0, time of slit exit] */

    /* here we should have dt*v = length (slit entrance -> exit) */
    t3 = fabs(FC_xzrot_dt(x,z,vx,vz,    t,   0 , 'z', FCVars) - FC_xzrot_dt(x,z,vx,vz,    t,   dt, 'z', FCVars))/length;
    if (fabs(t3-1) > 0.02) {
      if (verbose > 0 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER)
        printf("FermiChopper: %s: ABSORB Neutron propagation time v*dt/length=%g in slit does not match its length=%g (slit exit expected).\n",
          NAME_CURRENT_COMP, t3, length);
      #pragma acc atomic
      FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1;
      ABSORB;
    }

    /* dt= time shift to go to exit of slit package (or exit of cylinder in case of error) */
    /*
      |---------|
      | /       |
      |o        * (dt)
      |         |
      |---------|
     */

/*********************PROPAGATION INSIDE THE SLIT PACKAGE *******************/

    /* index of slit hit at entrance n1=-N/2 (xp1=-) ... N/2 (xp1=+) */
    n1 = floor(xp1/w);

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

    for (loopcounter=0; loopcounter<=FermiChopper_MAXITER; loopcounter++) {
      double dt_to_tangent=0; /* time shift to go to tangent intersection */
      double n2,n3;           /* slit indices */
      double xp2, zp2,xp3=0,vxp2=0,vzp1=0; /* position, velocity */
      double q;                   /* used for calculating new velocities after reflection */
      int    i;

      /* compute trajectory tangents: m1=Vz'+w.X'(t), m2=Vz'+w.X'(t+dt) */
      xp1 = FC_xrot    (x,z,          t,   FCVars);          /* X'(t)    current position */
      xp2 = FC_xzrot_dt(x,z,vx,vz,    t,   dt, 'x', FCVars); /* X'(t+dt) slit exit */
      zp2 = FC_xzrot_dt(x,z,vx,vz,    t,   dt, 'z', FCVars); /* Z'(t+dt) slit exit */

      /* slit index at the end of the slit: */
      n2 = floor(xp2/w);

      /* quick exit for absorbing walls when neutron changes slit index */
      if (n2 != n1 && (m <= 0 || R0 <= 0 || Qc <= 0)) {
        if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (change slit).\n",
            NAME_CURRENT_COMP);
	#pragma acc atomic
        FCVars.absorb_slitcoating =  FCVars.absorb_slitcoating + 1;
	ABSORB;
      }

      /* compute transversal velocity to determine their intersection */
      vxp1= FC_xrot    (vx+z*FCVars.omega,vz-x*FCVars.omega,
                                      t,   FCVars);          /* dX'(t)/dt slope at current position*/

      vxp2= FC_xrot    (vx+(z+vz*dt)*FCVars.omega,vz-(x+vx*dt)*FCVars.omega,
                                      t+dt,FCVars);          /* dX'(t+dt)/dt slope at slit exit */

      /* absolute time at tangent intersection, changed to time shift below */
      dt_to_tangent = (vxp1 - vxp2 ? (xp2 - xp1 - dt*vxp2)/(vxp1 - vxp2) : -1);

      /* If algorithm fails, take the middle of the interval*/
      if (dt_to_tangent < 0 || dt_to_tangent > dt) {
        if (verbose > 1)
          printf("FermiChopper: %s: WARNING could not determine tangent intersection dt=%g. Using middle.\n",
          NAME_CURRENT_COMP, dt_to_tangent);
        dt_to_tangent=dt*0.5;
      }

      /*
           *(dt_to_tangent, xp3)
      |---------|
      | /     \ |
  xp1 |o       \|(dt) xp2
      |         |
      |---------|
     */

      /* point coordinates at tangent intersection/middle point (max deviation from optical axis) */
      xp3 = FC_xzrot_dt(x,z,vx,vz, t, dt_to_tangent, 'x', FCVars); /* X'(t+dt_to_tangent) */

      /* slit index at the tangent intersection/middle point */
      n3 = floor(xp3/w);

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

      /* change slit means there is a reflection/intersection inside */
      if ( n2!=n1 || n3!=n1 ) {

        double t3a, t3b, distance_Wa, distance_Wb;
        if (m <= 0 || R0 <= 0 || Qc <= 0) {
          if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (change slit).\n",
            NAME_CURRENT_COMP);
	  #pragma acc atomic
          FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
	  ABSORB;
        }

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

        /* 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 */
        for (i=0; i< 2; i++) {
          /* first  attempt: [0,dt_to_tangent] (to max deviation location)
             second attempt: [0, dt]           (to slit exit) */

          double dt_search = (i == 0 ? dt_to_tangent : dt);
          t3a = FC_xintersect(x,z,vx,vz, t,dt_to_tangent, distance_Wa, FCVars);
          t3b = FC_xintersect(x,z,vx,vz, t,dt_to_tangent, distance_Wb, FCVars);
          if      (t3b < 0)             t3 = t3a;
          else if (t3a < 0 && t3b >= 0) t3 = t3b;
          else                          t3 = (t3a < t3b ? t3a : t3b);
          if (FermiChopper_TimeAccuracy < t3 && t3 < dt_search)
            break; /* we found the intersection location */
        }
        /* handle case where intersection search fails */
        if (t3 < FermiChopper_TimeAccuracy || t3 >= dt_to_tangent) {
          if (verbose > 1 && FCVars.warn_notreachslitwall < FermiChopper_MAXITER) {
            printf("FermiChopper: %s: Can not reach slit wall (iteration %i). dt=%8.3g t3=%8.3g (intersection:3).\n",
              NAME_CURRENT_COMP, loopcounter, dt_to_tangent, t3);
            if (t3 == -3)
              printf("        No sign change to determine intersection\n");
            else if (t3 == -2)
              printf("        Max iterations reached\n");
            else if (t3 < 0 || t3 >= dt_to_tangent)
              printf("        Error when solving intersection\n");
          }
          /* neutron can not reach slit wall. */
	  #pragma acc atomic
          FCVars.warn_notreachslitwall = FCVars.warn_notreachslitwall + 1;
          ABSORB;
        }

        /* Propagate to slit wall point (t+t3) on slit n3 wall */
        PROP_DT(t3); dt -= t3; /* dt: time remaining to slit exit after propagation */
        xp1 = FC_xrot(x,z, t, FCVars); /* X'(t+t3) : on slit wall */
        zp1 = FC_zrot(x,z, t, FCVars); /* Z'(t+t3) : on slit wall */
        X[2] = xp1; Z[2] = zp1;

        if (verbose > 2)
          printf("FermiChopper: %s: PROP_DT t3=%8.3g dt=%8.3g xyz=[%8.3g %8.3g %8.3g] (on wall). z'=%g\n",
           NAME_CURRENT_COMP, t3, dt, x,y,z, zp1);

        /* check if intersection point is still in the slit package, else exit loop */
        if (fabs(zp1) > length/2 || dt <= 0) {
          if (verbose > 2)
            printf("FermiChopper: %s: Neutron is outside slit pack (on slit wall).\n",
                NAME_CURRENT_COMP);
          break;
        }

    /*    here
      |---o-----|
      | /   \   |
      |/     \  |
      |       \ |(dt)
      |---------|
     */

        /* get velocity in rotating frame, on slit wall */
        vxp1 = FC_xrot(vx,vz, t, FCVars);
        vzp1 = FC_zrot(vx,vz, t, FCVars);

        q    = 2*V2Q*(fabs(vxp1));

        {
          /* double ref = 0;
	     double par[] = {R0, Qc, alpha, m, W};*/
          StdReflecFunc(q, par, &ref);
          if (ref>0) p *= ref;
          else {
            if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (on slit wall).\n",
              NAME_CURRENT_COMP);
            #pragma acc atomic
            FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
	    ABSORB;
          } /* Cutoff ~ 1E-10 */
        }
#ifndef OPENACC
        if (mcdotrace) {
          double xp2 = x; double zp2 = z;
          /* indicate position of neutron in mcdisplay */
          x = FC_xrot(x,z, t,FCVars); z= FC_zrot(x,z, t,FCVars); SCATTER; x=xp2; z=zp2;
        } else 
#endif
	SCATTER;

        /* reflect perpendicular velocity and compute new velocity in static frame */
        vxp1 *= -1;
        /* apply transposed Transformation matrix */
        vx = FC_xrot( vxp1,-vzp1, t,FCVars);
        vz = FC_zrot(-vxp1, vzp1, t,FCVars);

        /* recompute time to slit exit */
        /* solve Z'=-slit_input for time of exit of slit package */
        t3 = FC_zintersect(x,z,vx,vz, t,dt, -slit_input, FCVars);

        if(t3 < 0 || t3 > dt) {
          if (verbose > 1 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) {
            printf("FermiChopper: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g (intersection:4).\n",
              NAME_CURRENT_COMP, dt, t3);
            if (t3 == -3)
              printf("              No sign change to determine intersection\n");
            else if (t3 == -2)
              printf("              Max iterations reached\n");
            else if (t3 < 0)
              printf("              Error when solving intersection\n");
          }
	  #pragma acc atomic
          FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1;
          ABSORB;
          /* neutron can not reach slit output. */
        } else dt = t3; /* reduce time interval to [0, time of slit exit] */

      } /* end if changed slit index */
      else { /* neutron remains in the same slit: direct flight */
        if (dt > 0) PROP_DT(dt); /* go to slit exit */
        break;
      }
    } /* end for */

    xp1 = FC_xrot(x,z, t,FCVars);
    zp1 = FC_zrot(x,z, t,FCVars);
    X[3] = xp1; Z[3] = zp1; /* slit exit */

    if (fabs(xp1) >= nslit*w/2)
    {
      if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package, x'=%8.3g (slit exit).\n",
        NAME_CURRENT_COMP, xp1);
      #pragma acc atomic
      FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
      ABSORB;
    }

    if (loopcounter >= FermiChopper_MAXITER)
    {
      if (verbose > 1 && FCVars.absorb_maxiterations < FermiChopper_MAXITER)
      printf("FermiChopper: %s: Max iterations %i reached inside slit. Absorb.\n",
        NAME_CURRENT_COMP, FermiChopper_MAXITER);
      #pragma acc atomic
      FCVars.absorb_maxiterations =  FCVars.absorb_maxiterations + 1;
      ABSORB;
    }


    /********************* EXIT SLIT PACKAGE ********************************/

    /* propagation times to cylinder. Should use t2 to exit */
    if (!cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight)) {
      /* must be inside the cylinder */
      if (verbose > 1) printf("FermiChopper: %s: ABSORB Neutron has unexpectidely exited cylinder ! (exiting)\n",
        NAME_CURRENT_COMP);
      #pragma acc atomic
      FCVars.absorb_exitslitpack = FCVars.absorb_exitslitpack + 1;
      ABSORB;
    }

    if (t1 > 0)
    {
      if (verbose > 1 && FCVars.absorb_wrongdirection < FermiChopper_MAXITER)
      printf("FermiChopper: %s: Neutrons are leaving chopper\n"
             "              in the wrong direction. Absorb.\n", NAME_CURRENT_COMP);
      #pragma acc atomic
      FCVars.absorb_wrongdirection = FCVars.absorb_wrongdirection + 1;
      ABSORB;
    }

    if (t2 <= 0 && FCVars.absorb_nocontrol < FermiChopper_MAXITER)
    {
      if (verbose > 1)
      printf("FermiChopper: %s: Neutrons are leaving chopper\n"
             "              without any control. Absorb.\n", NAME_CURRENT_COMP);
      #pragma acc atomic
      FCVars.absorb_nocontrol = FCVars.absorb_nocontrol + 1;
      ABSORB;
    }

    /* propagate to cylinder surface (exit) */
    PROP_DT(t2);
    SCATTER;

    xp1 = FC_xrot(x,z, t,FCVars);
    zp1 = FC_zrot(x,z, t,FCVars);
    X[4] = xp1; Z[4] = zp1;

    if (verbose > 2)
      printf("FermiChopper: %s: t1=%8.3g PROP_DT t2=%8.3g xyz=[%8.3g %8.3g %8.3g] (OUT cyl). z'=%g\n",
           NAME_CURRENT_COMP, t1, t2, x,y,z, zp1);

    /* Check if the neutron left the cylinder by its top or bottom */
    if (fabs(y) >= yheight/2)
    {
      if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (exiting)\n",
        NAME_CURRENT_COMP, y);
      #pragma acc atomic
      FCVars.absorb_topbottom = FCVars.absorb_topbottom + 1;
      ABSORB;
    }

    /* must have X'< slit package width at package Z'=cylinder output */
    if (fabs(xp1) >= nslit*w/2)
    {
      if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package cylinder, xp1=%8.3g (exiting).\n",
        NAME_CURRENT_COMP, xp1);
      #pragma acc atomic
      FCVars.absorb_cylexit = FCVars.absorb_cylexit + 1;
      ABSORB;
    }

    /* Transmission coefficent */
    p = p*eff;          //finite cross section + transmission
    #pragma acc atomic
    FCVars.sum_N_pass = FCVars.sum_N_pass + p;

  } /* end if cylinder_intersect */

%}

SAVE
%{
  double mean_k, mean_v, mean_t, mean_w=0, mean_L=0;

  if (FCVars.sum_N) {
    mean_v = FCVars.sum_v/FCVars.sum_N;
    mean_t = FCVars.sum_t/FCVars.sum_N;
    mean_k = V2K*mean_v;
    if (mean_k) mean_L = 2*PI/mean_k;
    mean_w = VS2E*mean_v*mean_v;
    /* base opening time */
    double div, mean_phase;
    div        = atan2(w,length)/PI*180;
    mean_phase = fmod(mean_t*nu*360, 360);
    mean_phase+=radius/mean_v*360*nu;
    if (mean_phase > 180) mean_phase -= 360;

    if (!FCVars.sum_N_pass)
      printf("FermiChopper: %s: No neutron can pass the chopper.\n", NAME_CURRENT_COMP);
    if (!FCVars.sum_N_pass || verbose)
      printf("FermiChopper: %s\n"
           "              Mean velocity     v     = %g [m/s]\n"
           "              Mean wavelength   lambda= %g [Angs]\n"
           "              Mean energy       omega = %g [meV]\n"
           "              Mean arrival time t     = %g [s]\n"
           "              Mean phase              = %g [deg] (%s)\n"
           "              Slit pack divergence    = %g [deg] (full width)\n"
           "              Opening time      dt    = %g [s]\n"
           "              Intensity reaching FC   = %g [n/s]\n"
           "              Intensity passing  FC   = %g [n/s]\n"
           , NAME_CURRENT_COMP,
           mean_v, mean_L, mean_w, mean_t, mean_phase,
           (zero_time > 1 ? "set automatically" : "use phase=-(this value) to optimize"),
           2*div,
           (nu ? fabs(div/PI/nu) : 1),
           FCVars.sum_N,
           FCVars.sum_N_pass);
    if (!FCVars.sum_N_pass || verbose) {
      printf("FermiChopper: %s: Lost events anaylsis\n"
             "              Already inside:            %li\n"
             "              By Top/Bottom of cylinder: %li\n"
             "              At cylinder entrance:      %li\n"
             "              Hit cyl. entrance sides:   %li\n"
             "              Can't prop. to slit pack:  %li\n"
             "              At slit pack entrance:     %li\n"
             "              On absorbing slit coating: %li\n"
             "              Exiting slit pack:         %li\n"
             "              Too many iterations:       %li\n"
             "              Prop. in wrong direction : %li\n"
             "              Mad neutron (no control):  %li\n"
             "              At cylinder exit:          %li\n"
      , NAME_CURRENT_COMP,
      FCVars.absorb_alreadyinside,
      FCVars.absorb_topbottom,
      FCVars.absorb_cylentrance,
      FCVars.absorb_sideentrance,
      FCVars.absorb_notreachentrance,
      FCVars.absorb_packentrance,
      FCVars.absorb_slitcoating,

      FCVars.absorb_exitslitpack,
      FCVars.absorb_maxiterations,
      FCVars.absorb_wrongdirection,
      FCVars.absorb_nocontrol,
      FCVars.absorb_cylexit);

      if (FCVars.warn_notreachslitwall || FCVars.warn_notreachslitoutput)
        printf("Warning:      Can not reach slit wall:   %li\n"
               "Warning:      Can not reach slit output: %li\n",
        FCVars.warn_notreachslitwall,
        FCVars.warn_notreachslitoutput);
    }

  } else {
    printf("FermiChopper: %s: No neutron can reach the chopper.\n", NAME_CURRENT_COMP);
  }
%}

MCDISPLAY
%{
  double index_x=0;
  double index_z=0;
  double xpos, zpos;
  double Nz,Nx;
  double ymin=-yheight/2.0;
  double ymax= yheight/2.0;

  double omega=FCVars.omega;
  FCVars.omega=0;
  // FCVars.ph0  =0;
  Nz = (FCVars.C_slit ?  4 : 1);
  Nx = (nslit > 11    ? 11 : nslit);
  FCVars.C_slit *= -1;
  
  /* 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 */
  for (index_x = -Nx/2; index_x < Nx/2; index_x++) {
    for (index_z = -Nz/2; index_z < Nz/2; index_z++) {
      double xs1, zs1, zs2;
      double xp1, xp2, zp1, zp2;
      zs1 = length*index_z/Nz;
      zs2 = length*(index_z+1)/Nz;
      xs1 = w*nslit*index_x/Nx;
      xp1 = FC_xrot(xs1, zs1, 0, FCVars);
      xp2 = FC_xrot(xs1, zs2, 0, FCVars);
      zp1 = FC_zrot(xs1, zs1, 0, FCVars);
      zp2 = FC_zrot(xs1, zs2, 0, FCVars);
      multiline(5, xp1, ymin, zp1,
                   xp1, ymax, zp1,
                   xp2, ymax, zp2,
                   xp2, ymin, zp2,
                   xp1, ymin, zp1);
    }
  }
  /* cylinder inner sides containing slit package */
  double xp1, xp2, zp1, zp2;
  xpos = nslit*w/2;
  zpos = sqrt(radius*radius - xpos*xpos);
  xp1 = FC_xrot(xpos, -zpos, 0, FCVars);
  xp2 = FC_xrot(xpos, +zpos, 0, FCVars);
  zp1 = FC_zrot(xpos, -zpos, 0, FCVars);
  zp2 = FC_zrot(xpos, +zpos, 0, FCVars);
  multiline(5,  xp1, ymin, zp1,
                xp1, ymax, zp1,
                xp2, ymax, zp2,
                xp2, ymin, zp2,
                xp1, ymin, zp1);
  xpos *= -1;
  xp1 = FC_xrot(xpos, -zpos, 0, FCVars);
  xp2 = FC_xrot(xpos, +zpos, 0, FCVars);
  zp1 = FC_zrot(xpos, -zpos, 0, FCVars);
  zp2 = FC_zrot(xpos, +zpos, 0, FCVars);
  multiline(5,  xp1, ymin, zp1,
                xp1, ymax, zp1,
                xp2, ymax, zp2,
                xp2, ymin, zp2,
                xp1, ymin, zp1);
  FCVars.omega=omega;
%}
END