/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Beam_spy
*
* %I
*
* Written by: E. Farhi
* Date: Nov 2005
* Version: $Revision: 1.4 $
* Release: McStas 1.9
* Origin: Risoe
*
* Beam analyzer for previous component
*
* %D
* This component displays informations about the beam at the previous component
* position. No data file is produced.
* It behaves as the Monitor component, but No propagation to the Beam_spy is
* performed, and all events are analyzed.
* The component should be located at the same position as the previous one.
*
* %P
* Input parameters:
*
* (none)
*
* %L
* Monitor component
* %E
*******************************************************************************/

DEFINE COMPONENT Beam_spy
DEFINITION PARAMETERS ()
SETTING PARAMETERS ()
OUTPUT PARAMETERS (mean_vx, mean_vy, mean_vz, mean_dx, mean_dy, mean_dz, n_neutrons,
                   min_x, min_y, min_z, max_x, max_y, max_z, min_vx, min_vy, min_vz,
                   max_vx, max_vy, max_vz,  mean_x, mean_y, mean_z, mean_t)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
  double n_neutrons=0;
  double n_neutrons_p=0;
  double n_neutrons_p2=0;
  double mean_x =0, mean_y =0, mean_z =0;
  double mean_vx=0, mean_vy=0, mean_vz=0;
  double mean_dx=0, mean_dy=0, mean_dz=0;
  double mean_t =0, mean_v=0;
  double min_x = FLT_MAX, min_y = FLT_MAX, min_z = FLT_MAX;
  double max_x =-FLT_MAX, max_y =-FLT_MAX, max_z =-FLT_MAX;
  double min_vx= FLT_MAX, min_vy= FLT_MAX, min_vz= FLT_MAX;
  double max_vx=-FLT_MAX, max_vy=-FLT_MAX, max_vz=-FLT_MAX;
%}

TRACE
%{
  double v;

  mean_x  += p*x;  mean_y  += p*y;  mean_z  += p*z;
  mean_vx += p*vx; mean_vy += p*vy; mean_vz += p*vz;
  mean_t  += p*t;
  v = sqrt(vx*vx+vy*vy+vz*vz);
  if (v)
    { mean_dx += p*fabs(vx/v); mean_dy += p*fabs(vy/v); mean_dz += p*fabs(vz/v); mean_v += p*v; }
  if (x  < min_x)  min_x  = x;
  if (y  < min_y)  min_y  = y;
  if (z  < min_z)  min_z  = z;
  if (vx < min_vx) min_vx = vx;
  if (vy < min_vy) min_vy = vy;
  if (vz < min_vz) min_vz = z;
  if (x  > max_x)  max_x  = x;
  if (y  > max_y)  max_y  = y;
  if (z  > max_z)  max_z  = z;
  if (vx > max_vx) max_vx = vx;
  if (vy > max_vy) max_vy = vy;
  if (vz > max_vz) max_vz = z;
  n_neutrons++;
  n_neutrons_p  += p;
  n_neutrons_p2 += p*p;
%}

SAVE
%{
  double mean_k, mean_w=0, mean_L=0;
  Coords c;

  c = POS_A_CURRENT_COMP;
  /* display statitics */
  mean_x  /= n_neutrons_p;
  mean_y  /= n_neutrons_p;
  mean_z  /= n_neutrons_p;
  mean_vx /= n_neutrons_p;
  mean_vy /= n_neutrons_p;
  mean_vz /= n_neutrons_p;
  mean_dx /= n_neutrons_p;
  mean_dy /= n_neutrons_p;
  mean_dz /= n_neutrons_p;
  mean_t  /= n_neutrons_p;
  mean_v  /= n_neutrons_p;
  /* now estimates total ncount */

  mean_k = V2K*mean_v;
  if (mean_k) mean_L = 2*PI/mean_k;
  mean_w = VS2E*mean_v*mean_v;
  printf("Beam analysis for Component preceeding %s\n", NAME_CURRENT_COMP);
  printf("Absolute position AT (%g, %g, %g)\n", c.x, c.y, c.z);

  printf("  Beam size (full width in [m]):      ");
  printf("    dX=%g dY=%g dZ=%g\n", max_x-min_x, max_y-min_y, max_z-min_z);
  printf("  Beam center (in [m]):               ");
  printf("    X0=%g Y0=%g Z0=%g\n", mean_x, mean_y, mean_z);
  printf("  Beam divergence (full width in [deg]):");
  printf("  dVx=%g dVy=%g dVz=%g\n",
    atan(mean_dx)*RAD2DEG,
    atan(mean_dy)*RAD2DEG,
    atan(mean_dz)*RAD2DEG);
  printf("  Beam speed (in [m/s]):                ");
  printf("  Vx=%g Vy=%g Vz=%g\n", mean_vx, mean_vy, mean_vz);
  printf("  Beam mean energy:\n");
  printf("    speed=%g [m/s] energy=%g [meV]\n    wavelength=%g [Angs] wavevector=%g [Angs-1]\n", mean_v, mean_w, mean_L, mean_k);
  printf("  Mean arrival time: t=%g [s]\n", mean_t);
  DETECTOR_OUT_0D("Beam analyzer " NAME_CURRENT_COMP, n_neutrons, n_neutrons_p, n_neutrons_p2);
%}

MCDISPLAY
%{
  /* A bit ugly; hard-coded dimensions. */
  magnify("");
  line(0,0,0,0.2,0,0);
  line(0,0,0,0,0.2,0);
  line(0,0,0,0,0,0.2);
%}

END