/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Virtual_tripoli4_input.comp

* %I
* Written by: <a href="mailto:guillaume.campioni@cea.fr">Guillaume Campioni</a>
* Date: Sep 28th, 2001
* Version:  $Revision: 1.10 $
* Origin: <a href="http://www.serma.cea.fr/">SERMA</a>
* Release: McStas 1.8
* Modified by: E. Farhi: Added automatic orientation+beam info+intensity norm
*
* This component uses a file of recorded neutrons from the reactor monte carlo
* code TRIPOLI4.4 as a source of particles.
*
* %D
* This component generates neutron events from a file created using the
* TRIPOLI4 Monte Carlo code for nuclear reactors (as MCNP). It is used to
* calculate flux exiting from hot or cold neutron sources.
* Neutron position and velocity is set from the file. The neutron time is
* left at zero.
* Storage files from TRIPOLI4.4 contain several batches of particules, all
* of them having the same statictical weight.
*
* Note that axes orientation may be different between TRIPOLI4.4 and McStas.
* The component has the ability to center and orient the neutron beam to the Z-axis.
* It also may change the coordinate system from the Tripoli frame to the McStas one.
* The verbose mode is highly recommanded as it displays lots of useful informations,
* including the absolute intensity normalisation factor. All neutron fluxes in the
* instrument should be multiplied by this factor.
* The source total intensity is 1.054e18 for LLB/Saclay and 4.28e18 for ILL/Grenoble.
*
* Format of TRIPOLI4.4 event files is :
*
*   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 'source' from a Tripoli4 simulation event file for the ILL:
* COMPONENT source = Virtual_tripoli4_input(
*    file = "ILL_SFH.dat", intensity=4.28e18,
*    repeat_count  = 1, verbose = 1, autocenter="translate rotate")
*
* %P
* INPUT PARAMETERS
* file:         [str] Name of the Tripoli4.4 neutron input file
* repeat_count: [1]   Number of times the source must be generated.
*                           0 unactivates the component
* verbose:      [0|1] Displays additional informations if set to 1
* intensity:    [n/s] Initial total Source integrated intensity
* autocenter:   [str] String which may contain following keywords.
*                     "translate" or "center"  to center the beam area AT (0,0,0)
*                     "rotate"    or "orient"  to center the beam velocity along Z
*                     "change axes"            to change coordinate system definition
*                     Other words are ignored.
*
* OUTPUT PARAMETERS
* head:   [char]        header buffer
* nl:     [long]        nb of lines in header
* mean_x: [double]      \
* mean_y: [double]      source center coordinates
* mean_z: [double]      /
* angle2z:[double, rad] rotation angle required to orient beam along Z axis
* nbatch: [long]        number of read neutron batches
*
* %L
* <a href="http://www.nea.fr/html/dbprog/tripoli-abs.html">Tripoli</a>
* %L
* Virtual_tripoli4_output
* %L
* CAMPIONI Guillaume, Etude et Modelisation des Sources Froides de Neutrons, These de Doctorat, CEA Saclay/UJF (2004)
*
* %E
*
*******************************************************************************/

DEFINE COMPONENT Virtual_tripoli4_input
DEFINITION PARAMETERS (file=0, autocenter=0)
SETTING PARAMETERS (repeat_count=1, verbose=0, intensity=1)
OUTPUT PARAMETERS (hfile,head,nl,rep,begin_of_file,begin_batch,do_rotate,do_translate,do_axes,
                   bx,by,bz,angle2z,nbatch,nsize,mean_x, mean_y, mean_z)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx, sy, sz)

SHARE
%{
#ifndef TRIPOLI4_INPUT_DEFS
#define TRIPOLI4_INPUT_DEFS
/* define Tripoli File parsing functions */
#define T4_WORDSIZE   256
/* number of neutron events to read for file pre-analysis */
#define T4_ANALYSE    10000

#include <sys/stat.h>

/* functions to identify Tripoli keywords */
int begin_of_batch(const char * word )
{
  return strcmp(word,"BEGIN_OF_BATCH")==0?1:0;
}
int end_of_batch(const char * word)
{
  return strcmp(word,"END_OF_BATCH")==0?1:0;
}
int is_neutron(const char * word)
{
  return strcmp(word,"NEUTRON")==0?1:0;
}

/* tripoli_read_word:
   function that reads iteratively words in the file
   returns  EOF or last character
   word[0..l-1]
 */
int tripoli_read_word(FILE *hfile,char *word){
  char c;
  int  i=0;

  while((c=fgetc(hfile))!=EOF && isspace(c)){ /* skip spaces */ }

  if(c!=EOF){
    word[i++]=c;
    while((c=fgetc(hfile))!=EOF && i < T4_WORDSIZE) {
      if(isspace(c)) break;
      word[i++]=c;
    }
  }
  word[i]='\0';
  return c;
}

/* tripoli_create_neutron:
   function that reads values following the NEUTRON keyword
   and assigns neutron parameters.
   Warning: Tripoli Coord system (x y z) = (-x z y)
 */
int tripoli_create_neutron(FILE *hfile, double *x,double *y,double *z,
                                 double *vx, double *vy, double *vz,
                                 double *t,
                                 double *sx, double *sy, double *sz,
                                 double *p)
{
  double Mev2Joule=1.602e-13;
  double speed;
  int    exit_flag=0, ifield=0;
  char   word[T4_WORDSIZE];
  double field[8];

  while(ifield<8){
    if(tripoli_read_word(hfile,word)==EOF) { exit_flag=1; break; };
    field[ifield++]=strtod(word,NULL);
  }
  if(field[0] < 0 || exit_flag) { exit_flag=1; }
  else {
    speed=sqrt(2.* field[0]*Mev2Joule/MNEUTRON);
    *x=field[1]/100.0;
    *y=field[2]/100.0;
    *z=field[3]/100.0;
    *vx=field[4]*speed;
    *vy=field[5]*speed;
    *vz=field[6]*speed;
    *p=field[7];
    *sx=1.;*sy=*sz=0;
    *t=0;
  }
  return (exit_flag);
}
/* tripoli_get_header:
   function that gets header lines until the Tripoli start of a batch
 */
char **tripoli_get_header(FILE *hfile,long *nl){
  int i;
  char **head;
  char   word[T4_WORDSIZE];

  *nl=0;
  fgets(word,T4_WORDSIZE,hfile);
  while(strncmp(word,"BEGIN_OF_BATCH",14)!=0){
    fgets(word,T4_WORDSIZE,hfile);
    (*nl)++;
  }
  rewind(hfile);
  head=(char **)calloc(*nl,sizeof(char *));
  for(i=0;i<*nl;++i){
    fgets(word,T4_WORDSIZE,hfile);
    head[i]=(char *)malloc(T4_WORDSIZE);
    strcpy(head[i],word);
  }
  rewind(hfile);
  return head;
}

#endif

%}

DECLARE
%{
  int   rep=1;                    /* Neutron repeat count */
  FILE *hfile;                    /* Neutron input file handle */
  char **head;                    /*Tripol4 header*/
  long  nl;                       /* Number of lines in header*/
  long  first_batch=0;                /* position in file of first batch */
  long  begin_batch=0, end_batch=0;   /* current batch numbers */
  char  do_rotate=0, do_translate=0, do_axes=0;
  double mean_x =0, mean_y =0, mean_z =0;
  double bx,by,bz;
  double angle2z;
  long   nbatch=0;
  long   nsize =0;
%}

INITIALIZE
%{
  char  word[T4_WORDSIZE];
  char  exit_flag=0;  /* set to 1 if end of simulation */
  int   result_read;

  double mean_vx=0, mean_vy=0, mean_vz=0;
  double mean_dx=0, mean_dy=0, mean_dz=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;

  double n_neutrons=0;
  double n_neutrons_p=0;
  long   filesize  =0;

  struct stat stfile;

  /* Open neutron input file. */

  if (file) {
    stat(file,&stfile);
    filesize = stfile.st_size;
    hfile = fopen(file, "r");
  }

  if(!hfile)
  {
    fprintf(stderr, "Tripoli4_input: %s: Error: Cannot open input file.\n", NAME_CURRENT_COMP);
    exit(1);
  } else if (verbose)
    printf("Tripoli4_input: %s: opening Tripoli4 file '%s'\n", NAME_CURRENT_COMP, file);

  head = tripoli_get_header(hfile, &nl); /*  and reset to file start */

  while (tripoli_read_word(hfile,word) != EOF) {
    /* store position of first batch in file for repeat */
    first_batch=ftell(hfile);
    if (begin_of_batch(word)) {
      if (tripoli_read_word(hfile,word)!=EOF) {
        begin_batch=strtol(word,NULL,0);
        break;
      }
    }
    if (!strcmp(word,"SIZE") || !strcmp(word,"TAILLE")) {
      if (tripoli_read_word(hfile,word)!=EOF) nsize=strtol(word,NULL,0);
    }
  }

  if (verbose) printf("Analysing Tripoli4 file %s (%d events)\n", file, T4_ANALYSE);
  /* analyse Tripoli file: count neutrons, get beam center and mean speed
    * do that for the first 1e4 neutrons, and extrapolates to file size */
  while (!exit_flag && n_neutrons <= T4_ANALYSE) {
    result_read=tripoli_read_word(hfile,word);

    if (is_neutron(word)) {  /* if key word NEUTRON is found */
      double x,y,z,vx,vy,vz,t,p,sx,sy,sz;
      if (tripoli_create_neutron(hfile,&x,&y,&z,&vx,&vy,&vz,&t,&sx,&sy,&sz,&p)) {
        fprintf(stderr, "Tripoli4_input: %s: Error: Cannot get neutron %d in batch %d.\n",NAME_CURRENT_COMP,n_neutrons,begin_batch);
        exit_flag = 1;
      } else {
        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;
        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); }
        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;
      }
    }

    if (begin_of_batch(word)){/* if key word BEGIN_OF_BATCH is found */
      if (tripoli_read_word(hfile,word) != EOF) {
        begin_batch=strtol(word,NULL,0);
        nbatch++;
      }
    }

    if (result_read==EOF) {  /* normal end of file */
      exit_flag=1;
    }
  } /* end while */

  if (n_neutrons) {
    double n_count_extrapolated=0;
    double mean_v;
    long   end_analyse=0;
    double cx,cy,cz;
    end_analyse = ftell(hfile);
    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;
    /* now estimates total ncount */
    mean_v = sqrt(mean_vx*mean_vx+mean_vy*mean_vy+mean_vz*mean_vz);
    n_count_extrapolated = n_neutrons*(filesize-first_batch)/(end_analyse - first_batch);
    if (verbose) {
      double mean_k, mean_w=0, mean_L=0;

      mean_k = V2K*mean_v;
      if (mean_k) mean_L = 2*PI/mean_k;
      mean_w = VS2E*mean_v*mean_v;
      printf("Tripoli4 file %s\nContains about %g neutrons in %d batches\n", file, n_count_extrapolated, nbatch*(filesize-first_batch)/(end_analyse - first_batch));
      if (n_count_extrapolated > mcget_ncount())
        printf("    (will use only %.3g %% of file)\n", 100*mcget_ncount()/n_count_extrapolated);
      else
        printf("    (limiting simulation to %g neutrons)\n", n_count_extrapolated*repeat_count);
      printf("  Source 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("  Source 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);
    }
    if (autocenter) {
      if (strstr(autocenter, "rotate") || strstr(autocenter, "orient"))
        do_rotate    = 1;
      if (strstr(autocenter, "translate") || strstr(autocenter, "center"))
        do_translate = 1;
      if (strstr(autocenter, "change") &&
            (strstr(autocenter, "axis")   || strstr(autocenter, "axes")
          || strstr(autocenter, "system") || strstr(autocenter, "coord")) )
        do_axes      = 1;
    }
    /* compute the rotation matrix to make mean_v along the Z-axis */
    /* first normalize mean velocity (will be Z axis): c=v/|v| */
    cx = mean_vx/mean_v; cy = mean_vy/mean_v; cz = mean_vz/mean_v;

    /* compute angle to rotate in order to come back to Z axis */
    angle2z = acos(cz); /* in RAD */
    /* get rotation axis b=c x [0 0 1] */
    if (angle2z) {
      vec_prod(bx,by,bz, cx,cy,cz, 0,0,1);
    } else { /* already well oriented: nothing to do */
      do_rotate = 0;
    }
    if (verbose && (do_rotate || do_translate)) {
      printf("* Beam will be ");
      if (do_translate) printf("translated (in position) ");
      if (do_rotate)    printf("rotated (%.3g [deg] to Z-axis) ", angle2z*RAD2DEG);
      if (do_axes)      printf("axes-swaped Tripoli->McStas");
      printf("\n");
    }

  } else {
    fprintf(stderr, "Tripoli4_input: %s: Error: file %s neutrons does not contains any neutron. \n",NAME_CURRENT_COMP, file);
    exit(-1);
  }
  begin_batch = 0; nbatch=0;
  /* reposition at start of file (batch start) */
  if (fseek(hfile, first_batch,SEEK_SET)) {
    fprintf(stderr, "Tripoli4_input: %s: Error: Can not reset Tripoli4 file (fseek error at analyse). \n",NAME_CURRENT_COMP);
    exit_flag = 1;
  }

%}

TRACE
%{
  char exit_flag=0;  /* set to 1 if end of simulation */
  int  result_read;
  char word[T4_WORDSIZE];

  while(1) {

  result_read=tripoli_read_word(hfile,word);

  if (is_neutron(word)) {  /* if key word NEUTRON is found */
    if (tripoli_create_neutron(hfile,&x,&y,&z,&vx,&vy,&vz,&t,&sx,&sy,&sz,&p)) {
      fprintf(stderr, "Virtual_tripoli4_input: %s: Error: Cannot create neutron in batch %d.\n",NAME_CURRENT_COMP,begin_batch);
      fprintf(stderr,"  Finishing simulation\n");
      exit_flag = 1;
    } else {
      if (do_translate) { /* translate the beam to origin */
        x -= mean_x; y -= mean_y; z -= mean_z;
      }
      if (do_rotate) {    /* rotate the beam so that its main axis is along Z */
        double nvx, nvy, nvz;
        rotate(nvx,nvy,nvz, vx,vy,vz, angle2z, bx,by,bz);
        vx = nvx; vy=nvy; vz=nvz;
        rotate(nvx,nvy,nvz,  x, y, z, angle2z, bx,by,bz);
         x = nvx;  y=nvy;  z=nvz;
      }
      if (do_axes) {
        /* Tripoli/Vitess (xyz) - > McStas (zxy) */
        double tx, ty, tz;
        tx=x;  ty=y;  tz=z;   x=ty;  y=tz;  z=tx;
        tx=vx; ty=vy; tz=vz; vx=ty; vy=tz; vz=tx;
      }
    }
    SCATTER;
    break;
  }

  if (end_of_batch(word)) {  /* if key word END_OF_BATCH is found */
    if (tripoli_read_word(hfile,word) != EOF) {
    /* this is the number of finishing batch */
      end_batch=strtol(word,NULL,0);
      if (end_batch!=begin_batch && begin_batch) {
        /* finishing batch does not match the one we were reading ! */
        fprintf(stderr, "Tripoli4_input: %s: Warning: inconsistent batch numbers between END(%d) and BEGIN(%d).\n", NAME_CURRENT_COMP, end_batch,begin_batch);
      } else nbatch++;
    }
    else {
      fprintf(stderr, "Tripoli4_input: %s: Expect 'END_OF_BATCH', found '%s' (after batch %d).", word, begin_batch);
      exit_flag = 1;
    }
  }

  if (begin_of_batch(word)){/* if key word BEGIN_OF_BATCH is found */
    if (tripoli_read_word(hfile,word) != EOF) {
      begin_batch=strtol(word,NULL,0);
    }
    else {
      fprintf(stderr, "Tripoli4_input: %s: Expect 'BEGIN_OF_BATCH', found '%s' (after batch %d).", word, begin_batch);
      exit_flag = 1;
    }
  }

  if (result_read==EOF) {  /* normal end of file */
    rep++;
    if (rep < repeat_count) {
      /* reposition at start of file (batch start) */
      int ret = fseek(hfile, first_batch,SEEK_SET);
      if (ret) {
        fprintf(stderr, "Tripoli4_input: %s: Error: Can not repeat Tripoli4 file (fseek error at repeat %d). \n",NAME_CURRENT_COMP, rep);
        exit_flag = 1;
      }
      result_read = tripoli_read_word(hfile,word); /* should be batch start */
      begin_batch=strtol(word,NULL,0);
      if (verbose) printf("Tripoli4_input: %s: Start Batch Number %d (iteration %d)\n",NAME_CURRENT_COMP,begin_batch, rep);
    }
    else exit_flag=1;
  }

  if (exit_flag) { mcset_ncount(mcget_run_num()); ABSORB; }

} /* end while */


%}

FINALLY
%{
  if (head)  free(head);
  if (hfile) fclose(hfile);
  if (verbose) {
    printf("Tripoli4_input: %s: \n", NAME_CURRENT_COMP);
    printf("                %g neutrons generated\n", mcget_ncount());
    printf("                %d read batch (B) of intial %d neutrons (N)\n",nbatch, nsize);
    if (nbatch && nsize)
    printf("* Normalisation factor Intensity/N/B = %g [n/s]\n", intensity/nbatch/nsize);
  }
%}

END