/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Res_sample
*
* %I
* Written by: Kristian Nielsen
* Date: 1999
* Origin: Risoe
* Modified by: T. Weber, Nov 2020: updated for McStas 3 / OpenACC
*
* Sample for resolution function calculation.
*
* %D
* An inelastic sample with completely uniform scattering in both Q and
* energy. This sample is used together with the Res_monitor component and
* (optionally) the mcresplot front-end to compute the resolution function of
* triple-axis or inverse-geometry time-of-flight instruments.
*
* The shape of the sample is either a hollow cylinder or a rectangular box. The
* hollow cylinder shape is specified with an inner and outer radius.
* The box is specified with dimensions xwidth, yheight, zdepth.
*
* The scattered neutrons will have directions towards a given disk and
* energies betweed E0-dE and E0+dE.
* This target area may also be rectangular if specified focus_xw and focus_yh
* or focus_aw and focus_ah, respectively in meters and degrees.
* The target itself is either situated according to given coordinates (x,y,z), or
* setting the relative target_index of the component to focus at (next is +1).
* This target position will be set to its AT position. When targeting to centered
* components, such as spheres or cylinders, define an Arm component where to focus at.
*
* Example: Res_sample(thickness=0.001,radius=0.02,yheight=0.4,focus_r=0.05, E0=14.6,dE=2, target_x=0, target_y=0, target_z=1)
*
* %P
* INPUT PARAMETERS:
* thickness: [m]     Thickness of hollow cylinder in (x,z) plane 
* radius: [m]        Outer radius of hollow cylinder 
* yheight: [m]       vert. dimension of sample, as a height 
* focus_r: [m]       Radius of sphere containing target. 
* target_x: [m]       X-position of target to focus at [m]
* target_y: [m]      Y-position of target to focus at 
* target_z: [m]       Z-position of target to focus at [m]
* E0: [meV]          Center of scattered energy range 
* dE: [meV]          half width of scattered energy range 
*
* Optional parameters
* xwidth: [m]        horiz. dimension of sample, as a width 
* zdepth: [m]        depth of sample 
* focus_xw: [m]      horiz. dimension of a rectangular area 
* focus_yh: [m]      vert. dimension of a rectangular area 
* focus_aw: [deg]    horiz. angular dimension of a rectangular area 
* focus_ah: [deg]    vert. angular dimension of a rectangular area 
* target_index: [1]  relative index of component to focus at, e.g. next is +1 
*
* %E
*******************************************************************************/

DEFINE COMPONENT Res_sample

SETTING PARAMETERS (thickness=0, radius=0.01, focus_r=0.05, E0=14, dE=2,
  target_x=0, target_y=0, target_z=.5, focus_xw=0, focus_yh=0,
  focus_aw=0, focus_ah=0, xwidth=0, yheight=0.05, zdepth=0, int target_index=0)


SHARE %{
  struct res_sample_vars {
    char   isrect;                     /* true when sample is a box */
    double awdim, ahdim;               /* rectangular angular dimensions */
    double xwdim, yhdim;               /* rectangular metrical dimensions */
    double targetx, targety, targetz;  /* target coords */
  };
%}

/* Needed for resolution calculation */
USERVARS %{
  double res_pi;
  double res_ki_x;
  double res_ki_y;
  double res_ki_z;
  double res_kf_x;
  double res_kf_y;
  double res_kf_z;
  double res_rx;
  double res_ry;
  double res_rz;
%}

DECLARE %{
  struct res_sample_vars vars;
  char res_pi_var[20];
  char res_ki_x_var[20];
  char res_ki_y_var[20];
  char res_ki_z_var[20];
  char res_kf_x_var[20];
  char res_kf_y_var[20];
  char res_kf_z_var[20];
  char res_rx_var[20];
  char res_ry_var[20];
  char res_rz_var[20];
  int compindex;
%}


INITIALIZE %{
  if (!radius || !yheight) {
    if (!xwidth || !yheight || !zdepth)
      exit(fprintf(stderr,"Res_sample: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
    else
      vars.isrect = 1;
  }
  else {
    vars.isrect = 0;
  }

  /* now compute target coords if a component index is supplied */
  if (!target_index && !target_x && !target_y && !target_z)
    target_index = 1;
  if (target_index) {
    Coords ToTarget;
    ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
    ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
    coords_get(ToTarget, &vars.targetx, &vars.targety, &vars.targetz);
  }
  else {
    vars.targetx = target_x;
    vars.targety = target_y;
    vars.targetz = target_z;
  }

  if (!(vars.targetx || vars.targety || vars.targetz)) {
    printf("Res_sample: %s: The target is not defined. Using direct beam (Z-axis).\n",
      NAME_CURRENT_COMP);
    vars.targetz = 1;
  }

  /* different ways of setting rectangular area */
  vars.awdim = vars.ahdim = 0;
  if (focus_xw) vars.xwdim = focus_xw;
  if (focus_yh) vars.yhdim = focus_yh;
  if (focus_aw) vars.awdim = DEG2RAD*focus_aw;
  if (focus_ah) vars.ahdim = DEG2RAD*focus_ah;

  /* Initialize uservar strings */
  sprintf(res_pi_var,"res_pi_%i",_comp->_index);
  sprintf(res_ki_x_var,"res_ki_x_%i",_comp->_index);
  sprintf(res_ki_y_var,"res_ki_y_%i",_comp->_index);
  sprintf(res_ki_z_var,"res_ki_z_%i",_comp->_index);
  sprintf(res_kf_x_var,"res_kf_x_%i",_comp->_index);
  sprintf(res_kf_y_var,"res_kf_y_%i",_comp->_index);
  sprintf(res_kf_z_var,"res_kf_z_%i",_comp->_index);
  sprintf(res_rx_var,"res_rx_%i",_comp->_index);
  sprintf(res_ry_var,"res_ry_%i",_comp->_index);
  sprintf(res_rz_var,"res_rz_%i",_comp->_index);
  compindex=_comp->_index;
%}


TRACE %{
  double t0, t3;                /* Entry/exit time for outer cylinder */
  double t1, t2;                /* Entry/exit time for inner cylinder */
  double v;                     /* Neutron velocity */
  double E;
  double l_full;                /* Flight path length for non-scattered neutron */
  double dt0, dt1, dt2, dt;     /* Flight times through sample */
  double solid_angle=0;         /* Solid angle of target as seen from scattering point */
  double aim_x, aim_y, aim_z;   /* Position of target relative to scattering point */
  double scat_factor;           /* Simple cross-section model */
  int intersect = 0;
  double kix,kiy,kiz;
  double kfx,kfy,kfz;

  if(vars.isrect)
    intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
  else
    intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);

  if(intersect) {
    if(t0 < 0) ABSORB;
    if(vars.isrect) {
      t1 = t2 = t3;
      scat_factor = 2*zdepth;
    } /* box sample */
    else {  /* Hollow cylinder sample */
      /* Neutron enters at t=t0. */
      if(!thickness || !cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius-thickness, yheight))
        t1 = t2 = t3;
      scat_factor = 2*radius;
    }

    dt0 = t1-t0;                  /* Time in sample, ingoing */
    dt1 = t2-t1;                  /* Time in hole */
    dt2 = t3-t2;                  /* Time in sample, outgoing */

    v = sqrt(vx*vx + vy*vy + vz*vz);
    l_full = v * (dt0 + dt2);     /* Length of full path through sample */
    p *= l_full/scat_factor;      /* Scattering probability */
    dt = rand01()*(dt0+dt2);      /* Time of scattering (relative to t0) */
    if (dt > dt0)
      dt += dt1;

    PROP_DT(dt+t0);             /* Point of scattering */

    /* Store initial neutron state. */
    if(p == 0) ABSORB;
    kix=V2K*vx; kiy=V2K*vy; kiz=V2K*vz;
    particle_setvar_void(_particle, res_pi_var, &p);
    particle_setvar_void(_particle, res_ki_x_var, &(kix));
    particle_setvar_void(_particle, res_ki_y_var, &(kiy));
    particle_setvar_void(_particle, res_ki_z_var, &(kiz));
    particle_setvar_void(_particle, res_rx_var, &x);
    particle_setvar_void(_particle, res_ry_var, &y);
    particle_setvar_void(_particle, res_rz_var, &z);

    aim_x = vars.targetx - x;         /* Vector pointing at target (anal./det.) */
    aim_y = vars.targety - y;
    aim_z = vars.targetz - z;
 
   if(vars.awdim && vars.ahdim) {
      randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
        aim_x, aim_y, aim_z, vars.awdim, vars.ahdim, ROT_A_CURRENT_COMP);
    } else if(vars.xwdim && vars.yhdim) {
      randvec_target_rect(&vx, &vy, &vz, &solid_angle,
        aim_x, aim_y, aim_z, vars.xwdim, vars.yhdim, ROT_A_CURRENT_COMP);
    } else {
      randvec_target_circle(&vx, &vy, &vz, &solid_angle,
        aim_x, aim_y, aim_z, focus_r);
    }

    NORM(vx, vy, vz);
    E = E0 + dE*randpm1();
    v = sqrt(E)*SE2V;
    vx *= v;
    vy *= v;
    vz *= v;
    SCATTER;

     /* Store final neutron state. */
    kfx=V2K*vx; kfy=V2K*vy; kfz=V2K*vz;
    particle_setvar_void(_particle, res_kf_x_var, &(kfx));
    particle_setvar_void(_particle, res_kf_y_var, &(kfy));
    particle_setvar_void(_particle, res_kf_z_var, &(kfz));
  }
%}


MCDISPLAY %{
  if(vars.isrect) { /* Flat sample. */
    double xmin = -0.5*xwidth;
    double xmax =  0.5*xwidth;
    double ymin = -0.5*yheight;
    double ymax =  0.5*yheight;
    double len = zdepth/2;
    multiline(5, xmin, ymin, -len,
                 xmax, ymin, -len,
                 xmax, ymax, -len,
                 xmin, ymax, -len,
                 xmin, ymin, -len);
    multiline(5, xmin, ymin, len,
                 xmax, ymin, len,
                 xmax, ymax, len,
                 xmin, ymax, len,
                 xmin, ymin, len);
    line(xmin, ymin, -len, xmin, ymin, len);
    line(xmax, ymin, -len, xmax, ymin, len);
    line(xmin, ymax, -len, xmin, ymax, len);
    line(xmax, ymax, -len, xmax, ymax, len);
  }
  else {
    circle("xz", 0,  yheight/2.0, 0, radius);
    circle("xz", 0, -yheight/2.0, 0, radius);
    line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
    line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
    line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
    line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
    if (thickness) {
      double radius_i=radius-thickness;
      circle("xz", 0,  yheight/2.0, 0, radius_i);
      circle("xz", 0, -yheight/2.0, 0, radius_i);
      line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0);
      line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0);
      line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i);
      line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i);
    }
  }
%}

END