/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Virtual_tripoli4_output
*
* %I
* Written by: <a href="mailtoguillaume.campioni@cea.fr">Guillaume Campioni</a>
* Date: Sep 28th, 2001
* Version:  $Revision: 1.3 $
* Origin: <a href="http://www.cea-llb.fr">LLB</a>
* Release: McStas 1.8
*
* Detector-like component that writes neutron state parameters into a
* 'virtual source' neutron file when neutrons come from the source :
* Virtual_tripoli4_input.comp
*
* %D
* Detector-like component writing neutron state parameters to a
* virtual source neutron file when neutron are coming from a
* Virtual_tripoli4_input.comp.
* As McStas keeps the notion of batch with this component, it's possible
* to calculate standard deviations for the flux in differents energy groupes.
* The component geometry is the full plane, and saves the neutron state as
* it exits from the previous component.
* Format is the one used by TRIPOLI4.4 stock files :
*
*   NEUTRON energy position_X position_Y position_Z dir_X dir_Y dir_Z weight
*
* energy is in Mega eV
* positions are in cm and the direction vector is normalized to 1.
*
*
* EXAMPLE:
* To create a file collecting all neutron states with TRIPOLI4 format
* COMPONENT fichier_sortie = Virtual_tripoli4_output(
*    file = "exit_guide_result.dat",
*    virtual_tripoli4_input_comp = source_file_tripoli )
*  at the position where will be the Virtual_tripoli4_input.
*
* %P
* INPUT PARAMETERS
* file:         [str] name of the Tripoli4 neutron output file,
*                     or stdout if left to 0.
* Virtual_tripoli_input_comp   [comp] Name of the Tripoli4_input component
*                              from which is to be obtained batch information
*                              use '0' if not defined.
* batch         [1]   index of the Tripoli batch to generate, when no
*                     Tripoli4_input component is available in instrument.
*
* %L
* <a href="http://www.nea.fr/html/dbprog/tripoli-abs.html">Tripoli</a>
* Virtual_tripoli4_input
*
* %E
*******************************************************************************/


DEFINE COMPONENT Virtual_tripoli4_output
DEFINITION PARAMETERS (file=0, virtual_tripoli4_input_comp=0)
SETTING PARAMETERS (batch=1)
OUTPUT PARAMETERS (hfile,previous_batch)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)


DECLARE
%{
long  previous_batch=0;
FILE *hfile;
%}

INITIALIZE
%{
  int i;
  char ** head;
  long nl;
  /* Open neutron output file. */
  if(file)
    hfile = fopen(file, "w");
  else
    hfile = stdout;
  if(!hfile)
  {
    fprintf(stderr, "Tripoli4_output: %s: Error: Cannot open output file %s.\n", NAME_CURRENT_COMP, (file ? file : "stdout"));
    exit(1);
  }
  fprintf(hfile,"# Tripoli4  batch file %s written by MCSTAS (%s)\n", (file ? file : "stdout"), MCSTAS_VERSION);
  fprintf(hfile,"# Generated by component %s in instrument %s (%s)\n", NAME_CURRENT_COMP, mcinstrument_name, mcinstrument_source);
  fprintf(hfile,"# Format: Tripoli4 data text file with header\n");
#ifdef TRIPOLI4_INPUT_DEFS
  /* use header information from Virtual_tripoli4_input component instance */
  head = (MC_GETPAR(virtual_tripoli4_input_comp, head));/* header buffer */
  nl   = (MC_GETPAR(virtual_tripoli4_input_comp, nl));  /* nb of lines in header */
  for(i=0;i<nl;++i){
    fprintf(hfile,head[i]);
  }
#else
  if (batch <= 0)
  {
    fprintf(stderr, "Tripoli4_output: %s: Error: Negative batch number (%d). Using 1.\n", NAME_CURRENT_COMP, (long)batch);
    batch = 1;
  }
  /* generates a default Tripoli 4 header */
  fprintf(hfile,"# Data on the original simulation which generates this file : \n";
  fprintf(hfile,"NB_BATCH   1\n");
  fprintf(hfile,"TAILLE     Unknown\n");
  fprintf(hfile,"PARTICULE  NEUTRON\n");
  fprintf(hfile,"INTENSITE  %g\n",get_ncount());
  fprintf(hfile,"SIMULATION PROTECTION\n");
  fprintf(hfile,"Following lines are: type, energy in MeV, position = (x,y,z) in cm, direction = (u,v,w), statistical weight\n");
#endif

%}

TRACE
%{
  double Mev2Joule = 1.602e-13;
  double Mneutron  = 1.6749543e-27;
  long   current_batch=0;

  double energy, speed, weight;
  double ux,uy,uz,pos_x,pos_y,pos_z;

#ifdef TRIPOLI4_INPUT_DEFS
  current_batch = (MC_GETPAR(virtual_tripoli4_input_comp, num_batch));
#else
  current_batch = batch;
#endif
  if (current_batch != previous_batch)
  {
    if (previous_batch == 0) {
      /* write 1st batch header */
      fprintf(hfile,"------------------------------------------------------------------------------------------------------------\n");
      fprintf(hfile, "BEGIN_OF_BATCH  %d\n\n", current_batch);
    } else {
      fprintf(hfile, "\nEND_OF_BATCH  %d\n",   previous_batch);
      fprintf(hfile,"------------------------------------------------------------------------------------------------------------\n");
      fprintf(hfile,"------------------------------------------------------------------------------------------------------------\n");
      fprintf(hfile, "BEGIN_OF_BATCH  %d\n\n", current_batch);
    }
    previous_batch = current_batch;
  }
  /* writes Tripoli file line by line */
  speed = vx*vx + vy*vy + vz*vz;
  energy = speed / 2. * Mneutron / Mev2Joule;

  speed = sqrt(speed);
  ux = vx/speed;
  uy = vy/speed;
  uz = vz/speed;

  pos_x = x * 100.;
  pos_y = y * 100.;
  pos_z = z * 100.;

  weight = p;
//printf("NEUTRON FINAL MC\t%f  %f  %f  %f  %f  %f  %f  %f  %f  %f  %e \n",x,y,z,vx,vy,vz,t,sx,sy,sz,p);
//printf("NEUTRON FINAL TR\t%e \t%f \t%f \t%f \t%f \t%f \t%f \t%e \n", energy, pos_x, pos_y, pos_z, ux, uy, uz, weight);

  fprintf(hfile, "NEUTRON \t%e \t%f \t%f \t%f \t%f \t%f \t%f \t%e \n", energy, pos_x, pos_y, pos_z, ux, uy, uz, weight);

%}

SAVE
%{
  if (hfile) {
    fprintf(hfile, "\nEND_OF_BATCH  %d\n",   previous_batch);
    fprintf(hfile,"------------------------------------------------------------------------------------------------------------\n");
    fflush(hfile);
  }

%}

FINALLY
%{
  if (hfile) fclose(hfile);
%}

END