/****************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Pol_bender_tapering
*
* %I
* Written by: Erik Bergbaeck Knudsen (erkn@fysik.dtu.dk)
* Date: June 2012
* Origin: DTU Physics
*
* Polarising bender.
*
* %D
*
* A component modelling a set of identical blades bent to form a multcihannel (nslits)
* bender. The bender has a cylindrically curved entry plane (Will be extended to also
* allow a flat entry plane).
*
* The bender may also be tapered such that it can have focusing or defocusing porperties.
* This may be specified through the "channel_file"-parameter, which points to a data file
* in which the individual channels are defined by an entry- and an exit-width.
*
* The file format for channel file columns:
* > #d-spacing1 d-spacing2 l_spacer d-coating1 d-coating2
* > 1e-3  0.8e-3 1e-5  1e-9 1e-9
*
* Each blades is considered coated with a reflecting coating, whose thickness is defined in the
* latter two columsn of the channel data file.
*
* Coating reflectivities are read from
* another data file with separate reflectivities for up and down spin. When a neutron ray hits a blade
* the polarization vector associated with the ray is decomposed into spin-up and down
* probablities, which in turn determines the overall reflectivity for this ray on the mirror.
* Furthermore the probabilities are used to also reconstruct the polarization vector of the reflected
* ray.
*
* File format for the reflectivity file:
* > #parameter = m
* > q/m R(up) R(down)
*
* The effect of spacers inside the channels is modelled by absorption only. The depth of the spacers
* in a channel is considered constant for one channel (set in the channel data file), the number of
* them is constant across all channels and is an input parameter to the component.
*
* %P
*
* INPUT PARAMETERS:
*
* xwidth: [m]        Width at the bender entry. Currently unsupported.
* yheight: [m]       Height at the bender entry.
* length: [m]        Length of blades that make up the bender.
* d_substrate: [m]   Thickness of the substrate.
* entry_radius: [m]  Radius of curvature of the entrance window.
* radius: [m]        Curvature of a single plate/polarizer. +:curve left/-:right
* nslit:  [1]        Number of slits in the bender
* channel_file: [ ]  File name of file containing channel data. such as spacer widths etc. (see above for format)
* reflect_file: [ ]  File name of file containing reflectivity data parametrized either by q or m.
* sigma_abs: [barn]  Absorption per unit cell at v=2200 m/s of spacers. Defaults to literature value for Al.
* V0: [AA^3]         Volume of unit cell for spacers. Defaults to Al.
* nspacer: [ ]       Number of spacers per channel.
* G:     [m/s^2]     Magnitude of gravitation.
* debug: [ ]         If > 0 print out some internal parameters.
*
* %L
*
* %E
****************************************************************************/
DEFINE COMPONENT Pol_bender_tapering

SETTING PARAMETERS (string channel_file="channel.dat", xwidth=0, yheight, length,
        d_substrate, entry_radius=10, radius=100, G=9.8, int nslit=1, int debug=1,
        string reflect_file=NULL, sigma_abs=0.231, V0=66.4, nspacer=5)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
  %include "pol-lib"
  %include "ref-lib"
  %include "read_table-lib"
  #ifndef pol_reflect_spin_h
  #define pol_reflect_spin_h 1

  /***************************************************************************
   * Calculate new spin vector after reflection in a surface with reflectivities
   * Rup and Rdown referring to the direction of the vector (mx,my,mz)
   * m is assumed to have unit length, - if zero the default m=(0,1,0) is used.
   * returns the fraction of intensity thus reflected (if Rup+Rdown<1)
   **************************************************************************/
  double
  pol_reflect_spin_ref (double Rup, double Rdown, double* sx, double* sy, double* sz, double mx, double my, double mz) {
    double pm, ps1, ps2, pin;
    double s1x, s1y, s1z, s2x, s2y, s2z;
    if (mx == my && my == mz && mz == 0.0) {
      pm = *sy;
    } else {
      pm = scalar_prod (*sx, *sy, *sz, mx, my, mz);
      /*find two vectors perpendicular to m*/
      vec_prod (s1x, s1y, s1z, mx, my, mz, 1.0, 1.0, 1.0);
      NORM (s1x, s1y, s1z);
      vec_prod (s2x, s2y, s2z, mx, my, mz, s1x, s1y, s1z);
      vec_prod (s1x, s1y, s1z, mx, my, mz, s2x, s2y, s2z);

      ps1 = scalar_prod (*sx, *sy, *sz, s1x, s1y, s1z);
      ps2 = scalar_prod (*sx, *sy, *sz, s2x, s2y, s2z);
    }
    /*the fractions of "up" and "down" neutrons*/
    double phiu = ((pm) + 1) / 2.0;
    double phid = (1 - (pm)) / 2.0;

    /*the relative fractions generate the new fraction*/
    if (phiu * Rup + phid * Rdown) {
      pm = (phiu * Rup - phid * Rdown) / (phiu * Rup + phid * Rdown);
      pin = phiu * Rup + phid * Rdown;
      // printf("pin: %g\n",pin);
    } else {
      pm = 0;
      pin = 0;
    }
    if (mx == my && my == mz && mz == 0.0) {
      *sy = pm;
    } else {
      /*this may be wrong*/
      *sx = pm * mx + ps1 * s1x + ps2 * s2x;
      *sy = pm * my + ps1 * s1y + ps2 * s2y;
      *sz = pm * mz + ps1 * s1z + ps2 * s2z;
      NORM (*sx, *sy, *sz);
    }
    return pin;
  }
  /***************************************************************************
   * Calculate new spin vector after reflection in a surface with reflectivities
   * Rup and Rdown referring to the y-axis (up)
   *****************************************************************/
  double
  pol_reflect_spin (double Rup, double Rdown, double* sx, double* sy, double* sz) {
    return pol_reflect_spin_ref (Rup, Rdown, sx, sy, sz, 0.0, 0.0, 0.0);
  }

  #endif

  #ifndef to_channel_frame
   #define to_channel_frame(rot,offset) \
    do {\
      mccoordschange(offset,rot,&x,&y,&z,&vx,&vy,&vz,&sx,&sy,&sz);\
    }while(0);
  #endif

  #ifndef from_channel_frame
   #define from_channel_frame(rot,offset)\
    do {\
      double cx,cy,cz;\
      Rotation trot;\
      rot_transpose(rot,trot);\
      coords_get(offset,&cx,&cy,&cz);\
      x=x-cx;\
      y=y-cy;\
      z=z-cz;\
      mccoordschange_polarisation(trot,&x,&y,&z);\
      mccoordschange_polarisation(trot,&vx,&vy,&vz);\
      mccoordschange_polarisation(trot,&sx,&sy,&sz);\
    }while(0);
  #endif
%}

DECLARE
%{
  t_Table T;
  t_Table R;
  t_Table C;
  int by_q;
  double rw;
  double rw_q;
  double phi;
  int bounce;
%}

INITIALIZE
%{
  int status, i;

  if (xwidth && entry_radius) {
    fprintf (stderr, "Warning (%s): Both xwidth and entry_radius set. Entry_radius overrides xwidth\n", NAME_CURRENT_COMP);
    xwidth = 0;
  }

  if ((status = Table_Read (&(T), channel_file, 0)) == -1) {
    fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", channel_file, NAME_CURRENT_COMP);
    exit (-1);
  }

  /*compute total opening and exit face of bender*/
  rw = (T.rows + 1) * d_substrate;   /*the substrate thickness * number of blades*/
  rw_q = (T.rows + 1) * d_substrate; /*the substrate thickness * number of blades*/
  for (i = 0; i < T.rows; i++) {
    rw += T.data[i * T.columns + 0]; /*the channel spacing*/
    rw += T.data[i * T.columns + 3]; /*the positive side coating*/
    rw += T.data[i * T.columns + 4]; /*the negative side coating*/

    rw_q += T.data[i * T.columns + 1]; /*channel exit spacing*/
    rw_q += T.data[i * T.columns + 3]; /*the positive side coating*/
    rw_q += T.data[i * T.columns + 4]; /*the negative side coating*/
  }

  /*angular opening*/
  phi = rw / entry_radius;

  if ((status = Table_Read (&(R), reflect_file, 0)) == -1) {
    fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", reflect_file ? reflect_file : "", NAME_CURRENT_COMP);
    exit (-1);
  }
  char** header_parsed;
  header_parsed = Table_ParseHeader (R.header, "parameter");
  if (strstr (header_parsed[0], "m")) {
    by_q = 0;
  } else {
    by_q = 1;
  }

  /*Now do computations on the channel edge curvatures and store in a table
   *We need 2 cols for entry position, 2 for q (flat face exit pos), 2 for centre of curvature and 2 for exit pos
   * After this the channel details can be gotten by Table_index(C,coli,channel_no) and Table_Index(C,coli,channel_no+)
   * where the latter is the edge on the positive side.*/
  /*initialize the table*/
  Table_Init (&(C), 8, T.rows);
%}

TRACE
%{
  int status, channel_no;
  double t0, t1;
  double rx, dx1, dx2;
  double xo, yo, zo, vox, voy, voz;
  double outer_radius;
  double dx1_q, dx2_q;
  /*propagate to entry surface (cylinder!)*/
  /*check which channel we're in, if any*/
  bounce = 0;
  channel_no = 0;

  if (xwidth) {
    PROP_Z0;
    if (y < -0.5 * yheight || y > 0.5 * yheight || x < -0.5 * xwidth || x > 0.5 * xwidth) {
      /*Missed flat face bender. Leave neutron be*/
      channel_no = -1;
      RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
    }
    fprintf (stderr, "Error(%s): flat face bender tapering is not supported yet. Please use a large entry_raidus instead. Aborting.\n", NAME_CURRENT_COMP);
    exit (-1);
  } else if (entry_radius) {
    status = cylinder_intersect (&t0, &t1, x, y, z - (-entry_radius), vx, vy, vz, entry_radius, 1000);
    if (entry_radius > 0) {
      PROP_DT (t1);
    } else {
      PROP_DT (t0);
    }
    xo = x;
    yo = y;
    zo = z;
    vox = vx;
    voy = vy;
    voz = vz;
    if (!status || t1 < 0 || y < -0.5 * yheight || y > 0.5 * yheight || x < entry_radius * sin (-0.5 * rw / entry_radius)
        || x > entry_radius * sin (0.5 * rw / entry_radius)) {
      /*Missed curved face bender. leave neutron be*/
      channel_no = -1;
      RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
    }
  }
  if (channel_no != -1) {
    /*So we actually hit the bender face - proceed*/
    double cphi, sphi, d1, d2, p0x, p0z, p1x, p1z, p2x, p2z, c1x, c1z, c2x, c2z, q0x, q0z, q1x, q1z, q2x, q2z, r1x, r1z, r2x, r2z;
    if (entry_radius) {
      channel_no = 0;
      /*offsets of channel entry - initialize
       * outer edge (dx2) : bender end + 1 substrate thickness + coating (outer)
       * inner edge (dx1) : outer edge + channel width*/
      dx2 = -0.5 * rw + d_substrate + T.data[channel_no * T.columns + 3];
      if (x < dx2) {
        /*hit first substrate*/
        ABSORB;
      }
      dx1 = dx2 + T.data[channel_no * T.columns + 0];

      cphi = cos (length / radius);
      sphi = sin (length / radius);

      q0x = radius * (1 - cos (length / radius)); /*exit point of bender centre*/
      q0z = radius * (sin (length / radius));

      /*the plates will end on a cirle of this radius*/
      outer_radius = entry_radius + length; // sqrt(q0x*q0x + (q0z+entry_radius)*(q0z+entry_radius));

      p2x = entry_radius * sin (dx2 / entry_radius); /*channel entry coordinates*/
      p2z = entry_radius * (cos (dx2 / entry_radius) - 1);
      p1x = entry_radius * sin (dx1 / entry_radius);
      p1z = entry_radius * (cos (dx1 / entry_radius) - 1);

      dx2_q = -0.5 * rw_q + d_substrate + T.data[channel_no * T.columns + 3];
      /*add the offset due to blade curvature to dx2_q*/
      // dx2_q += asin(q0x/outer_radius)*outer_radius;
      dx1_q = dx2_q + T.data[channel_no * T.columns + 1];

      q2x = q0x - dx2_q * (-cphi); /*exit point of outer channel edge (without correction for channel blade length)*/
      q2z = q0z - dx2_q * (sphi);
      q1x = q0x - dx1_q * (-cphi); /*exit point of inner channel edge (without correction for channel blade length)*/
      q1z = q0z - dx1_q * (sphi);

      #ifdef MCDEBUG
      printf ("%sPTS: %i %g %g %g %g %g %g %g %g %g %g\n", NAME_CURRENT_COMP, channel_no, p1x, p1z, p2x, p2z, q1x, q1z, q2x, q2z, q0x, q0z);
      #endif

      while (!(x > entry_radius * sin (dx2 / entry_radius) && x < entry_radius * sin (dx1 / entry_radius))) {
        channel_no++;
        if (channel_no > T.rows) {
          // printf("oops I've missed all channels ");
          // printf("%lld %g %g %g %g %g %g %g\n",mcget_run_num(),x,y,z,vx,vy,vz,p);
          ABSORB;
        }
        /*offsets of channel entry - update
         * outer edge (dx2) : previous inner edge + 1 coating (inner) + substrate thickness + coating (outer)
         * inner edge (dx1) : outer edge + channel width*/
        dx2 = dx1 + T.data[channel_no * T.columns + 3] + d_substrate + T.data[channel_no * T.columns + 4];
        /*now check if we've hit the blade*/
        // if ( x< dx2 && x>dx1 ){
        /*hit a substrate*/
        //  ABSORB;
        //}
        dx1 = dx2 + T.data[channel_no * T.columns + 0];

        dx2_q = dx1_q + d_substrate + T.data[channel_no * T.columns + 3] + T.data[channel_no * T.columns + 4];
        dx1_q = dx2_q + T.data[channel_no * T.columns + 1];

        p2x = entry_radius * sin (dx2 / entry_radius);
        p2z = entry_radius * (cos (dx2 / entry_radius) - 1);
        p1x = entry_radius * sin (dx1 / entry_radius);
        p1z = entry_radius * (cos (dx1 / entry_radius) - 1);

        q2x = q0x - dx2_q * (-cphi); /*exit point of outer channel edge (without correction for channel blade length)*/
        q2z = q0z - dx2_q * (sphi);
        q1x = q0x - dx1_q * (-cphi); /*exit point of inner channel edge (without correction for channel blade length)*/
        q1z = q0z - dx1_q * (sphi);

        q1x = outer_radius * sin (dx1_q / outer_radius);
        q1z = outer_radius * (cos (dx1_q / outer_radius)) - entry_radius;
        q2x = outer_radius * sin (dx2_q / outer_radius);
        q2z = outer_radius * (cos (dx2_q / outer_radius)) - entry_radius;
        double h1, D1, h2, D2;
        /*helper variables to to compute intersection points of circles*/
        /*We need to put h1 and h2  to the side the bender is curving to - hence the sign operation*/
        D1 = sqrt ((q1x - p1x) * (q1x - p1x) + (q1z - p1z) * (q1z - p1z));
        h1 = (radius < 0 ? -1 : 1) * sqrt (radius * radius - 0.25 * D1 * D1);
        D2 = sqrt ((q2x - p2x) * (q2x - p2x) + (q2z - p2z) * (q2z - p2z));
        h2 = (radius < 0 ? -1 : 1) * sqrt (radius * radius - 0.25 * D2 * D2);

        c1x = p1x + 0.5 * (q1x - p1x) + h1 / D1 * (q1z - p1z);  /*these are the centers around which the channel edges actually rotate*/
        c1z = p1z + 0.5 * (q1z - p1z) + h1 / D1 * -(q1x - p1x); /*positive*/

        c2x = p2x + 0.5 * (q2x - p2x) + h2 / D2 * (q2z - p2z); /*negative*/
        c2z = p2z + 0.5 * (q2z - p2z) + h2 / D2 * -(q2x - p2x);

        /*now find the exit points to place exit plane*/
        r1x = c1x + cphi * (p1x - c1x) + sphi * (p1z - c1z);
        r1z = c1z - sphi * (p1x - c1x) + cphi * (p1z - c1z);

        r2x = c2x + cphi * (p2x - c2x) + sphi * (p2z - c2z);
        r2z = c2z - sphi * (p2x - c2x) + cphi * (p2z - c2z);
        double nx, ny, nz;
        vec_prod (nx, ny, nz, 0.0, 1.0, 0.0, r2x - r1x, 0, r2z - r1z);
        // printf("DEBUG: %d   %g %g %g %g     %g %g %g %g     %g %g %g %g   %g %g %g %g   %g %g %g %g   %g %g %g %g\n",channel_no,p1x,p1z,p2x,p2z,
        // r1x,r1z,r2x,r2z, q1x,q1z,q2x,q2z,outer_radius, q0x,q0z, asin(q0x/outer_radius)*outer_radius, c1x,c1z,c2x,c2z,dx1,dx2,dx1_q,dx2_q);
      }

      /*set a scatter upon entry*/
      SCATTER;

    } else {
      /*so it is a flat face bender*/
      channel_no = 0;
      dx2 = -0.5 * rw;
      dx1 = dx2 + T.data[channel_no * T.columns + 0] + T.data[channel_no * T.columns + 4] + T.data[channel_no * T.columns + 3];
      while (!(x > dx2 && x < dx1)) {
        channel_no++;
        dx2 = dx1;
        dx1 += T.data[channel_no * T.columns + 0] + T.data[channel_no * T.columns + 4] + T.data[channel_no * T.columns + 3];
      }
      rx = (dx1 + dx2) * 0.5;
      SCATTER;
    }

    /*now the channel number is known. Given the central curvature of bender we may compute the placement of inner and outer cylinders. */
    /*Some helper points and variables:
      q is the end point of the uncurved blade.
      r is the end point of the curved blade
      we get the offset of the curved endpoint from the central blade, and apply it to vectors
      parallel and perpendicular to QP.*/
    q1x = outer_radius * sin (dx1_q / outer_radius);
    q1z = outer_radius * (cos (dx1_q / outer_radius)) - entry_radius;
    q2x = outer_radius * sin (dx2_q / outer_radius);
    q2z = outer_radius * (cos (dx2_q / outer_radius)) - entry_radius;

    double L1, L2;
    L1 = sqrt ((q1x - p1x) * (q1x - p1x) + (q1z - p1z) * (q1z - p1z));
    L2 = sqrt ((q2x - p2x) * (q2x - p2x) + (q2z - p2z) * (q2z - p2z));
    r1x = p1x + q0z * (q1x - p1x) / L1 + q0x * (q1z - p1z) / L1;
    r1z = p1z + q0z * (q1z - p1z) / L1 + q0x * -(q1x - p1x) / L1;
    r2x = p2x + q0z * (q2x - p2x) / L2 + q0x * (q2z - p2z) / L2;
    r2z = p2z + q0z * (q2z - p2z) / L2 + q0x * -(q2x - p2x) / L2;

    double h1, D1, h2, D2;
    /*helper variables to to compute intersection points of circles*/
    /*We need to put h1 and h2  to the side the bender is curving to - hence the sign operation*/
    D1 = sqrt ((r1x - p1x) * (r1x - p1x) + (r1z - p1z) * (r1z - p1z));
    h1 = (radius < 0 ? -1 : 1) * sqrt (radius * radius - 0.25 * D1 * D1);
    D2 = sqrt ((r2x - p2x) * (r2x - p2x) + (r2z - p2z) * (r2z - p2z));
    h2 = (radius < 0 ? -1 : 1) * sqrt (radius * radius - 0.25 * D2 * D2);

    c1x = p1x + 0.5 * (r1x - p1x) + h1 / D1 * (r1z - p1z);  /*these are the centers around which the channel edges actually rotate*/
    c1z = p1z + 0.5 * (r1z - p1z) + h1 / D1 * -(r1x - p1x); /*positive*/

    c2x = p2x + 0.5 * (r2x - p2x) + h2 / D2 * (r2z - p2z); /*negative*/
    c2z = p2z + 0.5 * (r2z - p2z) + h2 / D2 * -(r2x - p2x);

    /*now find the exit points to place exit plane*/
    /*    r1x=c1x + cphi*(p1x-c1x) + sphi*(p1z-c1z);*/
    /*    r1z=c1z - sphi*(p1x-c1x) + cphi*(p1z-c1z);*/
    /**/
    /*    r2x=c2x + cphi*(p2x-c2x) + sphi*(p2z-c2z);*/
    /*    r2z=c2z - sphi*(p2x-c2x) + cphi*(p2z-c2z);*/
    double nx, ny, nz;
    vec_prod (nx, ny, nz, 0.0, 1.0, 0.0, r2x - r1x, 0, r2z - r1z);
    #ifdef MCDEBUG
    printf ("%sPTS2: %d %g %g %g %g\n", NAME_CURRENT_COMP, channel_no, r1x, r1z, r2x, r2z);
    #endif
    // printf("DEBUG: %d   %g %g %g %g     %g %g %g %g     %g %g %g %g   %g %g %g %g   %g %g %g %g\n",channel_no,p1x,p1z,p2x,p2z, r1x,r1z,r2x,r2z,
    // q1x,q1z,q2x,q2z,outer_radius, q0x,q0z, asin(q0x/outer_radius)*outer_radius, c1x,c1z,c2x,c2z);

    int exit = 0;
    int coat1, coat2;                               /*are the blades coated at all*/
    coat1 = (Table_Index (T, channel_no, 3) > 0.0); /*positive side coating*/
    coat2 = (Table_Index (T, channel_no, 4) > 0.0); /*negative side coating*/
    do {
      int i1, i2, ic, ie, cyl;
      double t10, t11, t20, t21, tc, te;

      /*initialize times to seomthing known*/
      t10 = t11 = 0;
      t20 = t21 = 0;
      te = 0;
      i1 = cylinder_intersect (&t10, &t11, x - c1x, y, z - c1z, vx, vy, vz, radius, yheight);
      i2 = cylinder_intersect (&t20, &t21, x - c2x, y, z - c2z, vx, vy, vz, radius, yheight);
      ie = plane_intersect (&te, x, y, z, vx, vy, vz, nx, ny, nz, r1x, 0, r1z);
      /*catch the different cases*/
      if ((!ie) && (!i1 && !i2)) {
        fprintf (stderr, "Pol_bender (%s): Neutron xyz=(%g %g %g), v=(%g %g %g) cannot exit the bender or does not intersect either edge cylinder\n",
                 NAME_CURRENT_COMP, x, y, z, vx, vy, vz);
        ABSORB; /*this should really not happen*/
      }

      /*find the smallest strictly positive intersection time of te, t21 and t10*/
      double tt = FLT_MAX;
      /*first we mask strictly nonpositive times with a very large number (FLT_MAX)*/
      t10 = t10 <= 0 ? FLT_MAX : t10;
      t21 = t21 <= 0 ? FLT_MAX : t21;
      te = te <= 0 ? FLT_MAX : te;

      /*if radius<0 the edge on the positve x side becomes the outer one.
       * This means we chould compare t11 with t20, so in that case swap, and the algorithm should work*/
      if (radius < 0) {
        t10 = t11;
        t21 = t20;
      }

      enum { TOPBOTTOM, ENDFACE, CYL, CYL2, CYL1, UNKNOWN } branch;

      if (i2 && (t21 < t10)) {
        tt = t21;
        cyl = 2;
        branch = CYL;
      } else if (i1) {
        tt = t10;
        cyl = 1;
        branch = CYL;
      }
      if (te < tt) {
        tt = te;
        cyl = 0;
        branch = ENDFACE;
      }
      if (i2 & 24) {
        /* bitwise and with 16+8 to catch exit codes from cylinder_intersect. These mean that
         * neutron _exits_ through the top or bottom of the cylinder. Only cyl. 2 needs to be checked
         * since by design the neutron is always outside cyl. 1.*/
        branch = TOPBOTTOM;
        tt = 0;
        ABSORB;
      }
      if (tt == FLT_MAX) {
        tt = 0;
        branch = UNKNOWN;
        fprintf (stderr, "Cannot determine reflection branch: t21=%g, t10=%g ,te=%g\n", t21, t10, te);
        ABSORB;
      }

      switch (branch) {
      case ENDFACE:
        /*exit through channel end - leave neutron be*/
        exit = 1;
        break;
      case TOPBOTTOM:
        /*exit through top or bottom - leave neutron be*/
        exit = 1;
        break;
      case CYL:
        /*we bounce on bender wall*/
        PROP_DT (tt);
        double nx, ny, nz, s, alpha, Q, Rup, Rdown;
        if (cyl == 1) {
          /*check if coated*/
          if (!coat1) {
            /*not coated on this side*/
            ABSORB;
          }
          double vox = vx;
          double voy = vy;
          double voz = vz;
          double v = sqrt (scalar_prod (vx, vy, vz, vx, vy, vz));
          /*bounce on cylinder 1*/
          nx = x - c1x;
          ny = 0;
          nz = z - c1z;
          NORM (nx, ny, nz);
          s = scalar_prod (vx, vy, vz, nx, ny, nz);
          vx = vx - s * 2 * nx;
          /*vy=vy - s*2*ny;*/
          vz = vz - s * 2 * nz;
          /*compute reflectivity based on coating (m=3.2?)*/
          alpha = acos (scalar_prod (vx, vy, vz, vox, voy, voz) / v / v);
          /* Now compute reflectivity. */
          Q = 2.0 * sin (alpha / 2.0) * sqrt (scalar_prod (vx, vy, vz, vx, vy, vz)) * V2K;

          if (!by_q) {
            Q = Q / 0.0217; /*convert to m-value by conversion factor 2.17 AA^-1*/
          }
          Rup = Table_Value (R, Q, 1);
          Rdown = Table_Value (R, Q, 2);
          p *= pol_reflect_spin (Rup, Rdown, &sx, &sy, &sz);
        } else {
          if (!coat2) {
            /*not coated on this side*/
            ABSORB;
          }
          double vox = vx;
          double voy = vy;
          double voz = vz;
          double v = sqrt (scalar_prod (vx, vy, vz, vx, vy, vz));
          /*bounce on cylinder 2*/
          nx = x - c2x;
          ny = 0;
          nz = z - c2z;
          NORM (nx, ny, nz);
          s = scalar_prod (vx, vy, vz, nx, ny, nz);
          vx = vx - s * 2 * nx;
          /*vy=vy - s*2*ny;*/
          vz = vz - s * 2 * nz;
          /*compute reflectivity based on coating (m=3.2?)*/
          alpha = acos (scalar_prod (vx, vy, vz, vox, voy, voz) / v / v);
          /* Now compute reflectivity. */
          Q = 2.0 * sin (alpha / 2.0) * sqrt (scalar_prod (vx, vy, vz, vx, vy, vz)) * V2K;

          if (!by_q) {
            Q = Q / 0.0217; /*convert to m-value by conversion factor 2.17 AA^-1*/
          }
          Rup = Table_Value (R, Q, 1);
          Rdown = Table_Value (R, Q, 2);
          p *= pol_reflect_spin (Rup, Rdown, &sx, &sy, &sz);
        }
        /*if trace is set - transform back to bender frame, set SCATTER, and transform back to channel frame again*/
        bounce++;
        SCATTER;

        break;
      }
    } while (!exit);

    /*add spacer absorption*/
    double v, pmul;
    v = sqrt (scalar_prod (vx, vy, vz, vx, vy, vz));
    pmul = exp (-(nspacer * Table_Index (T, channel_no, 2)) * (sigma_abs * 2200 * 100 / V0 / v));
    p *= pmul;

    SCATTER;
  }
%}

MCDISPLAY
%{

  /*first draw the outer edges*/
  double yh2 = yheight / 2.0;
  int channel_no = 0;
  int N = 16;
  if (channel_no != -1) {
    /*So we actually hit the bender face - proceed*/
    double cphi, sphi, d1, d2, p0x, p0z, p1x, p1z, p2x, p2z, c1x, c1z, c2x, c2z, q0x, q0z, q1x, q1z, q2x, q2z, r1x, r1z, r2x, r2z;
    double dx1, dx2, dx1_q, dx2_q;
    channel_no = 0;
    dx2 = -0.5 * rw + d_substrate + T.data[channel_no * T.columns + 3];
    dx1 = dx2 + T.data[channel_no * T.columns + 0];

    dx2_q = -0.5 * rw_q + d_substrate + T.data[channel_no * T.columns + 3];
    dx1_q = dx2_q + T.data[channel_no * T.columns + 1];

    for (channel_no = 0; channel_no < T.rows; channel_no++) {
      if (T.rows > 20 && channel_no != 0 && channel_no != T.rows - 1 && channel_no % 100 != 0) {
        /*update the relative coordinates of the channel entry and exits but don't actually draw anything*/
        dx2 = dx1 + T.data[channel_no * T.columns + 3] + d_substrate + T.data[channel_no * T.columns + 4];
        dx1 = dx2 + T.data[channel_no * T.columns + 0];

        dx2_q = dx1_q + d_substrate + T.data[channel_no * T.columns + 3] + T.data[channel_no * T.columns + 4];
        dx1_q = dx2_q + T.data[channel_no * T.columns + 1];
        continue;
      }

      if (entry_radius) {
        /*offsets of channel entry - initialize
         * outer edge (dx2) : bender end + 1 substrate thickness + coating (outer)
         * inner edge (dx1) : outer edge + channel width*/
        p2x = entry_radius * sin (dx2 / entry_radius); /*channel entry coordinates*/
        p2z = entry_radius * (cos (dx2 / entry_radius) - 1);
        p1x = entry_radius * sin (dx1 / entry_radius);
        p1z = entry_radius * (cos (dx1 / entry_radius) - 1);
      }

      q0x = radius * (1 - cos (length / radius)); /*exit point of bender centre*/
      q0z = radius * (sin (length / radius));

      cphi = cos (length / radius);
      sphi = sin (length / radius);

      q2x = q0x - dx2_q * (-cphi); /*exit point of outer channel edge (without correction for channel blade length)*/
      q2z = q0z - dx2_q * (sphi);
      q1x = q0x - dx1_q * (-cphi); /*exit point of inner channel edge (without correction for channel blade length)*/
      q1z = q0z - dx1_q * (sphi);

      double h1, D1, h2, D2; /*helper variables to to compute intersection points of circles*/
      D1 = sqrt ((q1x - p1x) * (q1x - p1x) + (q1z - p1z) * (q1z - p1z));
      h1 = (radius < 0 ? -1 : 1) * sqrt (radius * radius - 0.25 * D1 * D1);
      D2 = sqrt ((q2x - p2x) * (q2x - p2x) + (q2z - p2z) * (q2z - p2z));
      h2 = (radius < 0 ? -1 : 1) * sqrt (radius * radius - 0.25 * D2 * D2);

      c1x = p1x + 0.5 * (q1x - p1x) + h1 / D1 * (q1z - p1z);  /*these are the centers around which the channel edges actually rotate*/
      c1z = p1z + 0.5 * (q1z - p1z) + h1 / D1 * -(q1x - p1x); /*positive*/

      c2x = p2x + 0.5 * (q2x - p2x) + h2 / D2 * (q2z - p2z); /*negative*/
      c2z = p2z + 0.5 * (q2z - p2z) + h2 / D2 * -(q2x - p2x);

      /*now we know what the circles are - draw N circle segments at +h/2 and -h/2 to mark channels*/
      double dl = length / N;
      double ro1x, ro1z, ro2x, ro2z;
      ro1x = p1x;
      ro1z = p1z;
      ro2x = p2x;
      ro2z = p2z;
      multiline (5, ro1x, yh2, ro1z, ro1x, -yh2, ro1z, ro2x, -yh2, ro2z, ro2x, yh2, ro2z, ro1x, yh2, ro1z);

      cphi = cos (dl / radius);
      sphi = sin (dl / radius);
      /*now find the exit points to place exit plane of circle segment*/
      r1x = c1x + cphi * (p1x - c1x) + sphi * (p1z - c1z);
      r1z = c1z - sphi * (p1x - c1x) + cphi * (p1z - c1z);

      r2x = c2x + cphi * (p2x - c2x) + sphi * (p2z - c2z);
      r2z = c2z - sphi * (p2x - c2x) + cphi * (p2z - c2z);
      /*loop until exit*/
      while (dl <= length) {
        line (ro1x, yh2, ro1z, r1x, yh2, r1z);
        line (ro2x, yh2, ro2z, r2x, yh2, r2z);
        line (ro1x, -yh2, ro1z, r1x, -yh2, r1z);
        line (ro2x, -yh2, ro2z, r2x, -yh2, r2z);

        /*update dl*/
        dl += length / N;
        /*store old exit in ro1x etc.*/
        ro1x = r1x;
        ro1z = r1z;
        ro2x = r2x;
        ro2z = r2z;

        cphi = cos (dl / radius);
        sphi = sin (dl / radius);
        /*now find the exit points to place exit plane of circle segment*/
        r1x = c1x + cphi * (p1x - c1x) + sphi * (p1z - c1z);
        r1z = c1z - sphi * (p1x - c1x) + cphi * (p1z - c1z);

        r2x = c2x + cphi * (p2x - c2x) + sphi * (p2z - c2z);
        r2z = c2z - sphi * (p2x - c2x) + cphi * (p2z - c2z);
      }
      /*draw the last segment*/
      line (ro1x, yh2, ro1z, r1x, yh2, r1z);
      line (ro2x, yh2, ro2z, r2x, yh2, r2z);
      line (ro1x, -yh2, ro1z, r1x, -yh2, r1z);
      line (ro2x, -yh2, ro2z, r2x, -yh2, r2z);
      multiline (5, ro1x, yh2, ro1z, ro1x, -yh2, ro1z, ro2x, -yh2, ro2z, ro2x, yh2, ro2z, ro1x, yh2, ro1z);

      /*update the relative coordinates of the channel entry and exits*/
      dx2 = dx1 + T.data[channel_no * T.columns + 3] + d_substrate + T.data[channel_no * T.columns + 4];
      dx1 = dx2 + T.data[channel_no * T.columns + 0];

      dx2_q = dx1_q + d_substrate + T.data[channel_no * T.columns + 3] + T.data[channel_no * T.columns + 4];
      dx1_q = dx2_q + T.data[channel_no * T.columns + 1];
    }
  }
%}

END