/******************************************************************
*
*  McStas, version 3.1, 
*          
*          
*
* Component: Fermi_chop2a
*
* %Identification
* Written by: Garrett Granroth
* Date: 6 Feb 2005
* Origin: SNS Oak Ridge,TN
* 
* SNS Fermi Chopper as used in SNS_ARCS
*
* %D
* Models an SNS Fermi Chopper, used in the SNS_ARCS instrument model.
*
* %Parameters
* nput Parameters:
*  
* len: [m]      slit package length
* w: [m]        slit package width
* nu: [Hz]      frequency
* delta: [sec]  time from edge of chopper to center Phase angle
* tc: [sec]     time when desired neutron is at the center of the chopper
* ymin: [m]     Lower y bound
* ymax: [m]     Upper y bound
* nchan: [1]    number of channels in chopper
* bw: [m]       blade width
* blader: [m]   blade radius
*
*
* %End
*******************************************************************/
DEFINE COMPONENT Fermi_chop2a

SETTING PARAMETERS (len, w, nu, delta, tc, ymin, ymax, nchan, bw, blader)

SHARE
%{
  #ifndef FERMI_CHOP_DEFS
  #define FERMI_CHOP_DEFS
  /* routine to calculate acos in proper quadrant  range = 0 to 2PI*/
  #pragma acc routine
  double
  acos0_2pi (double x, double y) {
    if (y > 0.0) {
      return acos (x);
    }
    return 2.0 * PI - acos (x);
  }

  /*routine to calculate x and y positions of a neutron in a fermi chopper */
  #pragma acc routine
  void
  neutxypos (double* x, double* y, double phi, double inrad, double* c) {
    *x = c[0] + inrad * cos (phi);
    *y = c[1] + inrad * sin (phi);
  }

  /* routine to calculate the origin of a circle that describes the neutron path through the chopper */
  #pragma acc routine
  void
  calccenter (double* c, double* zr, double* xr) {
    double denom, A, B, C, D, a, b;
    denom = 2 * (-zr[0] * xr[2] + zr[0] * xr[1] + zr[1] * xr[2] + xr[0] * zr[2] - xr[0] * zr[1] - xr[1] * zr[2]);
    A = xr[1] - xr[2];
    B = xr[0] - xr[1];
    C = zr[2] - zr[1];
    D = zr[1] - zr[0];
    a = zr[0] * zr[0] - zr[1] * zr[1] + xr[0] * xr[0] - xr[1] * xr[1];
    b = zr[2] * zr[2] - zr[1] * zr[1] + xr[2] * xr[2] - xr[1] * xr[1];
    c[0] = 1.0 / denom * (A * a + B * b);
    c[1] = 1.0 / denom * (C * a + D * b);
  }

  #endif

  /* function that describes the shape of the blades */
  #pragma acc routine
  double
  blades (double zin, double rin, double off) {
    if (rin != 0.0)
      return rin * (1 - cos (asin (zin / fabs (rin)))) + off;
    else
      return 0;
  }

  /* function to calculate which channel the neturon is in and to check if it is in a blade
   *  or outside the slit package
   * return 0 if neutron does not transmit return 1 with channel number if neutron will pass*/
  #pragma acc routine
  int
  checkabsorb (double phi, int* chan_num, double inrad, double inw, double insw, double inbw, double blader, double off, double* c) {
    double xtmp, neuzr, neuxr;
    neutxypos (&neuzr, &neuxr, phi, inrad, c);
    // printf("xr:%g zr:%g phi: %g r: %g c[0]: %g c[1]: %g\n",neuxr,neuzr,phi,inrad,c[0],c[1]);
    if (fabs (neuxr) > inw / 2.0) // check if neutron x position is outside of slit package
      return 0;
    xtmp = neuxr + inw / 2.0;                                                // move origin to side of slit package
    *chan_num = ceil ((xtmp - blades (neuzr, blader, off)) / (inbw + insw)); // calculate channel number
                                                                             // check if neutron is in blade
    if (xtmp > *chan_num * (inbw + insw) + blades (neuzr, blader, off))
      return 0;
    return 1;
  }
%}

    DECLARE
%{

  double omega;
  double off;
  double splen;
  double rad;
  double sw;
  double tw;
%}
INITIALIZE
%{
  splen = len / 2.0;
  omega = 2.0 * PI * nu;
  off = blader * (1 - cos (asin (splen / fabs (blader)))); // the additional width needed to accomodate the curvature of the blade
  tw = (w + 2.0 * off);                                    // the total width needed to contain the slit package
  rad = sqrt (tw * tw / 4.0 + splen * splen);              // radius of cylinder containing slit package.
  sw = (w - bw) / nchan - bw;
  printf ("sw: %g rad: %g\n", sw, rad);
%}
TRACE
%{

  double t0, t1, dphi, dt2, tneuzr, tneuxr, nrad;
  double phivec[200], tpt[3], xpt[3], ypt[3], zpt[3], zr[3], xr[3], yr[3], theta[3], c[2];
  int chan_num, chan_num0, idx1, idx3;
  if (cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, rad, ymax - ymin)) {
    if (t0 < 0) /*Neutron started inside cylinder */
      ABSORB;
    dt2 = t1 - t0;
    PROP_DT (t0); /*propagate neutron to edge of chopper*/
    /*calculate neutron position and velocity in chopper frame
      calculate 3 points in the instrument frame and put them into the
      chopper frame inorder to determine the radius and center of a circle
       that describes the path of the neutron in the chopper frame. */
    tpt[1] = t;
    tpt[2] = t + dt2;
    tpt[0] = t + dt2 / 2.0;
    // set local 0 in time as tc and calculate angle of rotation for each point
    for (idx3 = 0; idx3 < 3; idx3++) {
      theta[idx3] = (tpt[idx3] - tc) * omega;
    }
    zpt[1] = -sqrt (rad * rad - x * x);
    xpt[1] = x;
    ypt[1] = y; /* point where neutron intersects chopper */
    zpt[2] = zpt[1] + vz * (dt2);
    xpt[2] = xpt[1] + vx * (dt2);
    ypt[2] = ypt[1] + vy * (dt2); /* point where neutron leaves the chopper */
    xpt[0] = xpt[1] + vx * (dt2 / 2.0);
    ypt[0] = ypt[1] + vy * (dt2 / 2.0);
    zpt[0] = zpt[1] + vz * (dt2 / 2.0); /*point half way between in time */
    /* do the rotation */
    for (idx3 = 0; idx3 < 3; idx3++) {
      rotate (xr[idx3], yr[idx3], zr[idx3], xpt[idx3], ypt[idx3], zpt[idx3], theta[idx3], 0, 1, 0);
    }
    calccenter (c, zr, xr);                                                          /* calculate the center */
    nrad = sqrt ((zr[0] - c[0]) * (zr[0] - c[0]) + (xr[0] - c[1]) * (xr[0] - c[1])); /*calculate the radius of curvature for the neutron path */
    /* calculate points along path of neutron through cylinder quit on absorption
     * or transmit neutron if 200 points are calculated
     * calculate phi for first and last points */
    phivec[0] = acos0_2pi ((zr[1] - c[0]) / nrad, xr[1] - c[1]);
    phivec[1] = acos0_2pi ((zr[2] - c[0]) / nrad, xr[2] - c[1]);
    neutxypos (&tneuzr, &tneuxr, phivec[0], nrad, c);
    /* reset phi[0] and phi[1] to match the length of the slit package rather than cylinder radius*/
    if (tneuzr < -splen) {
      phivec[0] = acos0_2pi ((-c[0] - splen) / nrad, -c[1]);
    }
    neutxypos (&tneuzr, &tneuxr, phivec[1], nrad, c);
    if (tneuzr > splen) {
      phivec[1] = acos0_2pi ((-c[0] + splen / 2.0) / nrad, -c[1]);
    }
    dphi = phivec[1] - phivec[0]; /* initial dphi */
    idx1 = 2;
    phivec[idx1] = phivec[0] + dphi / 2.0; /* calculate center point */
    if (!checkabsorb (phivec[idx1], &chan_num, nrad, tw, sw, bw, blader, off, c))
      ABSORB;
    chan_num0 = chan_num;
    while (idx1 < 129) {
      dphi = phivec[1] - phivec[idx1];
      idx1++;
      phivec[idx1] = phivec[0] + dphi / 2.0;
      if (!checkabsorb (phivec[idx1], &chan_num, nrad, tw, sw, bw, blader, off, c))
        ABSORB;
      if ((chan_num != chan_num0) || (chan_num > nchan))
        ABSORB;
      /*  If the current dphi is positive calculate points until a point is beyond phivec[1]
                     Check to see if the point is absorbed after each new point is generated stop if more than 129 iterations are performed
             */
      if (dphi > 0) {
        while ((phivec[idx1] < phivec[1]) && (idx1 < 129)) {
          /* printf("phivec[%i]: %g dphi: %g phivec[1]: %g\n", idx1,phivec[idx1],dphi,phivec[1]);*/
          idx1++;
          phivec[idx1] = phivec[idx1 - 1] + dphi;
          if (!checkabsorb (phivec[idx1], &chan_num, nrad, tw, sw, bw, blader, off, c))
            ABSORB;
          //   printf("chan_num0: %i chan_num: %i\n",chan_num0,chan_num);
          if ((chan_num != chan_num0) || (chan_num > nchan))
            ABSORB;
        }
        if (phivec[idx1] >= phivec[1])
          idx1--; // remove the point that is beyond phivec[1]
      }
      /*  If the current dphi is negative calculate points until a point is beyond phivec[1]
                     Check to see if the point is absorbed after each new point is generated stop if more than 129 iterations are performed
             */
      else if (dphi < 0) {
        while ((phivec[idx1] > phivec[1]) && (idx1 < 129)) {
          /* printf("phivec[%i]: %g\n", idx1,phivec[idx1]);*/
          idx1++;
          phivec[idx1] = phivec[idx1 - 1] + dphi;
          if (!checkabsorb (phivec[idx1], &chan_num, nrad, tw, sw, bw, blader, off, c))
            ABSORB;
          //   printf("chan_num0: %i chan_num: %i\n",chan_num0,chan_num);
          if ((chan_num != chan_num0) || (chan_num > nchan))
            ABSORB;
        }
        if (phivec[idx1] <= phivec[1])
          idx1--; // remove the point that is beyond phivec[1]
      } else
        ABSORB; /* dphi =0? */
    }
  } else /* The neutron failed to even hit the chopper */
    ABSORB;
%}
MCDISPLAY
%{
  double zstep, x1, x2, x3, x4, z1, z2;
  int idx, idx2;
  line (tw / 2.0, ymin, splen, tw / 2.0, ymax, splen);
  line (tw / 2.0, ymin, -splen, tw / 2.0, ymax, -splen);
  line (-tw / 2.0, ymin, splen, -tw / 2.0, ymax, splen);
  line (-tw / 2.0, ymin, -splen, -tw / 2.0, ymax, -splen);
  line (tw / 2.0, ymax, splen, tw / 2.0, ymax, -splen);
  line (-tw / 2.0, ymax, splen, -tw / 2.0, ymax, -splen);
  line (tw / 2.0, ymin, splen, tw / 2.0, ymin, -splen);
  line (-tw / 2.0, ymin, splen, -tw / 2.0, ymin, -splen);
  circle ("zx", 0, ymin, 0, rad);
  circle ("zx", 0, ymax, 0, rad);
  zstep = 2.0 * splen / 10.0;
  for (idx = 0; idx < nchan + 1; idx++) {
    for (idx2 = 0; idx2 < 10; idx2++) {
      z1 = idx2 * zstep - splen;
      z2 = (idx2 + 1) * zstep - splen;
      x1 = blades (z1, blader, off) + idx * (sw + bw) - tw / 2.0;
      x2 = blades (z2, blader, off) + idx * (sw + bw) - tw / 2.0;
      x3 = x1 + bw;
      x4 = x2 + bw;
      line (x1, ymin, z1, x2, ymin, z2);
      line (x1, ymax, z1, x2, ymax, z2);
      line (x3, ymin, z1, x4, ymin, z2);
      line (x3, ymax, z1, x4, ymax, z2);
      if (idx2 == 0) {
        line (x1, ymin, z1, x1, ymax, z1);
        line (x3, ymin, z1, x3, ymax, z1);
        line (x1, ymin, z1, x3, ymin, z1);
        line (x1, ymax, z1, x3, ymax, z1);
      }
      if (idx2 == 9) {
        line (x2, ymin, z2, x2, ymax, z2);
        line (x4, ymin, z2, x4, ymax, z2);
        line (x2, ymin, z2, x4, ymin, z2);
        line (x2, ymax, z2, x4, ymax, z2);
      }
    }
  }
%}
END