/*******************************************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Divergence_monitor
*
* %Identification
* Written by: Erik B Knudsen
* Based on neutron component by Kim Lefmann
* Date: Jun. '16
* Version: $Revision$
* Origin: DTU Physics
*
* Horizontal+vertical divergence monitor.
*
* %Description
* A 2D divergence sensitive monitor. The counts are distributed in
* (n times m) pixels.
*
* Example: Divergence_monitor(nh=20, nv=20, filename="Output.pos",
*           xwidth=0.1, yheight=0.1,
*           maxdiv_h=2, maxdiv_v=2)
*
* %Parameters
* INPUT PARAMETERS:
*
* xwidth: [m]           Width of detector. 
* yheight: [m]          Height of detector.
* nv: [1]               Number of pixel columns 
* nh: [1]               Number of pixel rows 
* nx: [1]               Vector definition of "forward" direction wrt. divergence, to be used e.g. when the monitor is rotated into the horizontal plane 
* ny: [1]               Vector definition of "forward" direction wrt. divergence, to be used e.g. when the monitor is rotated into the horizontal plane 
* nz: [1]               Vector definition of "forward" direction wrt. divergence, to be used e.g. when the monitor is rotated into the horizontal plane 
* maxdiv_v: [degrees]   Maximal vertical divergence detected 
* maxdiv_h: [degrees]   Maximal horizontal divergence detected 
* filename: [str]          Name of file in which to store the detector image text
* restore_xray: [1]     If set, the monitor does not influence the photon state 
* rad:  [1]             If set - divergence will be measured in radians.
* nowritefile: [1]      If set, monitor will skip writing to disk
*
* CALCULATED PARAMETERS:
*
* Div_N:    Array of photon ray counts
* Div_p:    Array of photon weight counts
* Div_p2:   Array of second moments
*
* %End
*******************************************************************************/
DEFINE COMPONENT Divergence_monitor

SETTING PARAMETERS (int nh=20, int nv=20, int rad=0, string filename=0, 
    xwidth=0.1, yheight=0.1, maxdiv_h=1, maxdiv_v=1, restore_xray=0, nx=0, ny=0, nz=1, int nowritefile=0)

    /* Xray  parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 
DECLARE
  %{
  DArray2d Div_N;
  DArray2d Div_p;
  DArray2d Div_p2;
  double xmin;
  double xmax;
  double ymin;
  double ymax;
  %}
INITIALIZE
  %{
  int i, j;

  xmax = xwidth / 2;
  xmin = -xmax;
  ymax = yheight / 2;
  ymin = -ymax;

  if ((xmin >= xmax) || (ymin >= ymax)) {
    printf ("ERROR: (%s): Null detection area! Exiting.\n", NAME_CURRENT_COMP);
    exit (-1);
  }

  Div_N = create_darr2d (nh, nv);
  Div_p = create_darr2d (nh, nv);
  Div_p2 = create_darr2d (nh, nv);

  NORM (nx, ny, nz);

  // Use instance name for monitor output if no input was given
  if (!strcmp (filename, "\0"))
    sprintf (filename, "%s", NAME_CURRENT_COMP);
  %}
TRACE
  %{
  int i, j;
  double h_div, v_div;
  double k, kn;

  PROP_Z0;
  if (x > xmin && x < xmax && y > ymin && y < ymax) {
    /* Find length of projection onto the [nx ny nz] axis */
    kn = scalar_prod (kx, ky, kz, nx, ny, nz);
    if (rad) {
      h_div = atan2 (kx, kn);
      v_div = atan2 (ky, kn);
    } else {
      h_div = RAD2DEG * atan2 (kx, kn);
      v_div = RAD2DEG * atan2 (ky, kn);
    }
    if (h_div < maxdiv_h && h_div > -maxdiv_h && v_div < maxdiv_v && v_div > -maxdiv_v) {
      i = floor ((h_div + maxdiv_h) * nh / (2.0 * maxdiv_h));
      j = floor ((v_div + maxdiv_v) * nv / (2.0 * maxdiv_v));
      double p2 = p * p;
      #pragma acc atomic
      Div_N[i][j] = Div_N[i][j] + 1;
      #pragma acc atomic
      Div_p[i][j] = Div_p[i][j] + p;
      #pragma acc atomic
      Div_p2[i][j] = Div_p2[i][j] + p2;
      SCATTER;
    }
  }
  if (restore_xray) {
    RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p);
  }
  %}
SAVE
  %{
  if (!nowritefile) {
    if (rad) {
      DETECTOR_OUT_2D ("Divergence monitor", "X divergence [rad]", "Y divergence [rad]", -maxdiv_h, maxdiv_h, -maxdiv_v, maxdiv_v, nh, nv, &Div_N[0][0],
                       &Div_p[0][0], &Div_p2[0][0], filename);
    } else {
      DETECTOR_OUT_2D ("Divergence monitor", "X divergence [deg]", "Y divergence [deg]", -maxdiv_h, maxdiv_h, -maxdiv_v, maxdiv_v, nh, nv, &Div_N[0][0],
                       &Div_p[0][0], &Div_p2[0][0], filename);
    }
  }
  %}
FINALLY
  %{
  destroy_darr2d (Div_N);
  destroy_darr2d (Div_p);
  destroy_darr2d (Div_p2);
  %}
MCDISPLAY
  %{

  multiline (5, (double)xmin, (double)ymin, 0.0, (double)xmax, (double)ymin, 0.0, (double)xmax, (double)ymax, 0.0, (double)xmin, (double)ymax, 0.0, (double)xmin,
             (double)ymin, 0.0);
%}

END