DEFINE COMPONENT MyPSD_monitor_4PI
DEFINITION PARAMETERS (nx=90, ny=90, string filename, string tth_file, restore_neutron=0)
SETTING PARAMETERS (radius)
OUTPUT PARAMETERS (PSD_N, PSD_p, PSD_p2)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)

DECLARE
%{
  double PSD_N[nx][ny];
  double PSD_p[nx][ny];
  double PSD_p2[nx][ny];
  double TTH_N[nx];
  double TTH_p[nx];
  double TTH_p2[nx];
%}
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;
    }
    TTH_N[i] = 0;
    TTH_p[i] = 0;
    TTH_p2[i] = 0;
  } 
%}
TRACE
%{
  double t0, t1, phi, theta;
  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)+0.5));
    if(i >= nx)
      i = nx-1;                      /* Special case for phi = PI. */
    else if(i < 0)
      i = 0;
    theta=asin(y/radius);
    j = floor(ny*(theta+PI/2)/PI+0.5);
    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;

    theta = RAD2DEG * acos( vz / sqrt(vx*vx + vy*vy + vz*vz) );
    i = floor((theta)*nx/(180));

    TTH_N[i]++;
    TTH_p[i] += p;
    TTH_p2[i] += p*p;
    SCATTER;
  }
  if (restore_neutron) {
    RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
  }
%}

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);
    
  DETECTOR_OUT_1D(
    "Diffractogram",
    "2theta",
    "Intensity",
    "TTH", 0, 180, nx,
    &TTH_N[0],&TTH_p[0],&TTH_p2[0],
    tth_file);
    
%}

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

END