/*******************************************************************************
*
* 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: Apr 2003
* Origin: <a href="http://www.ill.fr">ILL (France)</a>.
* Modified by: G. Venturi, Made Honeycomb from Guide_gravity (Apr 2003)
* Modified by: E. Farhi, corrected gravitation bug (Feb 2005)
* Modified by: E. Farhi, uniformize parameter names (Jul 2008)
*
* 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; nslit,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, nslit=1, d=0.0005)
*
* %P
* INPUT PARAMETERS:
*
* w1: [m]          Width at the guide entry
* w2: [m]          Width at the guide exit. If zero, sets w2=w1.
* 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]           Thickness of subdividing walls [0]
* nslit: [1]       Number of horizontal channels in the guide (>= 1) [1]
*                 (nslit-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

SETTING PARAMETERS (w1, w2=0, l,
R0=0.995, Qc=0.0218, alpha=4.38, m=1.0, W=0.003, nslit=1, d=0.0005,
mleftup=-1, mrightup=-1, mleftdown=-1, mrightdown=-1,mleft=-1, mright=-1, G=0)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
  %include "ref-lib"
  #ifndef Honeycomb_guide_Version
  #define Honeycomb_guide_Version "$Revision$"

  #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];
    double warnings;
  } 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_nslit, 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;
    aVars->warnings = 0;

    if (a_nslit <= 0) {
      fprintf (stderr, "%s: Fatal: no channel in this guide (kh or nslit=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_nslit - 1) * 2 * a_d) / (double)a_nslit; /*INPUT APOTHEM*/
    aVars->w2c = 0.5 * (a_w2 - (a_nslit - 1) * 2 * a_d) / (double)a_nslit; /*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;
    }
  }

  #pragma acc routine seq
  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 nslit, 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, NULL, aVars->A[7], B, C);
    if (ret && dt0 > 10e-10) {
      dt_min = dt0;
      side = 7;
    }

    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, NULL, 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;
  // #prama acc routine declare create(GVars)
%}

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

  if (!w2)
    w2 = w1;

  if (W < 0 || nslit <= 0 || R0 < 0 || Qc < 0) {
    fprintf (stderr, "%s:Guide_gravity: W nslit 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, nslit, d, Gx, Gy, Gz, mright, mleft, mleftup, mrightdown, mrightup, mleftdown);

  if (!G)
    for (i = 0; i < 7; GVars.A[i++] = 0)
      ;
  // #pragma acc update device(GVars)
%}

TRACE
%{
  double B, C, dt;
  int ret, bounces = 0;
  float n, m1, nv, mv;
  int nup = -1, mright1 = -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 = solve_2nd_order (&dt, NULL, 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 (m1 = 0; y > ((2 * m1 + 1) * y_min); m1++)
      ;
    mv = m1;
    if (m1 == 0) {
      for (m1 = 0; y < ((2 * m1 - 1) * (y_min)); m1--)
        ;
      mv = m1;
    }

    yhole = 2 * m1 * 1.732 * (GVars.w1c + d);
    if (y - yhole > 0)
      mright1 = 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 += mright1 * (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, nslit, h_chnum);

      /* only positive dt are valid */
      /* exit reflection loops if no intersection (neutron is after box) */
      if (side == 0 || dt <= 0) {
        if (GVars.warnings < 100)
          fprintf (stderr, "%s: warning: neutron has entered guide, but can not exit !\n", GVars.compcurname);
        GVars.warnings++;
        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[side] * vx + GVars.ny[side] * vy + GVars.nz[side] * 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=|nslit_i-nslit_f| = V2Q*|vf - v| = V2Q*2*n.v/|n| */
      q = 2 * V2Q * fabs (B) / GVars.norm_n[side];
      {
        double par[] = { R0, Qc, alpha, GVars.M[side], W };
        StdReflecFunc (q, par, &B);
      }
      if (B <= 0) {
        x += w_adj;
        y += h_adj;
        ABSORB;
      } else
        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;
%}

FINALLY
%{
  if (GVars.warnings > 100) {
    fprintf (stderr, "%s: warning: neutron has entered guide, but can not exit !\n", GVars.compcurname);
    fprintf (stderr, "%s: warning: This message has been repeated %g times\n", GVars.compcurname, GVars.warnings);
  }
%}


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;

  for (j = -nslit / 2; j <= nslit / 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 = -nslit; i <= nslit; 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