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