/*******************************************************************************
*
* McXtrace, x-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: PSD_monitor
*
* %Identification
* Written by: Erik B Knudsen
* Date: June 22, 2009
* Origin: Risoe
*
* Position-sensitive monitor.
*
* %Description
* Based on neutron component written by Kim Lefmann
* An (nx times ny) pixel PSD monitor. This component may also be used as a beam
* detector. If instead of xwidth, yheight a radius is given, the component has a circular footprint
* and integrates circularly (caking).
*
* Example: PSD_monitor(xwidth=0.1, yheight=0.1,
*          nx=90, ny=90, filename="Output.psd")
*
* %Parameters
* INPUT PARAMETERS:
*
* xwidth: [m]           Width of detector.
* yheight: [m]          Height of detector.
* radius:  [m]          Radius of circular detetor.
* nx: [1]               Number of pixel columns.
* ny: [1]               Number of pixel rows.
* nr: [1]               Number of radial pixels.
* filename: [str]       Name of file in which to store the detector image.
* restore_xray: [1]     If set, the monitor does not influence the xray state.
* nowritefile: [1]      If set, monitor will skip writing to disk
*
* %End
*******************************************************************************/

DEFINE COMPONENT PSD_monitor

SETTING PARAMETERS (string filename=0, xwidth=0.05, yheight=0.05, radius=0, restore_xray=1, 
  int nowritefile=0, int nx=90, int ny=90, int nr=0)

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 

DECLARE
%{
  DArray2d PSD_N;
  DArray2d PSD_p;
  DArray2d PSD_p2;
  double xmin;
  double xmax;
  double ymin;
  double ymax;
%}

INITIALIZE
%{
  int i, j;
  double *p1, *p2, *p3;

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

  if (((xmin >= xmax) || (ymin >= ymax)) && !radius) {
    fprintf (stderr, "ERROR (%s): Null detection area! Aborting.\n", NAME_CURRENT_COMP);
    exit (-1);
  }
  if (!radius) {
    PSD_N = create_darr2d (nx, ny);
    PSD_p = create_darr2d (nx, ny);
    PSD_p2 = create_darr2d (nx, ny);
  } else {
    PSD_N = create_darr2d (1, nr);
    PSD_p = create_darr2d (1, nr);
    PSD_p2 = create_darr2d (1, nr);
  }

  // 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, k;
  double e, p2;

  PROP_Z0;
  if (!radius) {
    if (x > xmin && x < xmax && y > ymin && y < ymax) {
      i = floor ((x - xmin) * nx / (xmax - xmin));
      j = floor ((y - ymin) * ny / (ymax - ymin));
      p2 = p * p;
      #pragma acc atomic
      PSD_N[i][j] += 1;
      #pragma acc atomic
      PSD_p[i][j] += p;
      #pragma acc atomic
      PSD_p2[i][j] += p2;
      SCATTER;
    }
  } else {
    double r = sqrt (x * x + y * y);
    if (r < radius) {
      i = floor (r * nr / radius);
      p2 = p * p;
      #pragma acc atomic
      PSD_N[0][i] += 1;
      #pragma acc atomic
      PSD_p[0][i] += p;
      #pragma acc atomic
      PSD_p2[0][i] += 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 (!radius) {
      if (nx == 1 && ny == 1) {
        char title[256];
        snprintf (title, 255, "Intensity monitor %s", NAME_CURRENT_COMP);
        DETECTOR_OUT_0D (title, (double)PSD_N[0][0], PSD_p[0][0], PSD_p2[0][0]);
      } else if (nx == 1) {
        DETECTOR_OUT_1D ("PSD_monitor", "Y Position[m]", "Intensity", "Y", ymin, ymax, ny, &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], filename);
      } else if (ny == 1) {
        DETECTOR_OUT_1D ("PSD_monitor", "X Position[m]", "Intensity", "X", xmin, xmax, nx, &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], filename);
      } else {
        DETECTOR_OUT_2D ("PSD monitor", "X position [m]", "Y position [m]", xmin, xmax, ymin, ymax, nx, ny, *PSD_N, *PSD_p, *PSD_p2, filename);
      }
    } else {
      DETECTOR_OUT_1D ("PSD_monitor", "Radial Position[m]", "Intensity", "R", 0, radius, nr, &PSD_N[0][0], &PSD_p[0][0], &PSD_p2[0][0], filename);
    }
  }
%}

FINALLY
%{
  destroy_darr2d (PSD_N);
  destroy_darr2d (PSD_p);
  destroy_darr2d (PSD_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