/*******************************************************************************
*
* 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