/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: PSD_monitor_4PI
*
* %I
* Written by: Kim Lefmann and Kristian Nielsen
* Date: April 17, 1998
* Version: $Revision: 1.23 $
* Origin: Risoe
* Release: McStas 1.6
*
* Spherical position-sensitive detector.
*
* %D
* An (n times m) pixel spherical PSD monitor using a cylindrical projection.
* Mostly for test and debugging purposes.
*
* Example: PSD_monitor_4PI(radius=0.1,
            nx=90, ny=90, filename="Output.psd")
*
* %P
* INPUT PARAMETERS:
*
* radius:   Radius of detector (m)
* nx:       Number of pixel columns (1)
* ny:       Number of pixel rows (1)
* filename: Name of file in which to store the detector image (text)
*
* OUTPUT PARAMETERS:
*
* PSD_N:    Array of neutron counts
* PSD_p:    Array of neutron weight counts
* PSD_p2:   Array of second moments
*
* %L
* <A HREF="http://neutron.risoe.dk/mcstas/components/tests/powder/">Test
* results</A> (not up-to-date).
*
* %E
*******************************************************************************/


DEFINE COMPONENT PSD_monitor_4PI
DEFINITION PARAMETERS (nx=90, ny=90, filename)
SETTING PARAMETERS (radius)
OUTPUT PARAMETERS (PSD_N, PSD_p, PSD_p2)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
  double PSD_N[nx][ny];
  double PSD_p[nx][ny];
  double PSD_p2[nx][ny];
%}
INITIALIZE
%{
  int i,j;

  for (i=0; i<nx; i++)
    for (j=0; j<ny; j++)
    {
      PSD_N[i][j] = 0;
      PSD_p[i][j] = 0;
      PSD_p2[i][j] = 0;
    }
%}
TRACE
%{
  double t0, t1, phi;
  int i,j;

  if(sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius) && t1 > 0)
  {
    if(t0 < 0)
      t0 = t1;
    /* t0 is now time of intersection with the sphere. */
    PROP_DT(t0);
    phi = atan2(x,z);
    i = floor(nx*(phi/(2*PI)));
    if(i >= nx)
      i = nx-1;                      /* Special case for phi = PI. */
    else if(i < 0)
      i = 0;
    j = floor(ny*(y/(2*radius)));
    if(j >= ny)
      j = ny-1;                      /* Special case for y = radius. */
    else if(j < 0)
      j = 0;
    PSD_N[i][j]++;
    PSD_p[i][j] += p;
    PSD_p2[i][j] += p*p;
    SCATTER;
  }
%}

SAVE
%{
  DETECTOR_OUT_2D(
    "4PI PSD monitor",
    "Longitude [deg]",
    "Lattitude [deg]",
    -180, 180, -90, 90,
    nx, ny,
    &PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0],
    filename);
%}

MCDISPLAY
%{
  magnify("");
  circle("xy",0,0,0,radius);
  circle("xz",0,0,0,radius);
  circle("yz",0,0,0,radius);
%}

END