/*******************************************************************************
* McXtrace, X-ray tracing package
*           Copyright, All rights reserved
*           Risoe National Laboratory, Roskilde, Denmark
*           Institut Laue Langevin, Grenoble, France
*
* Component: SANSQMonitor
*
* %I
* Written by: Søren Kynde (kynde@nbi.dk) and Martin Cramer Pedersen (mcpe@nbi.dk)
* Date: October 17, 2012
* Origin: KU-Science
*
* A circular detector measuring the radial average of intensity as a function
* of the momentum transform in the sample.
*
* %D
* A simple component simulating the scattering from a box-shaped, thin solution
* of monodisperse, spherical particles.
*
* %P
* RadiusDetector: [m]	Radius of the detector (in the xy-plane).
* DistanceFromSample: [m]  Distance from the sample to this component.
* LambdaMin: [AA]	Max sensitivity in lambda - used to compute the highest possible value of momentum transfer, q.
* NumberOfBins: []	Number of bins in the r (and q).
* RFilename: []	File used for storing I(r).
* qFilename: [] 	File used for storing I(q).
* Lambda0: []	If lambda is given, the momentum transfers of all rays are computed from this value. Otherwise, instrumental effects are neglected.
* restore_neutron: []	If set to 1, the component restores the original neutron.
*
* %E
*******************************************************************************/

DEFINE COMPONENT SANSQMonitor



SETTING PARAMETERS (NumberOfBins = 100,
  string RFilename = "RDetector", string qFilename = "QDetector",
  RadiusDetector, DistanceFromSample, LambdaMin = 1.0, Lambda0 = 0.0, restore_neutron = 0)



DECLARE
%{
  // Declarations
  double TwoThetaMax;
  double qMax;

  DArray1d Nofq;
  DArray1d Iofq;
  DArray1d IofqSquared;

  DArray1d NofR;
  DArray1d IofR;
  DArray1d IofRSquared;
%}

INITIALIZE
%{
  // Declarations
  int i;

  Nofq = create_darr1d (NumberOfBins);
  Iofq = create_darr1d (NumberOfBins);
  IofqSquared = create_darr1d (NumberOfBins);

  NofR = create_darr1d (NumberOfBins);
  IofR = create_darr1d (NumberOfBins);
  IofRSquared = create_darr1d (NumberOfBins);

  // Initializations
  for (i = 0; i < NumberOfBins; ++i) {
    Nofq[i] = 0.0;
    Iofq[i] = 0.0;
    IofqSquared[i] = 0.0;
  }

  TwoThetaMax = atan (RadiusDetector / DistanceFromSample);
  qMax = 4 * PI * sin (TwoThetaMax / 2.0) / LambdaMin;
%}

TRACE
%{
  // Declarations
  int i;
  double TwoTheta;
  double Lambda;

  double R;
  double RLow;
  double RHigh;

  double q;
  double qLow;
  double qHigh;

  double TwoThetaLow;
  double TwoThetaHigh;
  double AreaOfSlice;

  PROP_Z0;

  // Computation of R
  R = sqrt (pow (x, 2) + pow (y, 2));

  // Computation of q
  if (Lambda0 <= 0.0) {
    Lambda = 2.0 * PI / (V2K * sqrt (pow (vx, 2) + pow (vy, 2) + pow (vz, 2)));
  } else {
    Lambda = Lambda0;
  }

  TwoTheta = atan (R / DistanceFromSample);
  q = 4.0 * PI * sin (TwoTheta / 2.0) / Lambda;

  // Put neutron in the correct r-bin
  if (R < RadiusDetector) {
    i = floor (NumberOfBins * R / RadiusDetector);

    RLow = RadiusDetector / NumberOfBins * i;
    RHigh = RadiusDetector / NumberOfBins * (i + 1);

    TwoThetaLow = atan (RLow / DistanceFromSample);
    TwoThetaHigh = atan (RHigh / DistanceFromSample);

    AreaOfSlice = fabs ((cos (2.0 * TwoThetaLow) - cos (2.0 * TwoThetaHigh)) * 2.0 * PI);
    #pragma acc atomic
    NofR[i] = NofR[i] + 1;
    double p_A = p / AreaOfSlice;
    #pragma acc atomic
    IofR[i] = IofR[i] + p_A;
    double p2_A2 = p * p / (AreaOfSlice * AreaOfSlice);
    #pragma acc atomic
    IofRSquared[i] = IofRSquared[i] + p2_A2;
  }

  // Put neutron in the correct q-bin
  if (q < qMax) {
    i = floor (NumberOfBins * q / qMax);

    qLow = qMax / NumberOfBins * i;
    qHigh = qMax / NumberOfBins * (i + 1);

    TwoThetaLow = asin (qLow * Lambda / (4.0 * PI));
    TwoThetaHigh = asin (qHigh * Lambda / (4.0 * PI));

    AreaOfSlice = fabs ((cos (2.0 * TwoThetaLow) - cos (2.0 * TwoThetaHigh)) * 2.0 * PI);

    #pragma acc atomic
    Nofq[i] = Nofq[i] + 1;
    double p_A = p / AreaOfSlice;
    #pragma acc atomic
    Iofq[i] = Iofq[i] + p_A;
    double p2_A2 = p * p / (AreaOfSlice * AreaOfSlice);
    #pragma acc atomic
    IofqSquared[i] = IofqSquared[i] + p2_A2;

    SCATTER;
  }

  // Restore neutron if requested
  if (restore_neutron) {
    //		RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
  }
%}

SAVE
%{
  // Output I(r)
  DETECTOR_OUT_1D ("QMonitor - Radially averaged distribution", "Radius [m]", "I(r)", "r", 0.0, RadiusDetector, NumberOfBins, &NofR[0], &IofR[0], &IofRSquared[0],
                   RFilename);

  // Output I(q)
  DETECTOR_OUT_1D ("QMonitor - Distribution in q (Radially averaged)", "q [1 / AA]", "I(q)", "q", 0.0, qMax, NumberOfBins, &Nofq[0], &Iofq[0], &IofqSquared[0],
                   qFilename);
%}

FINALLY
%{
  destroy_darr1d (Nofq);
  destroy_darr1d (Iofq);
  destroy_darr1d (IofqSquared);

  destroy_darr1d (NofR);
  destroy_darr1d (IofR);
  destroy_darr1d (IofRSquared);
%}

MCDISPLAY
%{
  circle ("xy", 0, 0, 0, RadiusDetector);
%}

END