/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2006, 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
* Origin: <a href="http://www.serma.cea.fr/">SERMA</a>
* Modified by: E. Farhi: Added automatic orientation+beam info+intensity norm
*
* This component reads 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 changes the coordinate system from the Tripoli frame to the McStas one.
* The verbose mode is highly recommended 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. Such a renormalization is done
* when 'autocenter' contains the word 'rescale'.
* The source total intensity is 1.054e18 for LLB/Saclay (14 MW) and 4.28e18 for
* ILL/Grenoble (58 MW).
*
* 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.
*
* %BUGS
* We recommend NOT to use parallel execution (MPI) with this component. If you
* need to, set parameter 'smooth=1'.
*
* EXAMPLE:
* To create a 'source' from a Tripoli4 simulation event file for the ILL:
* COMPONENT source = Virtual_tripoli4_input(
*    filename = "ILL_SFH.dat", intensity=4.28e18,
*    verbose = 1, autocenter="translate rotate rescale")
*
* %P
* INPUT PARAMETERS
* filename: [str]         Name of the Tripoli4.4 neutron input file. Empty string "" unactivates component
* 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 "rescale" to adapt intensity to abs. units. Other words are ignored.
* T4_ANALYSE: [1]         Number of neutron events to read for file pre-analysis. Use 0 to analyze them all.
* T4_ANALYSE_EMIN: [meV]  Minimal energy to use for file pre-analysis 
* T4_ANALYSE_EMAX: [meV]  Maximal energy to use for file pre-analysis 
* radius: [m]             In the case the Tripoli4 batch file is a point, you may specify a disk emission area for the source.
* smooth: [0/1]           Smooth sparsed event files for file repetitions.
*
* OUTPUT PARAMETERS
* head: [char]            header buffer
* nl: [long]              nb of lines in header
* mean_x: [m]             source center coordinate
* mean_y: [m]             source center coordinate
* mean_z: [m]             source center coordinate
* angle2z: [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

SETTING PARAMETERS (string filename=0, string autocenter=0, repeat_count=1, verbose=0, intensity=1,
T4_ANALYSE=10000, radius=0, T4_ANALYSE_EMIN=0, T4_ANALYSE_EMAX=0,smooth=1)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
#ifndef TRIPOLI4_INPUT_DEFS
#define TRIPOLI4_INPUT_DEFS
/* define Tripoli File parsing functions */
#define T4_WORDSIZE   256

#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;
    // now change axis system
    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;
  }
  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, *s;
  char   word[T4_WORDSIZE];

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

#endif

%}

DECLARE
%{
  int rep;                       /* Neutron repeat count */
  int repeat_cnt;                /* Repeat count, MPI taken into account */
  FILE *hfile;                   /* Neutron input file handle */
  char **head;                   /* Tripol4 header*/
  long nl;                       /* Number of lines in header*/
  long first_batch;              /* position in file of first batch */
  long begin_batch, end_batch;   /* current batch numbers */
  char do_rotate, do_translate;
  double mean_x, mean_y, mean_z;
  double bx,by,bz;
  double angle2z;
  long nbatch;
  long nsize;
  double intensity_extrapolated;

  double min_x, min_y, min_z;
  double max_x, max_y, max_z;
  double min_vx, min_vy, min_vz;
  double max_vx, max_vy, max_vz;
  double n_count_extrapolated;
%}

INITIALIZE
%{
  rep=1;
  first_batch=0;
  begin_batch=0;
  end_batch=0;
  do_rotate=0;
  do_translate=0;
  mean_x=0;
  mean_y=0;
  mean_z=0;
  angle2z=0;
  nbatch=0;
  nsize=0;
  intensity_extrapolated=0;
  min_x=FLT_MAX;
  min_y=FLT_MAX;
  min_z=FLT_MAX;
  max_x=-FLT_MAX;
  max_y=-FLT_MAX;
  max_z=-FLT_MAX;
  min_vx=FLT_MAX;
  min_vy=FLT_MAX;
  min_vz=FLT_MAX;
  max_vx=-FLT_MAX;
  max_vy=-FLT_MAX;
  max_vz=-FLT_MAX;
  n_count_extrapolated=0;


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

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

  struct stat stfile;

  /* Open neutron input file. */
  /* (If empty file given, present warning in case of verbose, perform nothing) */

  if (!strcmp(filename,"")) {
    if (verbose) printf("Tripoli4_input: %s: Empty file given, doing nothing!\n", NAME_CURRENT_COMP);
  } else {
    printf("non-empty file given\n");
    if (filename) {
      stat(filename,&stfile);
      filesize = stfile.st_size;
      hfile = fopen(filename, "r");
    }

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

    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 (", filename);
      if (T4_ANALYSE) printf("%g events)\n", T4_ANALYSE);
      else printf("all events)\n");
    }

    /* 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 || !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 %g in batch %ld.\n",NAME_CURRENT_COMP,n_neutrons,begin_batch);
	  exit_flag = 1;
	} else {
	  double v,w;
	  /* if an energy filter is set, use it */
	  v = sqrt(vx*vx+vy*vy+vz*vz);
	  w = VS2E*v*v;
	  if ( (!T4_ANALYSE_EMIN && !T4_ANALYSE_EMAX)
	       || (T4_ANALYSE_EMIN <= w && w <= T4_ANALYSE_EMAX) ) {
	    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_v += p*v;
	    angle2z += p*vz/v;

	    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_p += p;
	  }
	  n_neutrons++;
	}
      }

      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) {
      long   end_analyse=0;
      double cx,cy,cz;
      double nbatch_extrapolated=0;
      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;
      mean_v  /= n_neutrons_p;
      angle2z /= n_neutrons_p;
      /* now estimates total ncount */
      if (T4_ANALYSE>0) {
	n_count_extrapolated = n_neutrons*fabs(filesize-first_batch)/(end_analyse - first_batch);
	nbatch_extrapolated  = nbatch    *fabs(filesize-first_batch)/(end_analyse - first_batch);
      } else {
	n_count_extrapolated = n_neutrons;
	nbatch_extrapolated  = nbatch;
      }
      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 %s%g neutrons in %ld batches of size %ld\n",
	       filename, T4_ANALYSE ? "about " : "",
	       n_count_extrapolated, (long)nbatch_extrapolated, nsize);
	if (n_count_extrapolated > mcget_ncount())
	  printf("    (will use only %.3g %% of file)\n", 100.0*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);
	if (fabs((max_x-min_x)*(max_y-min_y)) < 1e-6 && radius <= 0)
	  printf("WARNING: source area is a point. Use radius > 0 to make it a disk\n");
	printf("  Source center (in [m]):               ");
	printf("    X0=%g Y0=%g Z0=%g\n", mean_x, mean_y, mean_z);
	printf("  Beam velocity divergence (half 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, "rescale") && intensity && nbatch_extrapolated && nsize) {
	  intensity_extrapolated = fabs(intensity/nbatch_extrapolated/nsize);
	  printf("* Automatic Normalisation factor Intensity/N/B = %g [n/s]\n",
		 intensity_extrapolated);
	}
      }
      /* 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 */
      angle2z = acos(angle2z);
      /* 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 || smooth)) {
	printf("* Beam will be ");
	if (do_translate) printf("translated (in position) ");
	if (do_rotate)    printf("rotated (%.3g [deg] to Z-axis) ", angle2z*RAD2DEG);
	if (smooth)       printf("smoothed ");
	printf("\n");
      }

    } else {
      fprintf(stderr, "Tripoli4_input: %s: Error: file %s neutrons does not contains any neutron. \n",NAME_CURRENT_COMP, filename);
      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;
    }

#if defined (USE_MPI)
    if (!smooth && mpi_node_count > 1) {
      if (verbose)
	printf("Tripoli4_input: %s: smoothing (smooth=1) is recommended when running MPI execution\n", NAME_CURRENT_COMP);
    }
#endif

    double min_dv=fabs(max_vx-min_vx);
    if (min_dv > fabs(max_vy-min_vy)) min_dv = fabs(max_vy-min_vy);
    if (min_dv > fabs(max_vz-min_vz)) min_dv = fabs(max_vz-min_vz);
    min_vx = min_dv;
  } /* non-empty file */

  repeat_cnt = repeat_count;
  #if defined (USE_MPI)
    repeat_cnt = ceil(1.0*repeat_cnt/mpi_node_count);
  #endif
%}

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 %ld.\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 (radius > 0) /* force a disk source area */
      {
        double chi,r;
        chi=2*PI*rand01();                          /* Choose point on source */
        r=sqrt(rand01())*radius;                    /* with uniform distribution. */
        x=r*cos(chi);
        y=r*sin(chi);
      }
      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 (intensity_extrapolated) p *= fabs(intensity_extrapolated/repeat_cnt);
    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(%ld) and BEGIN(%ld).\n", NAME_CURRENT_COMP, end_batch,begin_batch);
      } else nbatch++;
    }
    else {
      fprintf(stderr, "Tripoli4_input: %s: Expect 'END_OF_BATCH', found '%s' (after batch %ld).", NAME_CURRENT_COMP, 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 %ld).", NAME_CURRENT_COMP, word, begin_batch);
      exit_flag = 1;
    }
  }

#if defined (USE_MPI)
    if (smooth && n_count_extrapolated)
#else
    if (smooth && rep > 1 && n_count_extrapolated)
#endif
    {
      /* apply smmothing */
      x += randnorm()*(max_x-min_x)/n_count_extrapolated/2;
      y += randnorm()*(max_y-min_y)/n_count_extrapolated/2;
      z += randnorm()*(max_z-min_z)/n_count_extrapolated/2;
      vx += randnorm()*min_vx/n_count_extrapolated/2;
      vy += randnorm()*min_vx/n_count_extrapolated/2;
      vx += randnorm()*min_vx/n_count_extrapolated/2;
    }

  if (result_read==EOF) {  /* normal end of file */
    rep++;
    if (rep <= repeat_cnt) {
      /* 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 %ld (iteration %d)\n",NAME_CURRENT_COMP,begin_batch, rep);
    }
    else exit_flag=1;
  }
  #ifdef USE_MPI
      /* We always repeat by the number of nodes in an MPI run */
      p /= mpi_node_count;
  #endif
  if (exit_flag) { mcset_ncount(mcget_run_num()); ABSORB; }

} /* end while */


%}

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

MCDISPLAY
%{
  
  double width =(max_x-min_x)/2;
  double height=(max_y-min_y)/2;
  multiline(5, -width/2.0, -height/2.0, 0.0,
                width/2.0, -height/2.0, 0.0,
                width/2.0,  height/2.0, 0.0,
               -width/2.0,  height/2.0, 0.0,
               -width/2.0, -height/2.0, 0.0);
%}

END