/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* %I
* Written by: Peter Willendrup, derived from TOF_lambda_monitor.comp
* Date: May 23, 2012
* Origin: DTU Physics
*
* Special "Brilliance" monitor.
*
* %D
* If used in the right setting, will output "instantaneous" and "mean" brilliances in units of Neutrons/cm^2/ster/AA/s. Conditions for proper units:
* <ul>
* <li>Use a with a source of area 1x1cm
* <li>The source must illuminate/focus to an area of 1x1cm a 1m distance
* <li>Parametrise the Brilliance_monitor with the frequency of the source
* <li>To not change the source TOF distribution, place the Brilliance monitor close to the source!
* </ul>
*
* with a source of area 1x1cm illuminating/focusing to an area of 1x1cm a 1m distance, this monitor will output "instantaneous" and "mean" brilliances in units of Neutrons/cm^2/ster/AA/s
*
* Here is an example of the use of the component. Note how the mentioned Unit conditions are implemented in instrument code.
*
*COMPONENT Source = ESS_moderator_long(
*    l_low = lambdamin, l_high = lambdamax, dist = 1, xw = 0.01, yh = 0.01,
*    freq = 14, T=50, tau=287e-6, tau1=0, tau2=20e-6,
*    n=20, n2=5, d=0.00286, chi2=0.9, I0=6.9e11, I2=27.6e10,
*    branch1=0, branch2=0.5, twopulses=0, size=0.01)
*  AT (0, 0, 0) RELATIVE Origin
*
*COMPONENT BRIL = Brilliance_monitor(nlam=196,nt=401,filename="bril.sim",
*	t_0=0,t_1=4000,lambda_0=lambdamin,
*	lambda_1=lambdamax, Freq=14)
*AT (0,0,0.000001) RELATIVE Source
*
* %P
* INPUT PARAMETERS:
*
* nlam: [1]             Number of bins in wavelength
* nt: [1]               Number of bins in TOF
* t_0: [us]             Minimum time
* t_1: [us]             Maximum time
* lambda_0: [AA]        Minimum wavelength detected
* lambda_1: [AA]        Maximum wavelength detected
* filename: [string]    Defines filenames for the detector images. Stored as:<br>Peak_&lt;filename&gt; and Mean_&lt;filename&gt;
* restore_neutron: [1]  If set, the monitor does not influence the neutron state
* Freq: [Hz]            Source frequency. Use freq=1 for reactor source
* srcarea:    Source area [cm^2]
* tofcuts: [1]          Flag to generate TOF-distributions as function of wavelength
* toflambda: [1]        Flag to generate TOF-lambda distribution output
* source_dist: [m]      Distance from source. Beware when transporting through neutron optics!
* xwidth: [m]           width of monitor
* yheight: [m]          height of monitor
* nowritefile: [1]      If set, monitor will skip writing to disk
*
* CALCULATED PARAMETERS:
*
* Div_N: []             Array of neutron counts
* Div_p: []             Array of neutron weight counts
* Div_p2: []            Array of second moments
*
* %E
*******************************************************************************/

DEFINE COMPONENT Brilliance_monitor



SETTING PARAMETERS (int nlam=101, int nt=1001, int nowritefile=0,
  lambda_0=0, lambda_1=20, int restore_neutron=0, Freq, int tofcuts=0,
  int toflambda=0, xwidth = 0.01, yheight=0.01, source_dist=1, string filename=0,
  t_0=0, t_1=20000, srcarea=1)



DECLARE
%{
  DArray2d BRIL_N;
  DArray2d BRIL_p;
  DArray2d BRIL_p2;

  DArray1d BRIL_mean;
  DArray1d BRIL_meanN;
  DArray1d BRIL_meanE;
  DArray1d BRIL_peak;
  DArray1d BRIL_peakN;
  DArray1d BRIL_peakE;

  DArray1d BRIL_shape;
  DArray1d BRIL_shapeN;
  DArray1d BRIL_shapeE;

  double tt_0;
  double tt_1;
  double xmin;
  double xmax;
  double ymin;
  double ymax;
  double ster;
  double prsec;
  double dlam;
  double dt;
%}

INITIALIZE
%{
  prsec = 1e-6;

  BRIL_N = create_darr2d(nt, nlam);
  BRIL_p = create_darr2d(nt, nlam);
  BRIL_p2 = create_darr2d(nt, nlam);
  BRIL_mean = create_darr1d(nlam);
  BRIL_meanN = create_darr1d(nlam);
  BRIL_meanE = create_darr1d(nlam);
  BRIL_peak = create_darr1d(nlam);
  BRIL_peakN = create_darr1d(nlam);
  BRIL_peakE = create_darr1d(nlam);

  BRIL_shape = create_darr1d(nt);
  BRIL_shapeN = create_darr1d(nt);
  BRIL_shapeE = create_darr1d(nt);


  tt_0 = t_0*prsec;
  tt_1 = t_1*prsec;
  dt = (t_1-t_0)*prsec/nt;
  dlam = (lambda_1-lambda_0)/(nlam-1);

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

  ster = xwidth * yheight/(source_dist*source_dist);

  // 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 div;
  double lambda;
  double Pnorm;

  PROP_Z0;
  lambda = (2*PI/V2K)/sqrt(vx*vx + vy*vy + vz*vz);
  if (x>xmin && x<xmax && y>ymin && y<ymax &&
      lambda > lambda_0 && lambda < lambda_1) {
    if (t < tt_1 && t > tt_0) {
      i = floor((lambda - lambda_0)*nlam/(lambda_1 - lambda_0));
      j = floor((t-tt_0)*nt/(tt_1-tt_0));
      Pnorm=p/dlam/ster/srcarea;
      double Pnorm2 = Pnorm*Pnorm;
      #pragma acc atomic
      BRIL_meanN[i] = BRIL_meanN[i]+1;
      #pragma acc atomic
      BRIL_mean[i] = BRIL_mean[i]+Pnorm;
      #pragma acc atomic      
      BRIL_meanE[i] = BRIL_meanE[i]+Pnorm2;
      
      Pnorm=Pnorm/Freq/dt;
      Pnorm2=Pnorm+Pnorm;
      #pragma acc atomic      
      BRIL_N[j][i] = BRIL_N[j][i]+1;
      #pragma acc atomic
      BRIL_p[j][i] = BRIL_p[j][i]+Pnorm;
      #pragma acc atomic
      BRIL_p2[j][i] = BRIL_p2[j][i]+Pnorm2;
    }
  }
  if (restore_neutron) {
    RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
  }
%}

SAVE
%{
if (!nowritefile) {
  /* First, dump the 2D monitor */

  /* For each Wavelength channel, find peak brilliance */
  int i,j,jmax;
  double Pnorm;
  char ff[256];
  char tt[256];

  for (i=0; i<nlam; i++) {
    Pnorm = -1;
    jmax = -1;
    for (j=0; j<nt; j++) {
  	  if (BRIL_p[j][i]>=Pnorm) {
	      Pnorm = BRIL_p[j][i];
	      jmax=j;
	    }
  	  BRIL_shape[j] = BRIL_p[j][i];
  	  BRIL_shapeN[j] = BRIL_N[j][i];
  	  BRIL_shapeE[j] = BRIL_p2[j][i];
    }
    if (tofcuts == 1) {
    	sprintf(ff, "Shape_%s_%g",filename,lambda_0+i*dlam);
    	sprintf(tt, "Peak shape at %g AA",lambda_0+i*dlam);
    	DETECTOR_OUT_1D(
  			tt,
  			"TOF [us]",
  			"Peak Brilliance",
  			"Shape", t_0, t_1, nt,
  			&BRIL_shapeN[0],&BRIL_shape[0],&BRIL_shapeE[0],
  			ff);
    }
    BRIL_peakN[i] = BRIL_N[jmax][i];
    BRIL_peak[i] = BRIL_p[jmax][i];
    BRIL_peakE[i] = BRIL_p2[jmax][i];
  }

  sprintf(ff, "Mean_%s",filename);
  DETECTOR_OUT_1D(
  	"Mean brilliance",
    "Wavelength [AA]",
    "Mean Brilliance",
    "Mean", lambda_0, lambda_1, nlam,
    &BRIL_meanN[0],&BRIL_mean[0],&BRIL_meanE[0],
    ff);
  sprintf(ff, "Peak_%s",filename);
  DETECTOR_OUT_1D(
  	"Peak brilliance",
    "Wavelength [AA]",
    "Peak Brilliance",
    "Peak", lambda_0, lambda_1, nlam,
    &BRIL_peakN[0],&BRIL_peak[0],&BRIL_peakE[0],
    ff);

  /* MPI related NOTE: Order is important here! The 2D-data used to generate wavelength-slices and calculate
     the peak brilliance should be done LAST, otherwise we will get a factor of MPI_node_count too much as
     scatter/gather has been performed on the arrays... */
  if (toflambda == 1) {
    sprintf(ff, "TOFL_%s",filename);
    DETECTOR_OUT_2D(
      "TOF-wavelength brilliance",
      "Time-of-flight [\\gms]", "Wavelength [AA]",
      t_0, t_1, lambda_0, lambda_1,
      nt, nlam,
      &BRIL_N[0][0],&BRIL_p[0][0],&BRIL_p2[0][0],
      filename);
  }
}
%}

FINALLY
%{
  destroy_darr2d(BRIL_N);
  destroy_darr2d(BRIL_p);
  destroy_darr2d(BRIL_p2);

  destroy_darr1d(BRIL_mean);
  destroy_darr1d(BRIL_meanN);
  destroy_darr1d(BRIL_meanE);
  destroy_darr1d(BRIL_peak);
  destroy_darr1d(BRIL_peakN);
  destroy_darr1d(BRIL_peakE);

  destroy_darr1d(BRIL_shape);
  destroy_darr1d(BRIL_shapeN);
  destroy_darr1d(BRIL_shapeE);
%}

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