/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2009, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: StatisticalChopper_Monitor
*
* %Identification
* Written by: C. Monzat/E. Farhi
* Date: 2009
* Origin: ILL
*
* Monitor designed to compute the autocorrelation signal for the Statistical Chopper
*
* %Description
* This component is a time sensitive monitor which calculates the cross
* correlation between the pseudo random sequence of a statistical chopper
* and the signal received. It mainly uses fonctions of component Monitor_nD.
* It is possible to use the various options of the Monitor_nD but the user MUST NOT
* specify "time" in the options. Auto detection of the time limits is possible if the
* user chooses tmin=>tmax.
*
* StatisticalChopper_Monitor(options ="banana bins=500, abs theta limits=[5,105],bins=1000")
*
* %Parameters
* INPUT PARAMETERS:
* comp: [StatisticalChopper]  quoted name of the component monitored
* tmin: [s]                   minimal time of detection
* tmax: [s]                   maximal time of detection
*
* CALCULATED PARAMETERS:
* delta_t: [s]                interval of time of the detection
* T_p: [s]                    period of the statistical chopper comp
*
* %L
* R. Von Jan and R. Scherm. The statistical chopper for neutron time-of-flight spectroscopy. Nuclear Instruments and Methods, 80 (1970) 69-76.
*
* %End
*******************************************************************************/
DEFINE COMPONENT StatisticalChopper_Monitor INHERIT Monitor_nD

SETTING PARAMETERS (string comp, tmin=0.0, tmax=0.0)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

DECLARE INHERIT Monitor_nD

INITIALIZE INHERIT Monitor_nD EXTEND
%{
  char op[CHAR_BUF_LENGTH];
  /* check if time is already in the variables to monitor */
  if (!strstr (Vars.option, "time")) {
    sprintf (op, "time %s", Vars.option);
    strcpy (Vars.option, op);
  } else {
    /* we must make sure time is the first variable */
    if (Vars.Coord_Number < 1 || strcmp (Vars.Coord_Var[1], "t"))
      exit (fprintf (stderr, "StatisticalChopper_Monitor: %s: First variable must be 'time', currently '%s (%s)'\n.", NAME_CURRENT_COMP, Vars.Coord_Label[1],
                     Vars.Coord_Var[1]));
    if (strcmp (Vars.Coord_Var[0], "I") && strcmp (Vars.Coord_Var[0], "p"))
      exit (fprintf (stderr, "StatisticalChopper_Monitor: %s: Must record signal intensity 'Intensity', currently '%s (%s)'\n.", NAME_CURRENT_COMP,
                     Vars.Coord_Label[0], Vars.Coord_Var[0]));
  }

  /* set up Monitor_nD time limits */
  if (tmin < tmax) {
    sprintf (op, "limits=[%g %g] %s", tmin, tmax, Vars.option);
  } else {
    sprintf (op, "auto %s", Vars.option);
  }
  strcpy (Vars.option, op);

  /* re-initialize Monitor_nD */
  Monitor_nD_Init (&DEFS, &Vars, xwidth, yheight, zdepth, xmin, xmax, ymin, ymax, zmin, zmax, 0, 0);
%}

TRACE INHERIT Monitor_nD

SAVE
%{
  // avoid name clash with Detector structure members and component parameters.
  #undef min
  #undef xmin
  #undef xmax
  #undef ymin
  #undef ymax

  /* the detector file written by the Monitor_nD is stored in a 'MCDETECTOR detector' structure */

  /* get back information from the StatisticalChopper */
  int m = *COMP_GETPAR3 (StatisticalChopper, comp, m);         /* number of appertures in the sequence */
  int nslit = *COMP_GETPAR3 (StatisticalChopper, comp, nslit); /* length of the sequence (number of possible slits around disk) */
  double nu = *COMP_GETPAR3 (StatisticalChopper, comp, nu);
  int* Sequence = *COMP_GETPAR3 (StatisticalChopper, comp, Sequence);

  double c = (double)(m - 1) / (double)(nslit - 1); /* duty cycle */
  int f;                                            /* number of periods in the detected time range */

  /* save results, but do not free pointers */
  detector = Monitor_nD_Save (&DEFS, &Vars);

  f = (int)floor ((detector.xmax - detector.xmin) * nu);
  if (f < 1)
    f = 1;

  if (detector.intensity && m >= 1 && f >= 1 && nslit > 1 && detector.m > 1 && c < 1 && nu > 0 && Sequence) {

    double S_sum = 0;
    int i, j, k; /* indices for loops */

    /* copy raw detector as correlation base */
    long correlation_m = detector.m / f; /* new time binning for autocorrelation */

    double* p0 = malloc (correlation_m * detector.n * detector.p * sizeof (double)); /* Arrays to store correlation monitor */
    double* p1 = malloc (correlation_m * detector.n * detector.p * sizeof (double));
    double* p2 = malloc (correlation_m * detector.n * detector.p * sizeof (double));
    if (p0 && p1 && p2) {

      /* initialize arrays to zero */
      for (i = 0; i < correlation_m * detector.n * detector.p; i++)
        p0[i] = p1[i] = p2[i] = 0;

      /* indices in loops for detector */
      /* p1[i][j] = p1[index] with index= (!detector.istransposed ? i*n*p + j : i+j*m); */

      /* compute scattering function S: Von Jan and Scherm, Eq (18) */
      for (i = 0; i < correlation_m; i++) {
        for (k = 0; k < detector.n * detector.p; k++) {
          long index_new = !detector.istransposed ? i * detector.n * detector.p + k : i + k * correlation_m; /* [i][k] index in correlation */
          for (j = 0; j < f * nslit; j++) {
            long i2 = (j * (detector.m - 1) / (f * nslit - 1) + i * (detector.m - 1) / (correlation_m - 1)) % detector.m; /* time index in raw detector */
            long index_old = !detector.istransposed ? i2 * detector.n * detector.p + k : i2 + k * detector.m;             /* full [i][k] index in raw detector */

            p1[index_new] += (Sequence[j % nslit] - c) * detector.p1[index_old];
          }
          p1[index_new] /= m * f * (1 - c);
          p1[index_new] += -detector.min / m;
          p0[index_new] = detector.events / correlation_m;
          S_sum += p1[index_new];
        }
      }

      /* Normalization of p1 compared to detector.p1: must have same integrated intensity */
      /* normalizations:
         coeff_ZS=(S_sum/detector.intensity)*f;
         coeff_ZS=(S_sum/(m*detector.intensity-detector.min/nu))*(m*f*(1-m/nslit)+m);
      */
      double coeff_ZS = (S_sum / (m * detector.intensity - detector.min / nu)) * (m * f * (1.0 - (double)m / nslit) + m);
      for (i = 0; i < correlation_m * detector.n * detector.p; i++) {
        p1[i] /= coeff_ZS;
      }

      /* Calculation of delta(S)^2: Von Jan and Scherm, Eq (19) */
      for (i = 0; i < correlation_m; i++) {
        for (k = 0; k < detector.n * detector.p; k++) {
          long index_new = !detector.istransposed ? i * detector.n * detector.p + k : i + k * correlation_m; /* [i][k] index in correlation */
          for (j = 0; j < f * nslit; j++) {
            long i2 = (j * (detector.m - 1) / (f * nslit - 1) + i * (detector.m - 1) / (correlation_m - 1)) % detector.m; /* time index in raw detector */
            long index_old = !detector.istransposed ? i2 * detector.n * detector.p + k : i2 + k * detector.m;             /* full [i][k] index in raw detector */
            p2[index_new] += (Sequence[j % nslit] - c) * (Sequence[j % nslit] - c) * detector.p1[index_old] * detector.p1[index_old];
          }
        }
      }
      /* Transformation from sigma to p2 before output */
      for (i = 0; i < correlation_m * detector.n * detector.p; i++) {
        p2[i] /= (double)(m * m * f * f * (1 - c) * (1 - c));
        if (p2[i] > 0)
          p2[i] = (p0[i] > 1 ? ((p0[i] - 1) * p2[i] * p2[i] + p1[i] * p1[i] / p0[i]) / p0[i] : p1[i]);
      }

      char file[CHAR_BUF_LENGTH];
      sprintf (file, "%s_corr", filename ? filename : NAME_CURRENT_COMP);

      if (detector.rank == 1)
        DETECTOR_OUT_1D ("Correlation monitor", detector.xlabel, detector.ylabel, detector.xvar, detector.xmin, detector.xmax, correlation_m, p0, p1, p2, file);
      else if (detector.rank == 2)
        DETECTOR_OUT_2D ("Correlation monitor", detector.xlabel, detector.ylabel, detector.xmin, detector.xmax, detector.ymin, detector.ymax, correlation_m,
                         detector.n, p0, p1, p2, file);
      free (p0);
      free (p1);
      free (p2);
    } else {
      fprintf (stderr, "WARNING: StatisticalChopper_monitor %s: Could not allocate arrays for monitor output. Skipped writing...\n", NAME_CURRENT_COMP);
    }
  }
%}

FINALLY INHERIT Monitor_nD

MCDISPLAY INHERIT Monitor_nD

END