/*******************************************************************************
*
* McStas, neutron ray-tracing pacxkage
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Guide_honeycomb
*
* %I
* Written by: G. Venturi
* Date: Aug 03 2001
* Version: $Revision: 1.6 $
* Origin: <a href="http://www.ill.fr">ILL (France)</a>.
* Release: McStas 1.6
* Modified by: E. Farhi, from Honeycomb_guide by K. Lefmann (buggy).
* Modified by: E. Farhi, focusing channels are now ok (Sept 4th, 2001).
* Modified by: E. Farhi, 2D channel array. Correct focus bug (Dec 14th 2002)
* Modified by: G. Venturi, Made Honeycomb from Guide_gravity (Apr 2003)
* Modified by: E. Farhi, corrected gravitation bug (Feb 2005)
*
* Neutron guide with gravity and honeycomb geometry. Can be channeled and focusing.
*
* %D
* Models a honeycomb guide tube centered on the Z axis. The entrance lies
* in the X-Y plane. Gravitation applies also when reaching the guide input
* window. The guide can be channeled (hexagonal channel section; k,d parameters).
* The guide coating specifications may be entered via different ways (global,
* or for each wall m-value).
* For details on the geometry calculation see the description in the McStas
* reference manual.
*
* Example: Guide_honeycomb(w1=0.1, w2=0.1, l=12,
*           R0=0.99, Qc=0.0219, alpha=6.07, m=1.0, W=0.003, k=1, d=0.0005)
*
* %P
* INPUT PARAMETERS:
*
* w1:      [m]    Width at the guide entry
* w2:      [m]    Width at the guide exit
* l:       [m]    length of guide
* R0:      [1]    Low-angle reflectivity
* Qc:      [AA-1] Critical scattering vector
* alpha:   [AA]   Slope of reflectivity
* m:       [1]    m-value of material. Zero means completely absorbing.
* W:       [AA-1] Width of supermirror cut-off
* d:       [m]    Thicxkness of subdividing walls [0]
* k:       [1]    Number of horizontal channels in the guide (>= 1) [1]
*                 (k-1 vertical dividing walls)
*
* Optional input parameters: (different ways for m-specifications)
*
* mright:    [1]    m-value of material for right.     vertical mirror
* mleft:     [1]    m-value of material for left.      vertical mirror
* mleftup:   [1]    m-value of material for leftup.    oblique mirror
* mrightup:  [1]    m-value of material for rightdown. oblique mirror
* mrightdown:[1]    m-value of material for leftup.    oblique mirror
* mleftdown: [1]    m-value of material for rightdown. oblique mirror
* G:         [m/s2] Gravitation norm. 0 value disables G effects.
*
* OUTPUT PARAMETERS
*
* GVars:              (1) internal variables
* GVars.N_reflection: (1) Array of the cumulated Number of reflections
*                   N_reflection[0] total nb of reflections
*                   N_reflection[1,2,3,4.5.6]  reflections
*                   N_reflection[7] total nb neutrons exiting guide
*                   N_reflection[8] total nb neutrons entering guide
*
* %D
* Example values: m=4 Qc=0.02 W=1/300 alpha=6.49 R0=1
*
* %E
*******************************************************************************/

DEFINE COMPONENT Guide_honeycomb
DEFINITION PARAMETERS ()
SETTING PARAMETERS (w1, w2, l,
  R0=0.99, Qc=0.0219, alpha=6.07, m=1.0, W=0.003, k=1, d=0.0005,
  mleftup=-1, mrightup=-1, mleftdown=-1, mrightdown=-1,mleft=-1, mright=-1, G=0)
OUTPUT PARAMETERS (GVars)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

SHARE
%{
#ifndef Honeycomb_guide_Version
#define Honeycomb_guide_Version "1.5"

#ifndef PROP_GRAV_DT
#error McStas : You need PROP_GRAV_DT (McStas >= 1.4.3) to run this component
#endif

/*
* G:       (m/s^2) Gravitation acceleration along y axis [-9.81]
* Gx:      (m/s^2) Gravitation acceleration along x axis [0]
* Gy:      (m/s^2) Gravitation acceleration along y axis [-9.81]
* Gz:      (m/s^2) Gravitation acceleration along z axis [0]
* mh:      (1)    m-value of material for left/right vert. mirrors
* mv:      (1)    m-value of material for top/bottom horz. mirrors
* mx:      (1)    m-value of material for left/right vert. mirrors
* my:      (1)    m-value of material for top/bottom horz. mirrors
*/

  typedef struct Honeycomb_guide_Vars
  {
    double gx;
    double gy;
    double gz;
    double nx[8], ny[8], nz[8];
    double wx[8], wy[8], wz[8];
    double A[8], norm_n2[8], norm_n[8];
    long   N_reflection[9];
    double w1c, w2c;
    double M[7];
    double nzC[7], norm_n2xy[7], Axy[7];
    char   compcurname[256];
  } Honeycomb_guide_Vars_type;

  void Honeycomb_guide_Init(Honeycomb_guide_Vars_type *aVars,
    MCNUM a_w1, MCNUM a_w2, MCNUM a_l, MCNUM a_R0,
    MCNUM a_Qc, MCNUM a_alpha, MCNUM a_m, MCNUM a_W, MCNUM a_k, MCNUM a_d,
    MCNUM a_Gx, MCNUM a_Gy, MCNUM a_Gz, MCNUM a_mright, MCNUM a_mleft, MCNUM a_mleftup,
    MCNUM a_mrightdown, MCNUM a_mrightup, MCNUM a_mleftdown)
  {
    int i;

    for (i=0; i<8; aVars->N_reflection[i++] = 0);
    for (i=0; i<7; aVars->M[i++] = 0);

    aVars->gx = a_Gx; /* The gravitation vector in the current component axis system */
    aVars->gy = a_Gy;
    aVars->gz = a_Gz;

    if (a_k <= 0) { fprintf(stderr,"%s: Fatal: no channel in this guide (kh or k=0).\n", aVars->compcurname); exit(-1); }
    if (a_d < 0)  { fprintf(stderr,"%s: Fatal: subdividing walls have negative thickness in this guide (d<0).\n", aVars->compcurname); exit(-1); }

    aVars->w1c = 0.5*(a_w1 - (a_k-1) *2*a_d)/(double)a_k;            /*INPUT APOTHEM*/
    aVars->w2c = 0.5*(a_w2 - (a_k-1) *2*a_d)/(double)a_k;            /*OUTPUT APOTHEM*/

    for (i=0; i <= 6;   aVars->M[i++]=a_m);
    if (a_mright     >= 0) aVars->M[1] =a_mright;
    if (a_mleft      >= 0) aVars->M[2] =a_mleft;
    if (a_mleftup    >= 0) aVars->M[3] =a_mleftup;
    if (a_mrightdown >= 0) aVars->M[4] =a_mrightdown;
    if (a_mrightup   >= 0) aVars->M[5] =a_mrightup;
    if (a_mleftdown  >= 0) aVars->M[6] =a_mleftdown;
    /* n: normal vectors to surfaces */
    aVars->nx[1] = -a_l;    aVars->ny[1] =  0;   aVars->nz[1] =  (aVars->w2c-aVars->w1c); /* 1:+X right     */
    aVars->nx[2] = +a_l;    aVars->ny[2] =  0;   aVars->nz[2] =  -aVars->nz[1];           /* 2:-X left      */

    aVars->nx[3] = +a_l*0.5;   aVars->ny[3] = -0.866025*a_l; aVars->nz[3] =  (aVars->w2c-aVars->w1c); /* 3:+Y leftup*/
    aVars->nx[4] = -a_l*0.5;   aVars->ny[4] = +0.866025*a_l; aVars->nz[4] = -(aVars->w2c-aVars->w1c); /* 4:+Y rightdown*/
    aVars->nx[5] = -a_l*0.5;   aVars->ny[5] = -0.866025*a_l; aVars->nz[5] =  (aVars->w2c-aVars->w1c); /* 5:+Y rightup   */
    aVars->nx[6] = +a_l*0.5;   aVars->ny[6] = +0.866025*a_l; aVars->nz[6] = -(aVars->w2c-aVars->w1c); /* 6:+Y leftdown */

    aVars->nx[7] =  0;   aVars->ny[7] =  0;   aVars->nz[7] =  a_l;
    aVars->nx[0] =  0;   aVars->ny[0] =  0;   aVars->nz[0] =  -a_l;
    /* w: a point on these surfaces */
    aVars->wx[1] =  (aVars->w1c);      aVars->wy[1] = 0;                           aVars->wz[1] = 0;   /* 1: right     */
    aVars->wx[2] =  -(aVars->w1c);     aVars->wy[2] = 0;                           aVars->wz[2] = 0;   /* 2: left      */
    aVars->wx[3] = -0.5*(aVars->w1c);  aVars->wy[3] = +0.866025*(aVars->w1c);      aVars->wz[3] = 0;   /* 3: leftup    */
    aVars->wx[4] = +0.5*(aVars->w1c);  aVars->wy[4] = -0.866025*(aVars->w1c);      aVars->wz[4] = 0;   /* 4: rightdown */
    aVars->wx[5] = +0.5*(aVars->w1c);  aVars->wy[5] = +0.866025*(aVars->w1c);      aVars->wz[5] = 0;   /* 5: rightup   */
    aVars->wx[6] = -0.5*(aVars->w1c);  aVars->wy[6] = -0.866025*(aVars->w1c);      aVars->wz[6] = 0;   /* 6: leftdown  */
    aVars->wx[7] =  0;                 aVars->wy[7] =  0;                          aVars->wz[7] = a_l; /* 7:+Z exit    */
    aVars->wx[0] =  0;                 aVars->wy[0] =  0;                          aVars->wz[0] = 0;   /* 0:Z0 input   */

    for (i=0; i <= 7; i++)
    {
      aVars->A[i] = scalar_prod(aVars->nx[i], aVars->ny[i], aVars->nz[i], aVars->gx, aVars->gy, aVars->gz)/2;
      aVars->norm_n2[i] = aVars->nx[i]*aVars->nx[i] + aVars->ny[i]*aVars->ny[i] + aVars->nz[i]*aVars->nz[i];
      if (aVars->norm_n2[i] <= 0)
     { fprintf(stderr,"%s: Fatal: normal vector norm %i is null/negative ! check guide dimensions.\n", aVars->compcurname, i); exit(-1); } /* should never occur */
      else
        aVars->norm_n[i] = sqrt(aVars->norm_n2[i]);
    }
    /* partial computations for sides, to save computing time */
    for (i=1; i <= 6; i++)
    {
      aVars->nzC[i]      = aVars->nz[i];
      aVars->norm_n2xy[i]= aVars->nx[i]*aVars->nx[i] + aVars->ny[i]*aVars->ny[i];
      aVars->Axy[i]      =(aVars->nx[i]*aVars->gx    + aVars->ny[i]*aVars->gy)/2;
    }
  }

  int Honeycomb_guide_Trace(double *dt,
        Honeycomb_guide_Vars_type *aVars,
        double cx, double cy, double cz,
        double cvx, double cvy, double cvz,
        double cxnum, int k, double cynum)
  {
    double B, C;
    int    ret=0;
    int    side=0;
    double n1;
    double dt0, dt_min=0;
    int i;
    double loc_num;
    int    i_slope=3;

    /* look if there is a previous intersection with guide sides */
    /* A = 0.5 n.g; B = n.v; C = n.(r-W); */
    /* 5=+Z side: n=(0, 0, -l) ; W = (0, 0, l) (at z=l, guide exit)*/
    B = aVars->nz[7]*cvz; C = aVars->nz[7]*(cz - aVars->wz[7]);
    ret = solve_2nd_order(&dt0, aVars->A[7], B, C);
    if (ret && dt0>10e-10)
    { dt_min = dt0; side=7; aVars->n_dot_v[7] = B; }

    loc_num = (3*cynum+cxnum)/2;
    for (i=6; i>0; i--)
    {
      if (i == 4) { i_slope=1; loc_num =(3*cynum-cxnum)/2; }
      if (i == 2) { i_slope=1; loc_num = cxnum;}

      if (aVars->nzC[i_slope] != 0) {
        n1=loc_num-1;
        loc_num++; loc_num++;
        aVars->nz[i] = aVars->nzC[i]*n1;
        aVars->A[i]  = aVars->Axy[i] + aVars->nz[i]*aVars->gz/2;
      }

      B = aVars->nx[i]*cvx + aVars->ny[i]*cvy + aVars->nz[i]*cvz;                            /* n.v */
      C = aVars->nx[i]*(cx-aVars->wx[i]) + aVars->ny[i]*(cy-aVars->wy[i]) + aVars->nz[i]*cz; /* n.(r-W) */

      ret = solve_2nd_order(&dt0, aVars->A[i], B, C);
      if (ret && dt0>10e-10 && (dt0<dt_min || !dt_min))
      { dt_min = dt0; side=i;
        if (aVars->nzC[i] != 0)
        { aVars->norm_n2[i] = aVars->norm_n2xy[i] + aVars->nz[i]*aVars->nz[i]; aVars->norm_n[i] = sqrt(aVars->norm_n2[i]); }
      }
     }

    *dt = dt_min;
    return (side);
  }
#endif
%}

DECLARE
%{
  Honeycomb_guide_Vars_type GVars;
%}

INITIALIZE
%{
  double Gx=0, Gy=9.81, Gz=0;
  Coords mcLocG;
  int i;

  if (W < 0 || k <= 0 || R0 < 0 || Qc < 0)
  { fprintf(stderr,"%s:Guide_gravity: W k R0 Qc must be >0.\n", NAME_CURRENT_COMP);
    exit(-1); }

  if (mcgravitation) G=-9.81;
  mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(Gx,G,Gz));
  coords_get(mcLocG, &Gx, &Gy, &Gz);

  strcpy(GVars.compcurname, NAME_CURRENT_COMP);
  Honeycomb_guide_Init(&GVars,
    w1, w2, l, R0,
    Qc, alpha, m, W, k, d,
    Gx, Gy, Gz,mright, mleft, mleftup, mrightdown, mrightup, mleftdown);

  if (!G) for (i=0; i<7; GVars.A[i++] = 0);

%}

TRACE
%{
  double B, C, dt;
  int    ret, bounces = 0;
  float    n,m,nv,mv;
  int    nup=-1, mright=-1;
  double xhole,yhole, y_min;
  double cn;
  double inside;
  double w_edge, w_adj; /* Channel displacement  */
  double h_edge, h_adj;
  double w_chnum,h_chnum,w_c,h_c; /* channel indexes */

  /* propagate to box input (with gravitation) in comp local coords */
  /* 0=Z0 side: n=(0, 0, l) ; W = (0, 0, 0) (at z=0, guide input)*/
  B = -l*vz; C = -l*z;

  ret = plane_intersect_Gfast(&dt, GVars.A[0], B, C);
  if (ret && dt>0)
  {
    PROP_GRAV_DT(dt, GVars.gx, GVars.gy, GVars.gz);
    GVars.N_reflection[8]++;
  }
  /* check if we are in the box input, else absorb */
  if(dt > 0 && fabs(x) <= w1/2 && fabs(y) <= w1/2)
  {
    for(n=0; x>((2*n+1)*(GVars.w1c+d)) ; n++); nv=n;
    if (n==0) { for(n=0; x<((2*n-1)*(GVars.w1c+d)); n--); nv=n; }

    xhole=2*n*(GVars.w1c+d);
    if(x-xhole>0)  nup=1;
    y_min=1.732*(GVars.w1c+d);

    for(m=0;y>((2*m+1)*y_min); m++); mv=m;
    if (m==0) { for (m=0; y<((2*m-1)*(y_min)); m--); mv=m; }

    yhole=2*m*1.732*(GVars.w1c+d);
    if(y-yhole>0)  mright=1;

    cn=1.1547*(GVars.w1c+d);
    inside=(fabs(y-yhole)+0.577351*fabs(x-xhole)-cn);
    if(inside>0)
    {
      xhole +=nup*(GVars.w1c+d);
      yhole +=mright*(1.732*(GVars.w1c+d));
    }

    /*   double w_chnum,h_chnum;*/ /* channel indexes  */
    /* Shift origin to center of channel hit (absorb if hit dividing walls) */
    w_c=xhole/(GVars.w1c+d);
    h_c=yhole/(1.732*(GVars.w1c+d));
    w_chnum=rint(w_c);
    h_chnum=rint(h_c);

    x -=xhole;
    y -=yhole;
    w_adj = xhole;
    h_adj = yhole;
    if(fabs(x) > GVars.w1c)
    {
      x += xhole; /* Re-adjust origin */
      y += yhole;
      ABSORB;
    }
    if(fabs(x*0.5+y*0.866025) > GVars.w1c)
    {
      x += xhole; /* Re-adjust origin */
      y += yhole;
      ABSORB;
    }
    if(fabs(-x*0.5+y*0.866025) > GVars.w1c)
    {
      x += xhole; /* Re-adjust origin */
      y += yhole;
      ABSORB;
    }

    /* neutron is now in the input window of the guide */
    /* do loops on reflections in the box */
    for(;;)
    {
      /* get intersections for all box sides */
      double q;
      int side=0;

      bounces++;
      /* now look for intersection with guide sides and exit */

      side = Honeycomb_guide_Trace(&dt, &GVars, x, y, z,
          vx, vy, vz, w_chnum, k, h_chnum);

      /* only positive dt are valid */
      /* exit reflection loops if no intersection (neutron is after box) */
      if (side == 0 || dt <= 0)
        { fprintf(stderr,"%s: warning: neutron has entered guide, but can not exit !\n", GVars.compcurname);
          x += w_adj; y += h_adj; ABSORB; } /* should never occur */

      /* propagate to dt */
      PROP_GRAV_DT(dt, GVars.gx, GVars.gy, GVars.gz);

      /* do reflection on speed for l/r/u/d sides */
      if (side == 7) /* neutron reaches end of guide: end loop and exit comp */
        { GVars.N_reflection[side]++; x += w_adj; y += h_adj; SCATTER; x -= w_adj; y -= h_adj; break; }
      /* else reflection on a guide wall */
      if(GVars.M[side] == 0 || Qc == 0)  /* walls are absorbing */
        { x += w_adj; y += h_adj; ABSORB; }
      /* change/mirror velocity: h_f = v - n.2*n.v/|n|^2 */
      B = GVars.nx[i]*vx + GVars.ny[i]*vy + GVars.nz[i]*vz;  /* n.v */
      GVars.N_reflection[side]++; /* GVars.norm_n2 > 0 was checked at INIT */
      dt = 2*B/GVars.norm_n2[side]; /* 2*n.v/|n|^2 */
      vx -= GVars.nx[side]*dt;
      vy -= GVars.ny[side]*dt;
      vz -= GVars.nz[side]*dt;

      /* compute q and modify neutron weight */
      /* scattering q=|k_i-k_f| = V2Q*|vf - v| = V2Q*2*n.v/|n| */
      q = 2*V2Q*fabs(B)/GVars.norm_n[side];
      B = R0;
      if(q > Qc)
      {
        double arg;
        if (W>0) arg = (q-GVars.M[side]*Qc)/W;
        else arg = (q-GVars.M[side]*Qc)*10000; /* W = 0.00001 */

        if(arg < 10) B *= .5*(1-tanh(arg))*(1-alpha*(q-Qc));
        else { x += w_adj; y += h_adj; ABSORB; }; /* Cutoff ~ 1E-10 */
       }
       if (B < 0) B=0;
       else if (B > 1) B=1;
       p *= B;
       x += w_adj; y += h_adj; SCATTER; x -= w_adj; y -= h_adj;
       GVars.N_reflection[0]++;
       /* go to the next reflection */
       if (bounces > 1000) ABSORB;
     } /* end for */
     x += w_adj; y += h_adj; /* Re-adjust origin after SCATTER */
  }
  else
    ABSORB;
%}


MCDISPLAY
%{
  int i,j;
  double a,b,c;
  double x0,x01,x1,x11,x2,x21;
  double y0,y01,y1,y11,y2,y21,y3,y31,y4,y41;

  magnify("xy");

  for(j = -k/2; j <= k/2; j++)
  {
    y0 =  j*(GVars.w1c+d)*1.732;
    y01=  j*(GVars.w2c+d)*1.732;
    y1 =y0 +(GVars.w1c+d)/1.732;
    y11=y01+(GVars.w2c+d)/1.732;
    y2 =y0 -(GVars.w1c+d)/1.732;
    y21=y01-(GVars.w2c+d)/1.732;
    y3 =y0 +(GVars.w1c+d)*1.1547;
    y31=y01+(GVars.w2c+d)*1.1547;
    y4 =y0 -(GVars.w1c+d)*1.1547;
    y41=y01-(GVars.w2c+d)*1.1547;

     for(i = -k; i <= k; i++)

    {

      a=i+j;
      b=a/2+0.1;
      c=rint(b);

      if(fabs(c-b)<0.3)

      {

       x0 =  i*(GVars.w1c+d);
       x01=  i*(GVars.w2c+d);
       x1 =x0 +(GVars.w1c+d);
       x11=x01+(GVars.w2c+d);
       x2 =x0 -(GVars.w1c+d);
       x21=x01-(GVars.w2c+d);


       multiline(5,
                x1, y1, 0.0,
                x11, y11, (double)l,
                x21, y11, (double)l,
                x2, y1, 0.0,
                x1, y1, 0.0);


       multiline(5,
                x1, y1, 0.0,
                x11, y11, (double)l,
                x01, y31, (double)l,
                x0, y3, 0.0,
                x1, y1, 0.0);

       multiline(5,
                x0, y3, 0.0,
                x01, y31, (double)l,
                x21, y11, (double)l,
                x2, y1, 0.0,
                x0, y3, 0.0);

       multiline(5,
                x2, y1, 0.0,
                x21, y11, (double)l,
                x21, y21, (double)l,
                x2, y2, 0.0,
                x2, y1, 0.0);

       multiline(5,
                x2, y2, 0.0,
                x21, y21, (double)l,
                x01, y41, (double)l,
                x0, y4, 0.0,
                x2, y2, 0.0);

       multiline(5,
                x0, y4, 0.0,
                x01, y41, (double)l,
                x11, y21, (double)l,
                x1, y2, 0.0,
                x0, y4, 0.0);

       }
     }
   }

%}

END