/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: PSD_monitor_rad
*
* %Identification
* Written by: Henrich Frielinghaus
* Date:       Sept 2004
* Origin:     FZ-Juelich/FRJ-2/IFF/KWS-2
* Version:    $Revision: 1.4 $
* Release:    McStas 1.9
* Modified:   Kim Lefmann
* Date:       7 June 2005
*
* Position-sensitive monitor with radially averaging.
*
* %D
* Radial monitor that allows for radial averaging.
* Comment: The intensity is given as two files:
*          1) a radial sum
*          2) a radial average (i.e. intensity per area)
*
* This monitor may be replaced by the Monitor_nD(options="auto radius") component
*
* Example: PSD_monitor_rad(rmax=0.2, nr=100, filename="Output.psd", filename_av="Output_av.psd")
*
* %P
* INPUT PARAMETERS:
*
* rmax:     Outer radius of detector (m)
* nr:       Number of concentric circles (1)
* filename: Name of file in which to store the detector image (text)
* filename_av: Name of file in which to store the averaged detector image (text)
*
* OUTPUT PARAMETERS:
*
* PSDr_N:      Array of neutron counts
* PSDr_p:      Array of neutron weight counts
* PSDr_p2:     Array of second moments
* PSDr_av_N:   Array of neutron counts, averaged
* PSDr_av_p:   Array of neutron weight counts, averaged
* PSDr_av_p2:  Array of second moments, averaged
*
* %E
*******************************************************************************/


DEFINE COMPONENT PSD_monitor_rad
DEFINITION PARAMETERS (nr=100, filename, filename_av)
SETTING PARAMETERS (rmax)
OUTPUT PARAMETERS (PSDr_N, PSDr_p, PSDr_p2, PSDr_av_p, PSDr_av_p2)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
  %{
    double PSDr_N[nr];
    double PSDr_p[nr];
    double PSDr_p2[nr];
    double PSDr_av_p[nr];
    double PSDr_av_p2[nr];
  %}
INITIALIZE
  %{
    int i;

    for (i=0; i<nr; i++)
     {
      PSDr_N[i]  = 0;
      PSDr_p[i]  = 0;
      PSDr_p2[i] = 0;
     }
  %}
TRACE
  %{
    int i;
    double radpos;

    PROP_Z0;

    radpos = sqrt(x*x+y*y);

    if (radpos < rmax)
    {
      i = floor(nr*radpos/rmax);
      PSDr_N[i]++;
      PSDr_p[i] += p;
      PSDr_p2[i]+= p*p;
      SCATTER;
    }
  %}
SAVE
  %{
    int i;
    for(i=0; i<nr; i++)
     {
      PSDr_av_p[i]  = PSDr_p[i]  / (PI*rmax*rmax/(nr*nr)*(2*i+1.0));
      PSDr_av_p2[i] = PSDr_p2[i]
          / (PI*rmax*rmax/(nr*nr)*(2*i+1.0)) / (PI*rmax*rmax/(nr*nr)*(2*i+1.0));
     }
    DETECTOR_OUT_1D(
        "PSD monitor radial sum",
        "Radius [m]",
        "Intensity",
        "r", 0.5*rmax/nr, rmax+0.5*rmax/nr, nr,
        &PSDr_N[0],&PSDr_p[0],&PSDr_p2[0],
        filename);
    DETECTOR_OUT_1D(
        "PSD monitor radial average",
        "Radius [m]",
        "Intensity/m^2",
        "r", 0.5*rmax/nr, rmax+0.5*rmax/nr, nr,
        &PSDr_N[0],&PSDr_av_p[0],&PSDr_av_p2[0],
        filename_av);
  %}

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

END