/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Collimator_radial
*
* %I
* Written by: Emmanuel Farhi <farhi@ill.fr>
* Date: July 2005
* Origin: ILL
*
* A radial Soller collimator.
*
* %D
* Radial Soller collimator with rectangular opening and specified length.
* The collimator is made of many rectangular channels stacked radially.
* Each channel is a set of transmitting layers (nslit), separated by an absorbing
* material (infinitely thin), the whole stuff is inside an absorbing housing.
*
* When specifying the number of channels (nchan), each channel has a total
* entrance width=radius*fabs(theta_max-theta_min)/nchan, but only the central
* portion 'xwidth' accepts neutrons. When xwidth=0, it is set to the full
* apperture so that all neutrons enter the channels (all walls are infinitely thin).
*
* When using zero as the number of channels (nchan), the collimator is continuous,
* whithout shadowing effect.
*
* The component should be positioned at the radius center.
* The component can be made oscillating (usual on diffractometers and TOF
* machines) with the 'roc' parameter.
* The neutron beam outside the collimator angular area is transmitted unaffected.
*
* When used as a focusing collimator, the focusing parameter should be set to 1.
*
* An example of a instrument that uses this collimator can be found in the SALSA instrument,
* in the example folder
*
* Example:
*   Channelled radial collimator with shadow parts
*     Collimator_radial(xwidth=0.015, yheight=.3, length=.35, divergence=40,transmission=1, theta_min=5, theta_max=165, nchan=128, radius=0.9)
*   A continuous radial collimator
*     Collimator_radial(yheight=.3, length=.35, divergence=40,transmission=1, theta_min=5, theta_max=165, radius=0.9)
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m]               Soller  window width, filled with nslit slits. Use 0 value for continuous collimator.
* yheight: [m]              Collimator height. If  yheight_inner is specified, then this is the outer cylinders height
* length: [m]               Length/Distance between inner and outer slits.
* divergence: [min of arc]  Divergence angle. May also be specified with the nslit parameter. A zero value unactivates component.
* theta_min: [deg]          Minimum Theta angle for the radial setting.
* theta_max: [deg]          Maximum Theta angle for the radial setting.
* nchan: [1]                Number of Soller channels in the theta range. Use 0 value for continuous collimator.
* radius: [m]               Radius of the collimator (to entry window).
*
* Optional parameters
* transmission: [1]         Maximum transmission of Soller (0<=t<=1).
* nslit: [1]                Number of blades composing each Soller. Overrides the divergence parameter.
* roc: [deg]                Amplitude of oscillation of collimator. 0=fixed.
* verbose:     []          Gives additional information.
* approx:      []          Use Soller triangular transmission approximation.
* focusing: [1]             When set allows you to use the collimators for focusing, rather than dispersing.
* yheight_inner: [1]         Defines the inner height of the collimator
*
* %E
*******************************************************************************/
DEFINE COMPONENT Collimator_radial

SETTING PARAMETERS (xwidth=0, yheight=.3, length=.35,
divergence=0,transmission=1,
theta_min=5, theta_max=165, nchan=0, radius=1.3, nslit=0,
roc=0, verbose=0, approx=0, focusing = 0, yheight_inner = 0)

DECLARE
%{
  double width_of_slit;
  double width_of_Soller;
  double slit_theta;
%}

INITIALIZE
%{
  width_of_slit = 0;
  width_of_Soller = 0;
  slit_theta = 0;

  if (radius <= 0)
    exit (printf ("Collimator_radial: %s: incorrect radius=%g\n", NAME_CURRENT_COMP, radius));
  if (length <= 0)
    exit (printf ("Collimator_radial: %s: invalid collimator length=%g\n", NAME_CURRENT_COMP, length));
  if (transmission <= 0 || transmission > 1)
    exit (printf ("Collimator_radial: %s: invalid transmission=%g\n", NAME_CURRENT_COMP, transmission));

  theta_max *= DEG2RAD;
  theta_min *= DEG2RAD;
  roc *= DEG2RAD;
  divergence *= MIN2RAD;

  if (xwidth && !nchan)
    nchan = ceil (radius * fabs (theta_max - theta_min) / xwidth);
  else if (!xwidth && nchan)
    xwidth = radius * fabs (theta_max - theta_min) / nchan;

  /* determine total width [m] of Soller channels, containing nslit in xwidth */
  if (nchan) {
    width_of_Soller = radius * fabs (theta_max - theta_min) / nchan;
  } else
    width_of_Soller = 0;

  if (!nchan || !xwidth || xwidth > width_of_Soller)
    nchan = xwidth = width_of_Soller = 0; /* continuous collimator */

  /* determine width [m] of slits */
  if (divergence) {
    width_of_slit = length * tan (divergence);
    if (xwidth) /* Soller */
      nslit = ceil (xwidth / width_of_slit);
    else if (!nchan) /* continuous collimator */
      nslit = ceil (radius * fabs (theta_max - theta_min) / width_of_slit);
  } else {
    if (!nchan && nslit) /* continuous collimator */
      width_of_slit = radius * fabs (theta_max - theta_min) / nslit;
    else if (nchan && nslit) /* Soller */
      width_of_slit = xwidth / nslit;
    divergence = atan2 (width_of_slit, length);
  }

  if (nslit <= 0)
    printf ("Collimator_radial: %s: number of channels must be positive nslit=%g.\n"
            "WARNING            Specify alternatively divergence, xwidth and nchan.\n",
            NAME_CURRENT_COMP, nslit);

  if (verbose && nslit && width_of_slit) {
    printf ("Collimator_radial: %s: divergence %g [min] %s"
            ". Total opening [%g:%g] [deg]\n",
            NAME_CURRENT_COMP, divergence * RAD2MIN, (roc ? "oscillating" : ""), theta_min * RAD2DEG, theta_max * RAD2DEG);
    if (approx)
      printf ("    Using triangular approximation model");
    else if (!nchan)
      printf ("    Using continuous model");
    else
      printf ("    Using %i Soller channels of width %g [cm]", (int)nchan, width_of_Soller * 100);

    printf (" with %i slits of width %g [mm] pitch %g [deg].\n", (int)nslit, width_of_slit * 1000, atan2 (width_of_slit, radius) * RAD2DEG);
  }
  if (!yheight_inner)
    yheight_inner = yheight;
%}

TRACE
%{
  double intersect = 0;
  double t0, t1, t2, t3;

  if (width_of_slit && nslit) {
    /* determine intersection with inner and outer cylinders */
    intersect = cylinder_intersect (&t0, &t3, x, y, z, vx, vy, vz, radius, yheight_inner);
    if (!intersect)
      ABSORB;
    else if (t3 > t0 && !focusing)
      t0 = t3;

    intersect = cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius + length, yheight);
    if (!intersect)
      ABSORB;
    else if (t2 > t1 && !focusing)
      t1 = t2;

    /* propagate/determine neutron position at ingoing cylinder */
    if ((t0 > 0 && t1 > t0 && !focusing) || (t1 > 0 && t0 > t1 && focusing)) {
      double input_chan = 0;
      double input_theta = 0, output_theta = 0;
      double roc_theta = 0;
      double input_slit = 0, output_slit = 0;
      if (!focusing)
        PROP_DT (t0);
      if (focusing)
        PROP_DT (t1);

      /* apply ROC oscillation with a linear random distribution */
      if (roc)
        roc_theta = roc * randpm1 () / 2;
      else
        roc_theta = 0;

      /* angle on the initial relevant cylinder */
      input_theta = atan2 (x, z) + roc_theta;

      /* check if we are within min/max collimator input bounds */
      if (input_theta >= theta_min && input_theta <= theta_max) {
        SCATTER;
        /* input Soller channel index */
        if (width_of_Soller) {
          input_chan = radius * (input_theta - theta_min) / width_of_Soller;
          input_chan = input_chan - floor (input_chan); /* position within Soller [0:1] */

          /* check if we hit an absorbing housing (between Sollers): ABSORB
           *   total Soller aperture is width_of_Soller,
           *   containg a slit pack  of aperture xwidth
           */
          if (input_chan < (1 - xwidth / width_of_Soller) / 2 || input_chan > (1 + xwidth / width_of_Soller) / 2)
            ABSORB;
        }

        /* determine input slit index */
        input_slit = floor (input_theta * radius / width_of_slit);

      } else /* neutron missed collimator input range */
        input_theta = 4 * PI;

      /* propagate to next cylinder */

      if (!focusing)
        PROP_DT (t1 - t0);
      if (focusing)
        PROP_DT (t0 - t1);

      /* angle on the outgoing cylinder */
      output_theta = atan2 (x, z) + roc_theta;

      /* check if we are within min/max collimator output bounds */
      if (output_theta >= theta_min && output_theta <= theta_max) {
        /* check if we come from sides:   ABSORB */
        if (input_theta > 2 * PI)
          ABSORB; /* input_theta=4*PI when missed input */

        if (approx) {
          double phi = atan2 (x, z) - atan2 (vx, vz); /* difference between positional slit angle and velocity */
          if (fabs (phi) > divergence)
            ABSORB; /* get outside transmission */
          else
            p *= (1.0 - phi / divergence);
        } else {
          /* check if we have changed slit: ABSORB */
          /* slits are considered radial so that their output size is:
             width_of_slit*(radius+length)/radius and it turns out this the same exp as for input
           */
          output_slit = floor (output_theta * radius / width_of_slit);
          if (input_slit != output_slit)
            ABSORB;
        }
        SCATTER;
        p *= transmission;

      } else if (focusing)
        ABSORB;
      /* else neutron missed collimator output range*/

    } /* else did not encounter collimator cylinders */
  } /* if nslit */
%}

MCDISPLAY
%{
  double Soller_theta;
  double height_inner = yheight / 2;
  double height_outer = yheight_inner / 2;
  double theta1, theta2;
  double x_in_l, z_in_l, x_in_r, z_in_r;
  double x_out_l, z_out_l, x_out_r, z_out_r;
  int i;

  /* display collimator radial geometry:
     in order to avoid too many lines, we shown main housing and channels
     but no slit */

  if (!nchan || nchan > 20)
    nchan = 20;
  if (nchan > 64)
    nchan = 64;
  Soller_theta = fabs (theta_max - theta_min) / nchan; /* angular width of Soller */

  /* draw all channels, which also show housing */

  for (i = 0; i < nchan; i++) {

    theta1 = i * Soller_theta + theta_min;
    theta2 = theta1 + Soller_theta;

    z_in_l = radius * cos (theta1);
    x_in_l = radius * sin (theta1);
    z_in_r = radius * cos (theta2);
    x_in_r = radius * sin (theta2);

    z_out_l = (radius + length) * cos (theta1);
    x_out_l = (radius + length) * sin (theta1);
    z_out_r = (radius + length) * cos (theta2);
    x_out_r = (radius + length) * sin (theta2);
    /* left -> right lines */
    line (x_in_l, height_outer, z_in_l, x_in_r, height_outer, z_in_r);
    line (x_out_l, height_inner, z_out_l, x_out_r, height_inner, z_out_r);
    line (x_out_l, -height_inner, z_out_l, x_out_r, -height_inner, z_out_r);
    line (x_in_l, -height_outer, z_in_l, x_in_r, -height_outer, z_in_r);
  }
  /* outer walls, at theta_min and theta_max */
  for (i = 0; i <= nchan; i += nchan) {
    theta1 = i * Soller_theta + theta_min;
    z_in_l = radius * cos (theta1);
    x_in_l = radius * sin (theta1);
    z_out_l = (radius + length) * cos (theta1);
    x_out_l = (radius + length) * sin (theta1);
    multiline (5, x_in_l, -height_outer, z_in_l, x_in_l, height_outer, z_in_l, x_out_l, height_inner, z_out_l, x_out_l, -height_inner, z_out_l, x_in_l,
               -height_outer, z_in_l);
  }
  /* absorbing channel walls in dark grey */
  mcdis_pencolor ("#404040");
  for (i = 0; i <= nchan; i++) {
    theta1 = i * Soller_theta + theta_min;
    z_in_l = radius * cos (theta1);
    x_in_l = radius * sin (theta1);
    z_out_l = (radius + length) * cos (theta1);
    x_out_l = (radius + length) * sin (theta1);
    polygon (4, x_in_l, -height_outer, z_in_l, x_in_l, height_outer, z_in_l, x_out_l, height_inner, z_out_l, x_out_l, -height_inner, z_out_l);
  }
%}

END