From farhi at ill.fr  Wed Jul 11 09:54:25 2001
From: farhi at ill.fr (Emmanuel Farhi)
Date: Wed, 11 Jul 2001 09:54:25 +0200
Subject: New instrument simulations available at McStas@ill web pages
Message-ID: <3B4C0631.CA768577@ill.fr>

An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20010711/e00ed10a/attachment.html>

From kim.lefmann at risoe.dk  Wed Jul 25 13:25:08 2001
From: kim.lefmann at risoe.dk (Kim Lefmann)
Date: Wed, 25 Jul 2001 13:25:08 +0200 (CEST)
Subject: Reflectivity of PG monochromator crystals
Message-ID: <Pine.LNX.3.96.1010725132156.12637B-100000@fys-lin-1.risoe.dk>

Dear simulators,

Do any of you happen to know any published (or unpublished) work on the
reflectivity of PG crystals for monochromators and analysers? I am
currently developing new and more realistic basic components - and I need
this information for the new monochromator component. 

Hope to see many of you in Munich in 6 weeks. 

Yours,
                  Kim

------------------------------
Kim Lefmann
Senior Scientist
Dept. Cond. Matt. Phys. & Chem.
Risoe National Laboratory

Phone: +45 46 77 47 26
Fax:   +45 46 77 47 90



From farhi at ill.fr  Thu Jul 26 13:57:46 2001
From: farhi at ill.fr (Emmanuel Farhi)
Date: Thu, 26 Jul 2001 13:57:46 +0200
Subject: Monitor_nD v0.15
Message-ID: <3B6005BA.1F60CB96@ill.fr>

An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20010726/40faa673/attachment.html>
-------------- next part --------------
/*******************************************************************************
*
* Component: Monitor_nD
*
* %Identification
* Written by: <a href="mailto:farhi at ill.fr">Emmanuel Farhi</a>
* Date: 14th Feb 2000.
* Origin: <a href="http://www.ill.fr">ILL (France)</a>
* Version: 0.15
* Modified by: EF, 29th Feb 2000 : added more options, monitor shape, theta, phi
* Modified by: EF, 10th Mar 2000 : use struct. 
* Modified by: EF, 25th May 2000 : correct Vars.Mon2D_Buffer alloc bug.
* Modified by: EF, 17th Oct 2000 : spin detector ok. (0.13.4)
* Modified by: EF, 19th Jan 2001 : 'auto limits' bug corrected (0.13.5)
* Modified by: EF, 01st Feb 2001 : PreMonitor for correlation studies (0.13.6)
* Modified by: EF, 02nd Feb 2001 : user variable (0.13.7)
* Modified by: EF, 2th  Feb 2001 : 3He gas absorption handling (0.13.8)
* Modified by: EF, 5th  Mar 2001 : intermediate savings (0.13.9)
* Modified by: EF, 5th  Apr 2001 : use global functions (0.14) compile faster
* Modified by: EF, 18th Apr 2001 : passes DETECTOR_OUT to mcdetector_out 0.14.1
* Modified by: EF, 23th Jul 2001 : log of signal, init arrays to 0, box (0.15)
*
* This component is a general Monitor that can output 0/1/2D signals 
* (Intensity vs. [something] and vs. [something] ...)
*
* %Description
* This component is a general Monitor that can output 0/1/2D signals 
* It can produce many 1D signals (one for any variable specified in 
* option list), or a single 2D output (two variables related).
* Also, an additional 'list' of neutrons can be produced.
* By default, monitor is square (in x/y plane). disk is also possible
* The 'cylinder' option will change that for banana shaped detector
* The 'sphere' option simulates spherical detector. The 'box' is a box.
* In normal configuration, the Monitor_nD measures the current parameters 
* of the neutron that is beeing detected. But a PreMonitor_nD component can 
* be used in order to study correlations between a neutron being detected in
* a Monitor_nD place, and given parameters that are monitored elsewhere 
* (at <b>PreMonitor_nD</b>).
* It is also possible, using the 'intermediate' keyword, to save monitor results
* every X percent of the simulation. The monitor can also act as a 3He gas
* detector, taking into account the detection efficiency.
* 
* <b>Possible options are</b>
* Variables to record:
*     kx ky kz k wavevector (Angs-1) [    usually axis are
*     vx vy vz v            (m/s)         x=horz., y=vert., z=on axis]
*     x y z radius          (m)      Distance, Position
*     t time                (s)      Time of Flight
*     energy omega          (meV)
*     lambda wavelength     (Angs)
*     p intensity flux      (n/s) or (n/cm^2/s)
*     ncounts               (1)
*     sx sy sz              (1)      Spin
*     vdiv ydiv dy          (deg)    vertical divergence (y)
*     hdiv divergence xdiv  (deg)    horizontal divergence (x)
*     angle                 (deg)    divergence from <z> direction
*     theta longitude       (deg)    longitude (x/z) [for sphere and cylinder]
*     phi   lattitude       (deg)    lattitude (y/z) [for sphere and cylinder]
*
*     user user1            will monitor the [Mon_Name]_Vars.UserVariable{1|2}
*     user2                 to be assigned in an other component (see below)
*
* <b>Other options are:</b>
*     auto {limits}             Automatically set detector limits
*     all  {limits or bins}     To set all limits or bins values
*     bins=[bins=20]            Number of bins in the detector along dimension
*     borders                   To also count off-limits neutrons
*                                (X < min or X > max)
*     file=string               Detector image file name. default is component
*                                name, plus date and variable extension.
*     limits=[min max]          Lower/Upper limits for axes
*                                (see up for the variable unit)
*     list=[counts=1000] or all For a long file of neutron characteristics
*                                with [counts] or all events
*     multiple                  For 1D monitors into multiple independant files
*     no or not                 Revert next option
*     outgoing                  Monitor outgoing beam in non flat (sph/cyl) det
*                                (default is incoming beam)
*     per cm2                   Intensity will be per cm^2
*     slit or absorb            Absorb neutrons that are out detector
*     unactivate                To unactivate detector (0D detector)
*     premonitor                Will monitor neutron parameters stored
*                                previously with <b>PreMonitor_nD</b>.
*     verbose                   To display additional informations
*     3He_pressure=[3 in bars]  The 3He gas pressure in detector.
*                               3He_pressure=0 means perfect detector (default)
*     intermediate=[5 in %]     Save results every n% steps in simulation
*     log                       Will monitor the log10 of the signal
*     parallel                  Use this option when the next component is at 
*                                the same position (parallel components)
* Detector shape options (specified as xwidth, yheigth, zthick or x/y/z/min/max)
*     box                       Box of size xwidth, yheigth, zthick
*     cylinder                  To get a cylindrical monitor (e.g. banana)
*                                (radius is xwidth, height is yheigth).
*     disk                      Disk flat xy monitor. radius is xwidth.
*     sphrere                   To get a spherical monitor (e.g. a 4PI)
*                                (radius is xwidth).
*     square                    Square flat xy monitor (xwidth, yheigth)
*
* <b>EXAMPLES:</b>
* MyMon = Monitor_nD(
*   xwidth = 0.1, yheigth = 0.1, zthick = 0,
*   options = "intensity per cm2 angle,limits=[-5 5] bins=10,with
*              borders, file = mon1");
*                           will monitor neutron angle from [z] axis, between -5
*                        and 5 degrees, in 10 bins, into "mon1.A" output 1D file
*   options = "sphere theta phi outgoing"   for a sphere PSD detector (out beam)
*                                    and saves into file "MyMon_[Date_ID].th_ph"
*   options = "angle radius auto limits"   is a 2D monitor with automatic limits
*   options = "list kx ky kz energy" records each neutron contribution in 1 file
*   options = "multiple kx ky kz energy and list all neutrons" 
*          makes 4 output 1D files and produces a complete list for all neutrons
*
*
* How to monitor anything in my simulation: with the 'user' option
*  In a component, you will put for instance in INITIALIZE and/or
*  TRACE sections (same is valid with user2)
*
*  struct MonitornD_Variables *Vars = &(MC_GETPAR(Guide_Mon, Vars));
*               with name of monitor that you will use in instr 
*  strcpy(Vars->UserName1,"My variable label");
*  Vars->UserVariable1++;
*
*  and in the instrument description:
*
*  COMPONENT Guide_Mon = Monitor_nD(
*   xwidth = 0.1, yheigth = 0.1, zthick = 0,
*   options="user1, limits=[0 15], bins=15")
*  AT (0,0,0) RELATIVE GuideEnd
*
* See also the example in <a href="PreMonitor_nD.html">PreMonitor_nD</a>.
*
* %Parameters
* INPUT PARAMETERS:
*
* xwidth:  (m)  Width/radius of detector (x)
* yheigth: (m)  Height of detector (y)
* zthick:  (m)  Thichness of detector (z)
* options: (str)  String that specifies the configuration of the monitor</br>
*                 The general syntax is "[x] options..." (see <b>Descr.</b>)
*
* Optional input parameters (override xwidth yheigth zthick):
* xmin:   (m)    Lower x bound of opening
* xmax:   (m)    Upper x bound of opening
* ymin:   (m)    Lower y bound of opening
* ymax:   (m)    Upper y bound of opening
* zmin:   (m)    Lower z bound of opening
* zmax:   (m)    Upper z bound of opening
*
* OUTPUT PARAMETERS:
*
* DEFS: structure containing Monitor_nD Defines (struct)
* Vars: structure containing Monitor_nD variables (struct)
*
* %Link 
* <a href="http://www.ill.fr/tas/mcstas/">McStas at ILL</a>
* <a href="PreMonitor_nD.html">PreMonitor_nD</a>
*
* %End
*******************************************************************************/

DEFINE COMPONENT Monitor_nD
DEFINITION PARAMETERS (options)
SETTING PARAMETERS (xwidth=.1, yheigth=.1, zthick=0, xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0)
/* these are protected C variables */
OUTPUT PARAMETERS (DEFS, Vars)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)

DECLARE
%{

#ifndef FLT_MAX
#define FLT_MAX 1E37
#endif
#ifndef Monitor_nD_Version

#define Monitor_nD_Version "0.15"
#define MONnD_COORD_NMAX  30  /* max number of variables to record */
  
  /* here we define some DEFINE constants */
  typedef struct MonitornD_Defines
  {
    char COORD_NONE  ;
    char COORD_X     ;
    char COORD_Y     ;
    char COORD_Z     ;
    char COORD_VX    ;
    char COORD_VY    ;
    char COORD_VZ    ;
    char COORD_T     ;
    char COORD_P     ;
    char COORD_SX    ;
    char COORD_SY    ;
    char COORD_SZ    ;
    char COORD_KX    ;
    char COORD_KY    ;
    char COORD_KZ    ;
    char COORD_K     ;
    char COORD_V     ; 
    char COORD_ENERGY; 
    char COORD_LAMBDA; 
    char COORD_RADIUS; 
    char COORD_HDIV  ; 
    char COORD_VDIV  ; 
    char COORD_ANGLE ; 
    char COORD_NCOUNT; 
    char COORD_THETA ; 
    char COORD_PHI   ; 
    char COORD_USER1 ; 
    char COORD_USER2 ;

    /* token modifiers */
    char COORD_VAR   ; /* next token should be a variable or normal option */
    char COORD_MIN   ; /* next token is a min value */
    char COORD_MAX   ; /* next token is a max value */
    char COORD_DIM   ; /* next token is a bin value */
    char COORD_FIL   ; /* next token is a filename */
    char COORD_EVNT  ; /* next token is a buffer size value */
    char COORD_3HE   ; /* next token is a 3He pressure value */
    char COORD_INTERM; /* next token is an intermediate save value (percent) */

    char TOKEN_DEL[32]; /* token separators */

    char SHAPE_SQUARE; /* shape of the monitor */
    char SHAPE_DISK  ; 
    char SHAPE_SPHERE; 
    char SHAPE_CYLIND; 
    char SHAPE_BOX   ; 
    
  } MonitornD_Defines_type;
  
  typedef struct MonitornD_Variables
  {
    double area;
    double Sphere_Radius     ;
    double Cylinder_Height   ;
    char   Flag_With_Borders ;   /* 2 means xy borders too */
    char   Flag_List         ;   /* 1 store 1 buffer, 2 is list all */
    char   Flag_Multiple     ;   /* 1 when n1D, 0 for 2D */
    char   Flag_Verbose      ;
    int    Flag_Shape        ;
    char   Flag_Auto_Limits  ;   /* get limits from first Buffer */
    char   Flag_Absorb       ;   /* monitor is also a slit */
    char   Flag_per_cm2      ;   /* flux is per cm2 */
    char   Flag_log          ;   /* log10 of the flux */
    char   Flag_parallel     ;   /* set neutron state back after detection (parallel components) */
    
    int    Coord_Number      ;   /* total number of variables to monitor, plus intensity (0) */
    long   Buffer_Block      ;   /* Buffer size for list or auto limits */
    long   Neutron_Counter   ;   /* event counter, simulation total counts is mcget_ncount() */
    long   Buffer_Counter    ;   /* index in Buffer size (for realloc) */
    long   Buffer_Size       ;
    char   Coord_Type[MONnD_COORD_NMAX];    /* type of variable */
    char   Coord_Label[MONnD_COORD_NMAX][30];       /* label of variable */
    char   Coord_Var[MONnD_COORD_NMAX][30]; /* short id of variable */
    int    Coord_Bin[MONnD_COORD_NMAX];             /* bins of variable array */
    double Coord_Min[MONnD_COORD_NMAX];             
    double Coord_Max[MONnD_COORD_NMAX];
    char   Monitor_Label[MONnD_COORD_NMAX*30];      /* Label for monitor */
    char   Mon_File[128];    /* output file name */

    double cx,cy,cz;
    double cvx, cvy, cvz;
    double csx, csy, csz;
    double cs1, cs2, ct, cp;
    double He3_pressure;
    char   Flag_UsePreMonitor    ;   /* use a previously stored neutron parameter set */
    char   UserName1[128];
    char   UserName2[128];
    double UserVariable1;
    double UserVariable2;
    double Intermediate;
    double IntermediateCnts;
    char   option[1024];
    
    int    Nsum;
    double psum, p2sum;
    int    **Mon2D_N;
    double **Mon2D_p;
    double **Mon2D_p2;
    double *Mon2D_Buffer;
    
    double mxmin,mxmax,mymin,mymax,mzmin,mzmax;
    
    char   compcurname[128];

  } MonitornD_Variables_type;
  
/* global functions to be defined once for faster compile */
/* code can anyway get big. may be included in an external .c file 
   to be included once */
  
/* this routine is used to save results at end of simulation, but also 
 * during simulation (SIGUSR... intermediate steps) */ 
void Monitor_nD_OutPut(aDEFS, aVars, dofree)
  MonitornD_Defines_type *aDEFS;
  MonitornD_Variables_type *aVars;
  char dofree;
  /* dofree = 0 : no free, =1 : free variables (last call) */
  {
    char   *fname;
    long    i,j;
    int    *p0m = NULL;
    double *p1m = NULL;
    double *p2m = NULL;
    char    Coord_X_Label[1024];
    double  min1d, max1d; 
    double  min2d, max2d;
    int     bin1d, bin2d;
    char    While_End = 0;
    long    While_Buffer = 0;
    double  XY, pp;
    double  Coord[MONnD_COORD_NMAX];
    long    Coord_Index[MONnD_COORD_NMAX];
    char    label[1024];
    
    if (dofree == 0)
    {
      if (aVars->Flag_Verbose) printf("Monitor_nD: %s save intermediate results (%.2f %%).\n", aVars->compcurname, 100*mcget_run_num()/mcget_ncount());
    }

    /* check Buffer flush when end of simulation reached */
    if ((aVars->Buffer_Counter < aVars->Buffer_Block) && aVars->Flag_Auto_Limits)
    {
      /* Get Auto Limits */
      if (aVars->Flag_Verbose) printf("Monitor_nD: %s getting %i Auto Limits from List (%li).\n", aVars->compcurname, aVars->Coord_Number, aVars->Buffer_Counter);
      for (i = 1; i <= aVars->Coord_Number; i++)
      {
        aVars->Coord_Min[i] = FLT_MAX;
        aVars->Coord_Max[i] = -FLT_MAX;

        for (j = 0; j < aVars->Buffer_Counter; j++)
        { 
                XY = aVars->Mon2D_Buffer[j*(aVars->Coord_Number+2) + (i-1)];  /* scanning variables in Buffer */
          if (XY < aVars->Coord_Min[i]) aVars->Coord_Min[i] = XY;
          if (XY > aVars->Coord_Max[i]) aVars->Coord_Max[i] = XY;

        }
      }
      aVars->Flag_Auto_Limits = 2;  /* pass to 2nd auto limits step */
      aVars->Buffer_Block = aVars->Buffer_Counter;
    
      while (!While_End)
      { /* we generate Coord[] and Coord_index[] from Buffer (auto limits) or passing neutron */
        if (While_Buffer < aVars->Buffer_Block)
        {
          /* first while loops (While_Buffer) */
          Coord[0] = aVars->Mon2D_Buffer[While_Buffer*(aVars->Coord_Number+2) + (aVars->Coord_Number)];
          pp = Coord[0];
            
          /* auto limits case : scan Buffer within limits and store in Mon2D */ 
          for (i = 1; i <= aVars->Coord_Number; i++)
          {
            /* scanning variables in Buffer */
            XY = (aVars->Coord_Max[i]-aVars->Coord_Min[i]);
            Coord[i] = aVars->Mon2D_Buffer[While_Buffer*(aVars->Coord_Number+2) + (i-1)];
            if (XY > 0) Coord_Index[i] = floor((aVars->Mon2D_Buffer[(i-1) + While_Buffer*(aVars->Coord_Number+2)]-aVars->Coord_Min[i])*aVars->Coord_Bin[i]/XY);
            else Coord_Index[i] = 0;
            if (aVars->Flag_With_Borders)
            {
              if (Coord_Index[i] < 0) Coord_Index[i] = 0;
              if (Coord_Index[i] >= aVars->Coord_Bin[i]) Coord_Index[i] = aVars->Coord_Bin[i] - 1;
            }
          } /* end for */
          While_Buffer++;
        } /* end if in Buffer */
        else /* (While_Buffer >= aVars->Buffer_Block) && (aVars->Flag_Auto_Limits == 2) */ 
        {
          aVars->Flag_Auto_Limits = 0;
          While_End = 1;
          if (!aVars->Flag_List && dofree) /* free Buffer not needed (no list to output) */
          { /* Dim : (aVars->Coord_Number+2)*aVars->Buffer_Block matrix (for p, dp) */ 
            free(aVars->Mon2D_Buffer);
          }
        }

        /* store n1d/2d section from Buffer */

        /* 1D and n1D case : aVars->Flag_Multiple */
        if (aVars->Flag_Multiple)
        { /* Dim : aVars->Coord_Number*aVars->Coord_Bin[i] vectors (intensity is not included) */ 
          for (i= 0; i < aVars->Coord_Number; i++)
          {
            j = Coord_Index[i+1];
            if (j >= 0 && j < aVars->Coord_Bin[i+1])
            {
              aVars->Mon2D_N[i][j]++;
              aVars->Mon2D_p[i][j] += pp;
              aVars->Mon2D_p2[i][j] += pp*pp;
            }
          }
        }
        else /* 2D case : aVars->Coord_Number==2 and !aVars->Flag_Multiple and !aVars->Flag_List */
        if ((aVars->Coord_Number == 2) && !aVars->Flag_Multiple)
        { /* Dim : aVars->Coord_Bin[1]*aVars->Coord_Bin[2] matrix */
          i = Coord_Index[1];
          j = Coord_Index[2];
          if (i >= 0 && i < aVars->Coord_Bin[1] && j >= 0 && j < aVars->Coord_Bin[2])
          {
            aVars->Mon2D_N[i][j]++;
            aVars->Mon2D_p[i][j] += pp;
            aVars->Mon2D_p2[i][j] += pp*pp;
          }
        } /* end store 2D/1D */
      } /* end while */
    } /* end Force Get Limits */

    if (aVars->Flag_Verbose) 
    {
      printf("Monitor_nD: %s is a %s.\n", aVars->compcurname, aVars->Monitor_Label);
      printf("Monitor_nD: version %s with options=%s\n", Monitor_nD_Version, aVars->option);
    }

    /* write oputput files (sent to file as p[i*n + j] vectors) */
    if (aVars->Coord_Number == 0) 
    {
      /* DETECTOR_OUT_0D(aVars->Monitor_Label, aVars->Nsum, aVars->psum, aVars->p2sum); */
      mcdetector_out_0D(aVars->Monitor_Label, aVars->Nsum, aVars->psum, aVars->p2sum, aVars->compcurname);
    }
    else
    if (strlen(aVars->Mon_File) > 0) 
    { 
      fname = (char*)malloc(strlen(aVars->Mon_File)+10*aVars->Coord_Number);
      if (aVars->Flag_List) /* List: DETECTOR_OUT_2D */
      {
        if (aVars->Flag_List == 2) aVars->Buffer_Size = aVars->Neutron_Counter;
        strcpy(fname,aVars->Mon_File);
        if (strchr(aVars->Mon_File,'.') == NULL) strcat(fname, "_list");

        min1d = 1; max1d = aVars->Coord_Number+2;
        min2d = 0; max2d = aVars->Buffer_Size; 
        bin1d = aVars->Coord_Number+2; bin2d = aVars->Buffer_Size;
        strcpy(Coord_X_Label,"");
        for (i= 1; i <= aVars->Coord_Number; i++)
        {
          if (min2d < aVars->Coord_Min[i]) min2d = aVars->Coord_Min[i];
          if (max2d < aVars->Coord_Max[i]) max2d = aVars->Coord_Max[i];
          strcat(Coord_X_Label, aVars->Coord_Var[i]);
          strcat(Coord_X_Label, " ");
          if (strchr(aVars->Mon_File,'.') == NULL)
                  { strcat(fname, "."); strcat(fname, aVars->Coord_Var[i]); }
        }
        strcat(Coord_X_Label, "I I_err");
        if (aVars->Flag_Verbose) printf("Monitor_nD: %s write monitor file %s List (%ix%li).\n", aVars->compcurname, fname,(aVars->Coord_Number+2),aVars->Buffer_Size);
        p0m = (int    *)malloc((aVars->Coord_Number+2)*aVars->Buffer_Size*sizeof(int));
        p1m = (double *)malloc((aVars->Coord_Number+2)*aVars->Buffer_Size*sizeof(double)); 
        if (min2d == max2d) max2d = min2d+1e-6;
        if (min1d == max1d) max1d = min1d+1e-6;
        if (dofree == 0)
        {
          sprintf(label, "%s (%.2f %%)", aVars->Monitor_Label, 100*mcget_run_num()/mcget_ncount());
        }
        else
          strcpy(label, aVars->Monitor_Label);
        if (p1m == NULL) /* use Raw Buffer line output */
        {
          if (aVars->Flag_Verbose) printf("Monitor_nD: %s cannot allocate memory for List transpose. Skipping.\n", aVars->compcurname);
          if (p0m != NULL) free(p0m);
          mcdetector_out_2D(
            label,
            aVars->Coord_Label[0],
            Coord_X_Label,
            min2d, max2d, 
            min1d, max1d, 
            bin2d,
            bin1d,
            (int *)aVars->Mon2D_Buffer,aVars->Mon2D_Buffer,aVars->Mon2D_Buffer,
            fname, aVars->compcurname); 
        
        }
        else /* transpose for column output */
        {
          for (i=0; i < aVars->Coord_Number+2; i++) 
            for (j=0; j < aVars->Buffer_Size; j++)
            {
              p1m[i*aVars->Buffer_Size+j] = aVars->Mon2D_Buffer[j*(aVars->Coord_Number+2) + i];
              p0m[i*aVars->Buffer_Size+j] = 1;
            }
            
          mcdetector_out_2D(
              label,
              Coord_X_Label,
              aVars->Coord_Label[0],
              min1d, max1d, 
              min2d, max2d, 
              bin1d,
              bin2d,
              p0m,p1m,aVars->Mon2D_Buffer,
              fname, aVars->compcurname); 
          free(p0m);
          free(p1m);
        }
      }
      if (aVars->Flag_Multiple) /* n1D: DETECTOR_OUT_1D */
      {
        for (i= 0; i < aVars->Coord_Number; i++)
        {
          
          strcpy(fname,aVars->Mon_File);
          if (strchr(aVars->Mon_File,'.') == NULL)
                  { strcat(fname, "."); strcat(fname, aVars->Coord_Var[i+1]); }
          sprintf(Coord_X_Label, "%s monitor", aVars->Coord_Label[i+1]);
          if (dofree == 0)
          {
            sprintf(label, "%s (%.2f %%)", Coord_X_Label, 100*mcget_run_num()/mcget_ncount());
          }
          else
            strcpy(label, Coord_X_Label);
          if (aVars->Flag_Verbose) printf("Monitor_nD: %s write monitor file %s 1D (%i).\n", aVars->compcurname, fname, aVars->Coord_Bin[i+1]);
          min1d = aVars->Coord_Min[i+1];
          max1d = aVars->Coord_Max[i+1];
          if (min1d == max1d) max1d = min1d+1e-6;
          p1m = (double *)malloc(aVars->Coord_Bin[i+1]*sizeof(double));
          p2m = (double *)malloc(aVars->Coord_Bin[i+1]*sizeof(double));
          if (p2m == NULL) /* use Raw Buffer line output */
          {
            if (aVars->Flag_Verbose) printf("Monitor_nD: %s cannot allocate memory for output. Using raw data.\n", aVars->compcurname);
            if (p1m != NULL) free(p1m);
            mcdetector_out_1D(
            label,
            aVars->Coord_Label[i+1],
            aVars->Coord_Label[0],
            aVars->Coord_Var[i+1],
            min1d, max1d, 
            aVars->Coord_Bin[i+1],
            aVars->Mon2D_N[i],aVars->Mon2D_p[i],aVars->Mon2D_p2[i],
            fname, aVars->compcurname); 
          }
          else
          {
            if (aVars->Flag_log != 0)
            {
              XY = FLT_MAX;
              for (j=0; j < aVars->Coord_Bin[i+1]; j++) /* search min of signal */
                if ((XY > aVars->Mon2D_p[i][j]) && (aVars->Mon2D_p[i][j] > 0)) XY = aVars->Mon2D_p[i][j];
              if (XY <= 0) XY = -log(FLT_MAX)/log(10); else XY = log(XY)/log(10)-1;
            }
            
            for (j=0; j < aVars->Coord_Bin[i+1]; j++)
            {
              if (aVars->Flag_log == 0)
              {
                p1m[j] = aVars->Mon2D_p[i][j];
                p2m[j] = aVars->Mon2D_p2[i][j];
              }
              else
              {
                if (aVars->Mon2D_p[i][j] > 0)
                {
                  p1m[j] = log(aVars->Mon2D_p[i][j])/log(10);
                  p2m[j] = aVars->Mon2D_p2[i][j]/aVars->Mon2D_p[i][j];
                }
                else
                {
                  p1m[j] = XY;
                  p2m[j] = 0;
                }
              }
            }
            mcdetector_out_1D(
              label,
              aVars->Coord_Label[i+1],
              aVars->Coord_Label[0],
              aVars->Coord_Var[i+1],
              min1d, max1d, 
              aVars->Coord_Bin[i+1],
              aVars->Mon2D_N[i],p1m,p2m,
              fname, aVars->compcurname); 
            
          }
          if (p1m != NULL) free(p1m);
          if (p2m != NULL) free(p2m);
            
        }
      }
      else
      if (aVars->Coord_Number == 2)  /* 2D: DETECTOR_OUT_2D */
      {
        strcpy(fname,aVars->Mon_File);

        p0m = (int    *)malloc(aVars->Coord_Bin[1]*aVars->Coord_Bin[2]*sizeof(int));
        p1m = (double *)malloc(aVars->Coord_Bin[1]*aVars->Coord_Bin[2]*sizeof(double));
        p2m = (double *)malloc(aVars->Coord_Bin[1]*aVars->Coord_Bin[2]*sizeof(double));
        if (p2m == NULL) 
        {
          if (aVars->Flag_Verbose) printf("Monitor_nD: %s cannot allocate memory for 2D array (%i). Skipping.\n", aVars->compcurname, 3*aVars->Coord_Bin[1]*aVars->Coord_Bin[2]*sizeof(double));
          if (p0m != NULL) free(p0m);
          if (p1m != NULL) free(p1m);
        }
        else
        {
          if (aVars->Flag_log != 0)
          {
            XY = FLT_MAX;
            for (i= 0; i < aVars->Coord_Bin[1]; i++)
              for (j= 0; j < aVars->Coord_Bin[2]; j++) /* search min of signal */
                if ((XY > aVars->Mon2D_p[i][j]) && (aVars->Mon2D_p[i][j]>0)) XY = aVars->Mon2D_p[i][j];
            if (XY <= 0) XY = -log(FLT_MAX)/log(10); else XY = log(XY)/log(10)-1;
          }
          for (i= 0; i < aVars->Coord_Bin[1]; i++)
          {
            for (j= 0; j < aVars->Coord_Bin[2]; j++)
            {
              p0m[j + i*aVars->Coord_Bin[2]] = aVars->Mon2D_N[i][j];
              if (aVars->Flag_log == 0)
              {
                p1m[j + i*aVars->Coord_Bin[2]] = aVars->Mon2D_p[i][j];
                p2m[j + i*aVars->Coord_Bin[2]] = aVars->Mon2D_p2[i][j];
              }
              else
              {
                if (aVars->Mon2D_p[i][j] > 0)
                {
                  p1m[j + i*aVars->Coord_Bin[2]] = log(aVars->Mon2D_p[i][j])/log(10);
                  p2m[j + i*aVars->Coord_Bin[2]] = aVars->Mon2D_p2[i][j]/aVars->Mon2D_p[i][j];
                }
                else
                {
                  p1m[j + i*aVars->Coord_Bin[2]] = XY;
                  p2m[j + i*aVars->Coord_Bin[2]] = 0;
                }
              }
            }
          }
          if (strchr(aVars->Mon_File,'.') == NULL)
                  { strcat(fname, "."); strcat(fname, aVars->Coord_Var[1]);
              strcat(fname, "_"); strcat(fname, aVars->Coord_Var[2]); }
          if (aVars->Flag_Verbose) printf("Monitor_nD: %s write monitor file %s 2D (%ix%i).\n", aVars->compcurname, fname, aVars->Coord_Bin[1], aVars->Coord_Bin[2]);

          min1d = aVars->Coord_Min[1];
          max1d = aVars->Coord_Max[1];
          if (min1d == max1d) max1d = min1d+1e-6;
          min2d = aVars->Coord_Min[2];
          max2d = aVars->Coord_Max[2];
          if (min2d == max2d) max2d = min2d+1e-6;
          if (dofree == 0)
          {
            sprintf(label, "%s (%.2f %%)", aVars->Monitor_Label, 100*mcget_run_num()/mcget_ncount());
          }
          else
            strcpy(label, aVars->Monitor_Label);

          mcdetector_out_2D(
            label,
            aVars->Coord_Label[1],
            aVars->Coord_Label[2],
            min1d, max1d,  
            min2d, max2d,  
            aVars->Coord_Bin[1],
            aVars->Coord_Bin[2],
            p0m,p1m,p2m,
            fname, aVars->compcurname);

          if (p0m != NULL) free(p0m);
          if (p1m != NULL) free(p1m);
          if (p2m != NULL) free(p2m); 
        }
      }
      free(fname);
    }

    /* Now Free memory Mon2D.. */
    if ((aVars->Flag_Auto_Limits || aVars->Flag_List) && aVars->Coord_Number)
    { /* Dim : (aVars->Coord_Number+2)*aVars->Buffer_Block matrix (for p, dp) */ 
      if (aVars->Mon2D_Buffer != NULL && dofree) free(aVars->Mon2D_Buffer);
    }
       
    /* 1D and n1D case : aVars->Flag_Multiple */
    if (aVars->Flag_Multiple && aVars->Coord_Number && dofree)
    { /* Dim : aVars->Coord_Number*aVars->Coord_Bin[i] vectors */ 
      for (i= 0; i < aVars->Coord_Number; i++)
      { 
        free(aVars->Mon2D_N[i]);
        free(aVars->Mon2D_p[i]);
        free(aVars->Mon2D_p2[i]);
      }
      free(aVars->Mon2D_N);
      free(aVars->Mon2D_p);
      free(aVars->Mon2D_p2);
    }
    

    /* 2D case : aVars->Coord_Number==2 and !aVars->Flag_Multiple and !aVars->Flag_List */
    if ((aVars->Coord_Number == 2) && !aVars->Flag_Multiple && dofree)
    { /* Dim : aVars->Coord_Bin[1]*aVars->Coord_Bin[2] matrix */
      for (i= 0; i < aVars->Coord_Bin[1]; i++)
      { 
        free(aVars->Mon2D_N[i]);
        free(aVars->Mon2D_p[i]);
        free(aVars->Mon2D_p2[i]);
      }
      free(aVars->Mon2D_N); 
      free(aVars->Mon2D_p);
      free(aVars->Mon2D_p2);
    }
  }
  
void Monitor_nD_Init(aDEFS, aVars, m_xwidth, m_yheight, m_zthick, m_xmin, m_xmax, m_ymin, m_ymax, m_zmin, m_zmax)
  MonitornD_Defines_type *aDEFS;
  MonitornD_Variables_type *aVars;
  double m_xwidth, m_yheight, m_zthick;
  double m_xmin, m_xmax, m_ymin, m_ymax, m_zmin, m_zmax;
  {
    long carg = 1;
    char *option_copy, *token;
    char Flag_New_Token = 1;
    char Flag_End       = 1;
    char Flag_All       = 0;
    char Flag_No        = 0;
    char Set_aVars_Coord_Type;
    char Set_Coord_Flag = 0;
    char Set_aVars_Coord_Label[30];
    char Set_aVars_Coord_Var[30];
    char Short_Label[MONnD_COORD_NMAX][30];
    char Set_Coord_Mode;
    long i=0, j=0;
    double lmin, lmax, XY;
    long    t;
    
    t = (long)time(NULL);

    aDEFS->COORD_NONE   =0;
    aDEFS->COORD_X      =1;
    aDEFS->COORD_Y      =2;
    aDEFS->COORD_Z      =3;
    aDEFS->COORD_VX     =4;
    aDEFS->COORD_VY     =5;
    aDEFS->COORD_VZ     =6;
    aDEFS->COORD_T      =7;
    aDEFS->COORD_P      =8;
    aDEFS->COORD_SX     =9;
    aDEFS->COORD_SY     =10;
    aDEFS->COORD_SZ     =11;
    aDEFS->COORD_KX     =12;
    aDEFS->COORD_KY     =13;
    aDEFS->COORD_KZ     =14;
    aDEFS->COORD_K      =15;
    aDEFS->COORD_V      =16;
    aDEFS->COORD_ENERGY =17;
    aDEFS->COORD_LAMBDA =18;
    aDEFS->COORD_RADIUS =19;
    aDEFS->COORD_HDIV   =20;
    aDEFS->COORD_VDIV   =21;
    aDEFS->COORD_ANGLE  =22;
    aDEFS->COORD_NCOUNT =23;
    aDEFS->COORD_THETA  =24;
    aDEFS->COORD_PHI    =25;
    aDEFS->COORD_USER1  =26;
    aDEFS->COORD_USER2  =26;

/* token modifiers */
    aDEFS->COORD_VAR    =0;    /* next token should be a variable or normal option */
    aDEFS->COORD_MIN    =1;    /* next token is a min value */
    aDEFS->COORD_MAX    =2;    /* next token is a max value */
    aDEFS->COORD_DIM    =3;    /* next token is a bin value */
    aDEFS->COORD_FIL    =4;    /* next token is a filename */
    aDEFS->COORD_EVNT   =5;    /* next token is a buffer size value */
    aDEFS->COORD_3HE    =6;    /* next token is a 3He pressure value */
    aDEFS->COORD_INTERM =7;    /* next token is an intermediate save value (%) */

    strcpy(aDEFS->TOKEN_DEL, " =,;[](){}:");  /* token separators */

    aDEFS->SHAPE_SQUARE =0;    /* shape of the monitor */
    aDEFS->SHAPE_DISK   =1;
    aDEFS->SHAPE_SPHERE =2;
    aDEFS->SHAPE_CYLIND =3;
    aDEFS->SHAPE_BOX    =4;

    aVars->Sphere_Radius     = 0;
    aVars->Cylinder_Height   = 0;
    aVars->Flag_With_Borders = 0;   /* 2 means xy borders too */
    aVars->Flag_List         = 0;   /* 1 store 1 buffer, 2 is list all */
    aVars->Flag_Multiple     = 0;   /* 1 when n1D, 0 for 2D */
    aVars->Flag_Verbose      = 0;
    aVars->Flag_Shape        = aDEFS->SHAPE_SQUARE;
    aVars->Flag_Auto_Limits  = 0;   /* get limits from first Buffer */
    aVars->Flag_Absorb       = 0;   /* monitor is also a slit */
    aVars->Flag_per_cm2      = 0;   /* flux is per cm2 */
    aVars->Flag_log          = 0;   /* log10 of the flux */
    aVars->Flag_parallel     = 0;   /* set neutron state back after detection (parallel components) */
    aVars->Coord_Number      = 0;   /* total number of variables to monitor, plus intensity (0) */
    aVars->Buffer_Block      = 1000;     /* Buffer size for list or auto limits */
    aVars->Neutron_Counter   = 0;   /* event counter, simulation total counts is mcget_ncount() */
    aVars->Buffer_Counter    = 0;   /* index in Buffer size (for realloc) */
    aVars->Buffer_Size       = 0;
    aVars->UserVariable1     = 0;
    aVars->UserVariable2     = 0;
    aVars->He3_pressure      = 0;
    aVars->IntermediateCnts  = 0;
    
    Set_aVars_Coord_Type = aDEFS->COORD_NONE;
    Set_Coord_Mode = aDEFS->COORD_VAR;
    
    /* handle size parameters */
    /* normal use is with xwidth, yheight, zthick */
    /* if xmin,xmax,ymin,ymax,zmin,zmax are non 0, use them */
    if (fabs(m_xmin-m_xmax) == 0) 
      { aVars->mxmin = -fabs(m_xwidth)/2; aVars->mxmax = fabs(m_xwidth)/2; }
    else
      { if (m_xmin < m_xmax) {aVars->mxmin = m_xmin; aVars->mxmax = m_xmax;}
        else {aVars->mxmin = m_xmax; aVars->mxmax = m_xmin;}
      }
    if (fabs(m_ymin-m_ymax) == 0) 
      { aVars->mymin = -fabs(m_yheight)/2; aVars->mymax = fabs(m_yheight)/2; }
    else
      { if (m_ymin < m_ymax) {aVars->mymin = m_ymin; aVars->mymax = m_ymax;}
        else {aVars->mymin = m_ymax; aVars->mymax = m_ymin;}
      }
    if (fabs(m_zmin-m_zmax) == 0) 
      { aVars->mzmin = -fabs(m_zthick)/2; aVars->mzmax = fabs(m_zthick)/2; }
    else
      { if (m_zmin < m_zmax) {aVars->mzmin = m_zmin; aVars->mzmax = m_zmax; }
        else {aVars->mzmin = m_zmax; aVars->mzmax = m_zmin; }
      }
      
    if (fabs(zmax-zmin) == 0)
      aVars->Flag_Shape        = aDEFS->SHAPE_SQUARE;
    else
      aVars->Flag_Shape        = aDEFS->SHAPE_BOX;
      
    /* parse option string */ 
    
    option_copy = (char*)malloc(strlen(aVars->option));
    if (option_copy == NULL)
    {
      printf("Monitor_nD: %s cannot allocate option_copy (%i). Fatal.\n", aVars->compcurname, strlen(aVars->option));
      exit(-1);
    }
    
    
    if (strlen(aVars->option))
    {
      Flag_End = 0;
      strcpy(option_copy, aVars->option);
    }
    
    if (strstr(aVars->option, "cm2") || strstr(aVars->option, "cm^2")) aVars->Flag_per_cm2 = 1;
  
    if (aVars->Flag_per_cm2) strcpy(aVars->Coord_Label[0],"Intensity [n/cm^2/s]");
    else strcpy(aVars->Coord_Label[0],"Intensity [n/s]");
    strcpy(aVars->Coord_Var[0],"p");
    aVars->Coord_Type[0] = aDEFS->COORD_P;
    aVars->Coord_Bin[0] = 1;
    aVars->Coord_Min[0] = 0;
    aVars->Coord_Max[0] = FLT_MAX;
    
    /* default file name is comp name+dateID */
    sprintf(aVars->Mon_File, "%s_%i", aVars->compcurname, t); 
  
    carg = 1;
    while((Flag_End == 0) && (carg < 128))
    {
      if (Flag_New_Token) /* to get the previous token sometimes */
      {
        if (carg == 1) token=(char *)strtok(option_copy,aDEFS->TOKEN_DEL);
        else token=(char *)strtok(NULL,aDEFS->TOKEN_DEL);
        if (token == NULL) Flag_End=1;
      }
      Flag_New_Token = 1;
      if ((token != NULL) && (strlen(token) != 0))
      {
      /* first handle option values from preceeding keyword token detected */
        if (Set_Coord_Mode == aDEFS->COORD_MAX)
        { 
          if (!Flag_All) 
            aVars->Coord_Max[aVars->Coord_Number] = atof(token); 
          else
            for (i = 0; i <= aVars->Coord_Number; aVars->Coord_Max[i++] = atof(token));
          Set_Coord_Mode = aDEFS->COORD_VAR; Flag_All = 0;
        }
        if (Set_Coord_Mode == aDEFS->COORD_MIN)
        { 
          if (!Flag_All) 
            aVars->Coord_Min[aVars->Coord_Number] = atof(token); 
          else
            for (i = 0; i <= aVars->Coord_Number; aVars->Coord_Min[i++] = atof(token));
          Set_Coord_Mode = aDEFS->COORD_MAX; 
        }
        if (Set_Coord_Mode == aDEFS->COORD_DIM)
        { 
          if (!Flag_All) 
            aVars->Coord_Bin[aVars->Coord_Number] = atoi(token); 
          else
            for (i = 0; i <= aVars->Coord_Number; aVars->Coord_Bin[i++] = atoi(token));
          Set_Coord_Mode = aDEFS->COORD_VAR; Flag_All = 0;
        }
        if (Set_Coord_Mode == aDEFS->COORD_FIL)
        { 
          if (!Flag_No) strcpy(aVars->Mon_File,token); 
          else { strcpy(aVars->Mon_File,""); aVars->Coord_Number = 0; Flag_End = 1;}
          Set_Coord_Mode = aDEFS->COORD_VAR;
        }
        if (Set_Coord_Mode == aDEFS->COORD_EVNT)
        { 
          if (!strcmp(token, "all") || Flag_All) aVars->Flag_List = 2;
          else { i = atoi(token); if (i) aVars->Buffer_Counter = i; }
          Set_Coord_Mode = aDEFS->COORD_VAR; Flag_All = 0;
        }
        if (Set_Coord_Mode == aDEFS->COORD_3HE)
        { 
            aVars->He3_pressure = atof(token); 
            Set_Coord_Mode = aDEFS->COORD_VAR; Flag_All = 0; 
        }
        if (Set_Coord_Mode == aDEFS->COORD_INTERM)
        { 
            aVars->Intermediate = atof(token); 
            Set_Coord_Mode = aDEFS->COORD_VAR; Flag_All = 0; 
        }

        /* now look for general option keywords */
        if (!strcmp(token, "borders")) 
        { if (Flag_No) { aVars->Flag_With_Borders = 0; Flag_No = 0; }
          else aVars->Flag_With_Borders = 1; }
        if (!strcmp(token, "verbose")) 
        { if (Flag_No) { aVars->Flag_Verbose = 0; Flag_No = 0; }
          else aVars->Flag_Verbose      = 1; }
        if (!strcmp(token, "log")) 
        { if (Flag_No) { aVars->Flag_log = 0; Flag_No = 0; }
          else aVars->Flag_log      = 1; }
        if (!strcmp(token, "multiple")) 
        { if (Flag_No) { aVars->Flag_Multiple = 0; Flag_No = 0; }
          else aVars->Flag_Multiple = 1; }
        if (!strcmp(token, "list")) 
        { if (Flag_No) { aVars->Flag_List = 0; Flag_No = 0; }
          else aVars->Flag_List = 1; 
          Set_Coord_Mode = aDEFS->COORD_EVNT; }

        if (!strcmp(token, "limits") || !strcmp(token, "min")) Set_Coord_Mode = aDEFS->COORD_MIN;
        if (!strcmp(token, "slit") || !strcmp(token, "absorb")) 
        { if (Flag_No) { aVars->Flag_Absorb = 0; Flag_No = 0; }
          else aVars->Flag_Absorb = 1; }
        if (!strcmp(token, "max")) Set_Coord_Mode = aDEFS->COORD_MAX;
        if (!strcmp(token, "bins")) Set_Coord_Mode = aDEFS->COORD_DIM;
        if (!strcmp(token, "file")) 
        { Set_Coord_Mode = aDEFS->COORD_FIL;
          if (Flag_No) { strcpy(aVars->Mon_File,""); aVars->Coord_Number = 0; Flag_End = 1;}}
        if (!strcmp(token, "unactivate")) { Flag_End = 1; aVars->Coord_Number = 0; }
        if (!strcmp(token, "all")) 
        { if (Flag_No) { Flag_All = 0; Flag_No = 0; }
          else Flag_All = 1; }
        if (!strcmp(token, "sphere")) 
        { if (Flag_No) { aVars->Flag_Shape = aDEFS->SHAPE_SQUARE; Flag_No = 0; }
          else aVars->Flag_Shape = aDEFS->SHAPE_SPHERE; }
        if (!strcmp(token, "cylinder")) 
        { if (Flag_No) { aVars->Flag_Shape = aDEFS->SHAPE_SQUARE; Flag_No = 0; }
          else aVars->Flag_Shape = aDEFS->SHAPE_CYLIND; }
        if (!strcmp(token, "square")) aVars->Flag_Shape = aDEFS->SHAPE_SQUARE; 
        if (!strcmp(token, "disk")) aVars->Flag_Shape = aDEFS->SHAPE_DISK; 
        if (!strcmp(token, "box")) aVars->Flag_Shape = aDEFS->SHAPE_BOX; 
        if (!strcmp(token, "parallel")) aVars->Flag_parallel = 1; 
        if (!strcmp(token, "auto")) 
        { if (Flag_No) { aVars->Flag_Auto_Limits = 0; Flag_No = 0; }
          else aVars->Flag_Auto_Limits = 1; }
        if (!strcmp(token, "premonitor"))
        { if (Flag_No) { aVars->Flag_UsePreMonitor = 0; Flag_No = 0; }
        else aVars->Flag_UsePreMonitor == 1; } 
        if (!strcmp(token, "3He_pressure"))
        { if (!Flag_No)  Set_Coord_Mode = aDEFS->COORD_3HE; 
          aVars->He3_pressure = 3; }
        if (!strcmp(token, "intermediate"))
        { if (!Flag_No)  Set_Coord_Mode = aDEFS->COORD_INTERM; 
          aVars->Intermediate = 5; }
        if (!strcmp(token, "no") || !strcmp(token, "not")) Flag_No = 1;
  
        /* now look for variable names to monitor */
        Set_aVars_Coord_Type = aDEFS->COORD_NONE; lmin = 0; lmax = 0;

        if (!strcmp(token, "x")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_X; strcpy(Set_aVars_Coord_Label,"x [m]"); strcpy(Set_aVars_Coord_Var,"x"); lmin = aVars->mxmin; lmax = aVars->mxmax; }
        if (!strcmp(token, "y")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_Y; strcpy(Set_aVars_Coord_Label,"y [m]"); strcpy(Set_aVars_Coord_Var,"y"); lmin = aVars->mymin; lmax = aVars->mymax; }
        if (!strcmp(token, "z")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_Z; strcpy(Set_aVars_Coord_Label,"z [m]"); strcpy(Set_aVars_Coord_Var,"z"); lmin = aVars->mzmin; lmax = aVars->mzmax; }
        if (!strcmp(token, "k") || !strcmp(token, "wavevector")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_K; strcpy(Set_aVars_Coord_Label,"|k| [Angs-1]"); strcpy(Set_aVars_Coord_Var,"k"); lmin = 0; lmax = 10; }
        if (!strcmp(token, "v"))
          { Set_aVars_Coord_Type = aDEFS->COORD_V; strcpy(Set_aVars_Coord_Label,"Velocity [m/s]"); strcpy(Set_aVars_Coord_Var,"v"); lmin = 0; lmax = 10000; }
        if (!strcmp(token, "t") || !strcmp(token, "time")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_T; strcpy(Set_aVars_Coord_Label,"TOF [s]"); strcpy(Set_aVars_Coord_Var,"t"); lmin = 0; lmax = .1; }
        if ((aVars->Coord_Number > 0) && (!strcmp(token, "p") || !strcmp(token, "intensity") || !strcmp(token, "flux")))
          { Set_aVars_Coord_Type = aDEFS->COORD_P;  
            if (aVars->Flag_per_cm2) strcpy(Set_aVars_Coord_Label,"Intensity [n/cm^2/s]");
            else strcpy(Set_aVars_Coord_Label,"Intensity [n/s]"); 
            strcpy(Set_aVars_Coord_Var,"I"); 
            lmin = 0; lmax = FLT_MAX; }

        if (!strcmp(token, "vx")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_VX; strcpy(Set_aVars_Coord_Label,"vx [m/s]"); strcpy(Set_aVars_Coord_Var,"vx"); lmin = -1000; lmax = 1000; }
        if (!strcmp(token, "vy")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_VY; strcpy(Set_aVars_Coord_Label,"vy [m/s]"); strcpy(Set_aVars_Coord_Var,"vy"); lmin = -1000; lmax = 1000; }
        if (!strcmp(token, "vz")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_VZ; strcpy(Set_aVars_Coord_Label,"vz [m/s]"); strcpy(Set_aVars_Coord_Var,"vz"); lmin = -10000; lmax = 10000; }
        if (!strcmp(token, "kx")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_KX; strcpy(Set_aVars_Coord_Label,"kx [Angs-1]"); strcpy(Set_aVars_Coord_Var,"kx"); lmin = -1; lmax = 1; }
        if (!strcmp(token, "ky")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_KY; strcpy(Set_aVars_Coord_Label,"ky [Angs-1]"); strcpy(Set_aVars_Coord_Var,"ky"); lmin = -1; lmax = 1; }
        if (!strcmp(token, "kz")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_KZ; strcpy(Set_aVars_Coord_Label,"kz [Angs-1]"); strcpy(Set_aVars_Coord_Var,"kz"); lmin = -10; lmax = 10; }
        if (!strcmp(token, "sx")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_SX; strcpy(Set_aVars_Coord_Label,"sx [1]"); strcpy(Set_aVars_Coord_Var,"sx"); lmin = -1; lmax = 1; }
        if (!strcmp(token, "sy")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_SY; strcpy(Set_aVars_Coord_Label,"sy [1]"); strcpy(Set_aVars_Coord_Var,"sy"); lmin = -1; lmax = 1; }
        if (!strcmp(token, "sz")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_SZ; strcpy(Set_aVars_Coord_Label,"sz [1]"); strcpy(Set_aVars_Coord_Var,"sz"); lmin = -1; lmax = 1; }

        if (!strcmp(token, "energy") || !strcmp(token, "omega")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_ENERGY; strcpy(Set_aVars_Coord_Label,"Energy [meV]"); strcpy(Set_aVars_Coord_Var,"E"); lmin = 0; lmax = 100; }
        if (!strcmp(token, "lambda") || !strcmp(token, "wavelength")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_LAMBDA; strcpy(Set_aVars_Coord_Label,"Wavelength [Angs]"); strcpy(Set_aVars_Coord_Var,"L"); lmin = 0; lmax = 100; }
        if (!strcmp(token, "radius")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_RADIUS; strcpy(Set_aVars_Coord_Label,"Radius [m]"); strcpy(Set_aVars_Coord_Var,"R"); lmin = 0; lmax = m_xmax; }
        if (!strcmp(token, "angle")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_ANGLE; strcpy(Set_aVars_Coord_Label,"Angle [deg]"); strcpy(Set_aVars_Coord_Var,"A"); lmin = -5; lmax = 5; }
        if (!strcmp(token, "hdiv")|| !strcmp(token, "divergence") || !strcmp(token, "xdiv") || !strcmp(token, "dx")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_HDIV; strcpy(Set_aVars_Coord_Label,"Hor. Divergence [deg]"); strcpy(Set_aVars_Coord_Var,"HD"); lmin = -5; lmax = 5; }
        if (!strcmp(token, "vdiv") || !strcmp(token, "ydiv") || !strcmp(token, "dy")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_VDIV; strcpy(Set_aVars_Coord_Label,"Vert. Divergence [deg]"); strcpy(Set_aVars_Coord_Var,"VD"); lmin = -5; lmax = 5; }
        if (!strcmp(token, "theta") || !strcmp(token, "longitude")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_THETA; strcpy(Set_aVars_Coord_Label,"Longitude [deg]"); strcpy(Set_aVars_Coord_Var,"th"); lmin = -180; lmax = 180; }
        if (!strcmp(token, "phi") || !strcmp(token, "lattitude")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_PHI; strcpy(Set_aVars_Coord_Label,"Lattitude [deg]"); strcpy(Set_aVars_Coord_Var,"ph"); lmin = -180; lmax = 180; }
        if (!strcmp(token, "ncounts")) 
          { Set_aVars_Coord_Type = aDEFS->COORD_NCOUNT; strcpy(Set_aVars_Coord_Label,"Neutrons [1]"); strcpy(Set_aVars_Coord_Var,"N"); lmin = 0; lmax = 1e10; }
              if (!strcmp(token, "user") || !strcmp(token, "user1"))
          { Set_aVars_Coord_Type = aDEFS->COORD_USER1; strncpy(Set_aVars_Coord_Label,aVars->UserName1,32); strcpy(Set_aVars_Coord_Var,"U1"); lmin = -1e10; lmax = 1e10; }
              if (!strcmp(token, "user2"))
          { Set_aVars_Coord_Type = aDEFS->COORD_USER2; strncpy(Set_aVars_Coord_Label,aVars->UserName2,32); strcpy(Set_aVars_Coord_Var,"U2"); lmin = -1e10; lmax = 1e10; }

        /* now stores variable keywords detected, if any */ 
        if (Set_aVars_Coord_Type != aDEFS->COORD_NONE)
        {
          if (aVars->Coord_Number < MONnD_COORD_NMAX) aVars->Coord_Number++;
          else if (aVars->Flag_Verbose) printf("Monitor_nD: %s reached max number of variables (%i).\n", aVars->compcurname, MONnD_COORD_NMAX);
          aVars->Coord_Type[aVars->Coord_Number] = Set_aVars_Coord_Type;
          strcpy(aVars->Coord_Label[aVars->Coord_Number], Set_aVars_Coord_Label); 
          strcpy(aVars->Coord_Var[aVars->Coord_Number], Set_aVars_Coord_Var);
          if (lmin > lmax) { XY = lmin; lmin=lmax; lmax = XY; }
          aVars->Coord_Min[aVars->Coord_Number] = lmin;
          aVars->Coord_Max[aVars->Coord_Number] = lmax;
          aVars->Coord_Bin[aVars->Coord_Number] = 20;
          Set_Coord_Mode = aDEFS->COORD_VAR;
          Flag_All = 0;
          Flag_No = 0;
        }
      carg++;
      } /* end if token */
    } /* end while carg */
    free(option_copy);
    if (carg == 128) printf("Monitor_nD: %s reached max number of tokens (%i). Skipping.\n", aVars->compcurname, 128);
    
    if ((aVars->Flag_Shape == aDEFS->SHAPE_BOX) && (fabs(aVars->mzmax - aVars->mzmin) == 0)) aVars->Flag_Shape = aDEFS->SHAPE_SQUARE;

    /* now setting Monitor Name from variable labels */
    strcpy(aVars->Monitor_Label,"");
    for (i = 0; i <= aVars->Coord_Number; i++)
    {
      Set_aVars_Coord_Type = aVars->Coord_Type[i];
      if ((Set_aVars_Coord_Type == aDEFS->COORD_THETA)
       || (Set_aVars_Coord_Type == aDEFS->COORD_PHI)
       || (Set_aVars_Coord_Type == aDEFS->COORD_X)
       || (Set_aVars_Coord_Type == aDEFS->COORD_Y)
       || (Set_aVars_Coord_Type == aDEFS->COORD_Z)
       || (Set_aVars_Coord_Type == aDEFS->COORD_RADIUS))
       strcpy(Short_Label[i],"Position"); 
      else
      if ((Set_aVars_Coord_Type == aDEFS->COORD_VX)
       || (Set_aVars_Coord_Type == aDEFS->COORD_VY)
       || (Set_aVars_Coord_Type == aDEFS->COORD_VZ)
       || (Set_aVars_Coord_Type == aDEFS->COORD_V))
       strcpy(Short_Label[i],"Velocity"); 
      else
      if ((Set_aVars_Coord_Type == aDEFS->COORD_KX)
       || (Set_aVars_Coord_Type == aDEFS->COORD_KY)
       || (Set_aVars_Coord_Type == aDEFS->COORD_KZ)
       || (Set_aVars_Coord_Type == aDEFS->COORD_K))
       strcpy(Short_Label[i],"Wavevector"); 
      else
      if ((Set_aVars_Coord_Type == aDEFS->COORD_SX)
       || (Set_aVars_Coord_Type == aDEFS->COORD_SY)
       || (Set_aVars_Coord_Type == aDEFS->COORD_SZ))
       strcpy(Short_Label[i],"Spin");
      else
      if ((Set_aVars_Coord_Type == aDEFS->COORD_HDIV)
       || (Set_aVars_Coord_Type == aDEFS->COORD_VDIV)
       || (Set_aVars_Coord_Type == aDEFS->COORD_ANGLE))
       strcpy(Short_Label[i],"Divergence");
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_ENERGY)
       strcpy(Short_Label[i],"Energy"); 
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_LAMBDA)
       strcpy(Short_Label[i],"Wavelength"); 
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_NCOUNT)
       strcpy(Short_Label[i],"Neutron counts");
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_T)
          strcpy(Short_Label[i],"Time Of Flight");
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_P)
          strcpy(Short_Label[i],"Intensity");
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_USER1)
          strncpy(Short_Label[i],aVars->UserName1,32);
      else
      if (Set_aVars_Coord_Type == aDEFS->COORD_USER2)
          strncpy(Short_Label[i],aVars->UserName2,32);
      else
          strcpy(Short_Label[i],"Unknown"); 
          
      strcat(aVars->Monitor_Label, " ");
      strcat(aVars->Monitor_Label, Short_Label[i]);
    } /* end for Short_Label */
    strcat(aVars->Monitor_Label, " Monitor");
    if (aVars->Flag_Shape == aDEFS->SHAPE_SQUARE) strcat(aVars->Monitor_Label, " (Square)");
    if (aVars->Flag_Shape == aDEFS->SHAPE_DISK)   strcat(aVars->Monitor_Label, " (Disk)");
    if (aVars->Flag_Shape == aDEFS->SHAPE_SPHERE) strcat(aVars->Monitor_Label, " (Sphere)");
    if (aVars->Flag_Shape == aDEFS->SHAPE_CYLIND) strcat(aVars->Monitor_Label, " (Cylinder)");
    if (aVars->Flag_Shape == aDEFS->SHAPE_BOX)    strcat(aVars->Monitor_Label, " (Box)");
    if (((aVars->Flag_Shape == aDEFS->SHAPE_CYLIND) || (aVars->Flag_Shape == aDEFS->SHAPE_SPHERE) || (aVars->Flag_Shape == aDEFS->SHAPE_BOX))
        && strstr(aVars->option, "outgoing"))
    {
      aVars->Flag_Shape *= -1;
      strcat(aVars->Monitor_Label, " [out]");
    }
    if (aVars->Flag_UsePreMonitor == 1)
    {
        strcat(aVars->Monitor_Label, " at ");
        strncat(aVars->Monitor_Label, aVars->UserName1,32);
    }
    if (aVars->Flag_log == 1)
    {
        strcat(aVars->Monitor_Label, " (log) ");
    }
    
    /* aVars->Coord_Number  0   : intensity
     * aVars->Coord_Number  1:n : detector variables */
    
    /* now allocate memory to store variables in TRACE */
    if ((aVars->Coord_Number != 2) && !aVars->Flag_Multiple && !aVars->Flag_List) 
    { aVars->Flag_Multiple = 1; aVars->Flag_List = 0; } /* default is n1D */
    
   /* list and auto limits case : aVars->Flag_List or aVars->Flag_Auto_Limits 
    * -> Buffer to flush and suppress after aVars->Flag_Auto_Limits
    */
    if ((aVars->Flag_Auto_Limits || aVars->Flag_List) && aVars->Coord_Number)
    { /* Dim : (aVars->Coord_Number+2)*aVars->Buffer_Block matrix (for p, dp) */ 
      aVars->Mon2D_Buffer = (double *)malloc((aVars->Coord_Number+2)*aVars->Buffer_Block*sizeof(double));
      if (aVars->Mon2D_Buffer == NULL)
      { printf("Monitor_nD: %s cannot allocate aVars->Mon2D_Buffer (%li). No list and auto limits.\n", aVars->compcurname, aVars->Buffer_Block*(aVars->Coord_Number+2)*sizeof(double)); aVars->Flag_List = 0; aVars->Flag_Auto_Limits = 0; }
      else
      {
        for ( i=0; i < (aVars->Coord_Number+2)*aVars->Buffer_Block; aVars->Mon2D_Buffer[i++] = (double)0);
      }
      aVars->Buffer_Size = aVars->Buffer_Block;
    }
    
    /* 1D and n1D case : aVars->Flag_Multiple */
    if (aVars->Flag_Multiple && aVars->Coord_Number)
    { /* Dim : aVars->Coord_Number*aVars->Coord_Bin[i] vectors */ 
      aVars->Mon2D_N  = (int    **)malloc((aVars->Coord_Number)*sizeof(int *));
      aVars->Mon2D_p  = (double **)malloc((aVars->Coord_Number)*sizeof(double *));
      aVars->Mon2D_p2 = (double **)malloc((aVars->Coord_Number)*sizeof(double *));
      if ((aVars->Mon2D_N == NULL) || (aVars->Mon2D_p == NULL) || (aVars->Mon2D_p2 == NULL))
      { printf("Monitor_nD: %s n1D cannot allocate aVars->Mon2D_N/p/2p (%i). Fatal.\n", aVars->compcurname, (aVars->Coord_Number)*sizeof(double *)); exit(-1); }
      for (i= 1; i <= aVars->Coord_Number; i++)
      { 
        aVars->Mon2D_N[i-1]  = (int    *)malloc(aVars->Coord_Bin[i]*sizeof(int));
        aVars->Mon2D_p[i-1]  = (double *)malloc(aVars->Coord_Bin[i]*sizeof(double));
        aVars->Mon2D_p2[i-1] = (double *)malloc(aVars->Coord_Bin[i]*sizeof(double));
        if ((aVars->Mon2D_N == NULL) || (aVars->Mon2D_p == NULL) || (aVars->Mon2D_p2 == NULL))
        { printf("Monitor_nD: %s n1D cannot allocate %s aVars->Mon2D_N/p/2p[%li] (%i). Fatal.\n", aVars->compcurname, aVars->Coord_Var[i], i, (aVars->Coord_Bin[i])*sizeof(double *)); exit(-1); }
        else 
        {
          for (j=0; j < aVars->Coord_Bin[i]; j++ )
          { aVars->Mon2D_N[i-1][j] = (int)0; aVars->Mon2D_p[i-1][j] = (double)0; aVars->Mon2D_p2[i-1][j] = (double)0; }
        }
      }
    }
    else /* 2D case : aVars->Coord_Number==2 and !aVars->Flag_Multiple and !aVars->Flag_List */
    if ((aVars->Coord_Number == 2) && !aVars->Flag_Multiple)
    { /* Dim : aVars->Coord_Bin[1]*aVars->Coord_Bin[2] matrix */ 
      aVars->Mon2D_N  = (int    **)malloc((aVars->Coord_Bin[1])*sizeof(int *));
      aVars->Mon2D_p  = (double **)malloc((aVars->Coord_Bin[1])*sizeof(double *));
      aVars->Mon2D_p2 = (double **)malloc((aVars->Coord_Bin[1])*sizeof(double *));
      if ((aVars->Mon2D_N == NULL) || (aVars->Mon2D_p == NULL) || (aVars->Mon2D_p2 == NULL))
      { printf("Monitor_nD: %s 2D cannot allocate %s aVars->Mon2D_N/p/2p (%i). Fatal.\n", aVars->compcurname, aVars->Coord_Var[1], (aVars->Coord_Bin[1])*sizeof(double *)); exit(-1); }
      for (i= 0; i < aVars->Coord_Bin[1]; i++)
      { 
        aVars->Mon2D_N[i]  = (int    *)malloc(aVars->Coord_Bin[2]*sizeof(int));
        aVars->Mon2D_p[i]  = (double *)malloc(aVars->Coord_Bin[2]*sizeof(double));
        aVars->Mon2D_p2[i] = (double *)malloc(aVars->Coord_Bin[2]*sizeof(double));
        if ((aVars->Mon2D_N == NULL) || (aVars->Mon2D_p == NULL) || (aVars->Mon2D_p2 == NULL))
        { printf("Monitor_nD: %s 2D cannot allocate %s aVars->Mon2D_N/p/2p[%li] (%i). Fatal.\n", aVars->compcurname, aVars->Coord_Var[1], i, (aVars->Coord_Bin[2])*sizeof(double *)); exit(-1); }
        else
        {
          for (j=0; j < aVars->Coord_Bin[2]; j++ )
          { aVars->Mon2D_N[i][j] = (int)0; aVars->Mon2D_p[i][j] = (double)0; aVars->Mon2D_p2[i][j] = (double)0; }
        }
      }
    }      
      /* no Mon2D allocated for 
       * (aVars->Coord_Number != 2) && !aVars->Flag_Multiple && aVars->Flag_List */
    
    aVars->psum  = 0;
    aVars->p2sum = 0;
    aVars->Nsum  = 0;
    
    aVars->area  = fabs(aVars->mxmax - aVars->mxmin)*fabs(aVars->mymax - aVars->mymin)*1E4; /* in cm**2 for square and box shapes */
    aVars->Sphere_Radius = fabs(aVars->mxmax - aVars->mxmin);
    if ((abs(aVars->Flag_Shape) == aDEFS->SHAPE_DISK) || (abs(aVars->Flag_Shape) == aDEFS->SHAPE_SPHERE))
    {
      aVars->area = PI*aVars->Sphere_Radius*aVars->Sphere_Radius; /* disk shapes */
    }
    if (aVars->area == 0) aVars->Coord_Number = 0;
    if (aVars->Coord_Number == 0 && aVars->Flag_Verbose)  printf("Monitor_nD: %s is unactivated (0D)\n", aVars->compcurname);
    aVars->Cylinder_Height = fabs(aVars->mymax - aVars->mymin);
    
    if (aVars->Intermediate < 0) aVars->Intermediate = 0;
    if (aVars->Intermediate > 1) aVars->Intermediate /= 100;
    aVars->IntermediateCnts = aVars->Intermediate*mcget_ncount();
  } 
  
void Monitor_nD_Trace(aDEFS, aVars)
  MonitornD_Defines_type *aDEFS;
  MonitornD_Variables_type *aVars;
{

  double  XY=0;
  double  t0 = 0;
  double  t1 = 0;
  long    i,j;
  double  pp;
  double  Coord[MONnD_COORD_NMAX];
  long    Coord_Index[MONnD_COORD_NMAX];
  char    While_End   =0;
  long    While_Buffer=0;
  int     intersect = 0;
  char    Set_aVars_Coord_Type = aDEFS->COORD_NONE;
  
  
  pp = aVars->cp;
  
  if (aVars->Coord_Number > 0)
    {
      /* aVars->Flag_Auto_Limits */
      if ((aVars->Buffer_Counter >= aVars->Buffer_Block) && (aVars->Flag_Auto_Limits == 1))
      {
        /* auto limits case : get limits in Buffer for each variable */ 
              /* Dim : (aVars->Coord_Number+2)*aVars->Buffer_Block matrix (for p, dp) */ 
        if (aVars->Flag_Verbose) printf("Monitor_nD: %s getting %i Auto Limits from List (%li).\n", aVars->compcurname, aVars->Coord_Number, aVars->Buffer_Counter);
        for (i = 1; i <= aVars->Coord_Number; i++)
        {
          aVars->Coord_Min[i] = FLT_MAX;
          aVars->Coord_Max[i] = -FLT_MAX;
          for (j = 0; j < aVars->Buffer_Block; j++)
          { 
            XY = aVars->Mon2D_Buffer[j*(aVars->Coord_Number+2) + (i-1)];  /* scanning variables in Buffer */
            if (XY < aVars->Coord_Min[i]) aVars->Coord_Min[i] = XY;
            if (XY > aVars->Coord_Max[i]) aVars->Coord_Max[i] = XY;
          }
        }
        aVars->Flag_Auto_Limits = 2;  /* pass to 2nd auto limits step */
      }

      /* manage realloc for list all if Buffer size exceeded */
      if ((aVars->Buffer_Counter >= aVars->Buffer_Block) && (aVars->Flag_List == 2))
      {
        aVars->Mon2D_Buffer  = (double *)realloc(aVars->Mon2D_Buffer, (aVars->Coord_Number+2)*(aVars->Neutron_Counter+aVars->Buffer_Block)*sizeof(double));
        if (aVars->Mon2D_Buffer == NULL)
              { printf("Monitor_nD: %s cannot reallocate aVars->Mon2D_Buffer[%li] (%li). Skipping.\n", aVars->compcurname, i, (aVars->Neutron_Counter+aVars->Buffer_Block)*sizeof(double)); aVars->Flag_List = 1; }
        else { aVars->Buffer_Counter = 0; aVars->Buffer_Size = aVars->Neutron_Counter+aVars->Buffer_Block; }
      }

      while (!While_End)
      { /* we generate Coord[] and Coord_index[] from Buffer (auto limits) or passing neutron */
        if (aVars->Flag_Auto_Limits == 2)
        {
          if (While_Buffer < aVars->Buffer_Block)
          {
            /* first while loops (While_Buffer) */
            /* auto limits case : scan Buffer within limits and store in Mon2D */ 
            Coord[0] = aVars->Mon2D_Buffer[While_Buffer*(aVars->Coord_Number+2) + (aVars->Coord_Number)];
            
            for (i = 1; i <= aVars->Coord_Number; i++)
            {
              /* scanning variables in Buffer */
              XY = (aVars->Coord_Max[i]-aVars->Coord_Min[i]);
              Coord[i] = aVars->Mon2D_Buffer[While_Buffer*(aVars->Coord_Number+2) + (i-1)];
              if (XY > 0) Coord_Index[i] = floor((aVars->Mon2D_Buffer[(i-1) + While_Buffer*(aVars->Coord_Number+2)]-aVars->Coord_Min[i])*aVars->Coord_Bin[i]/XY);
              else Coord_Index[i] = 0;
              if (aVars->Flag_With_Borders)
              {
                if (Coord_Index[i] < 0) Coord_Index[i] = 0;
                if (Coord_Index[i] >= aVars->Coord_Bin[i]) Coord_Index[i] = aVars->Coord_Bin[i] - 1;
              }
            } /* end for */
            While_Buffer++;
          } /* end if in Buffer */
          else /* (While_Buffer >= aVars->Buffer_Block) && (aVars->Flag_Auto_Limits == 2) */ 
          {
            aVars->Flag_Auto_Limits = 0;
            if (!aVars->Flag_List) /* free Buffer not needed (no list to output) */
            { /* Dim : (aVars->Coord_Number+2)*aVars->Buffer_Block matrix (for p, dp) */ 
              free(aVars->Mon2D_Buffer);
            }
          }
        }
        else /* aVars->Flag_Auto_Limits == 0 or 1 */
        {
          for (i = 0; i <= aVars->Coord_Number; i++)
          { /* handle current neutron : last while */

            XY = 0;
            Set_aVars_Coord_Type = aVars->Coord_Type[i];
            if (Set_aVars_Coord_Type == aDEFS->COORD_X) XY = aVars->cx; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_Y) XY = aVars->cy; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_Z) XY = aVars->cz; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_VX) XY = aVars->cvx; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_VY) XY = aVars->cvy; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_VZ) XY = aVars->cvz; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_KX) XY = V2K*aVars->cvx; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_KY) XY = V2K*aVars->cvy; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_KZ) XY = V2K*aVars->cvz; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_SX) XY = aVars->csx; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_SY) XY = aVars->csy; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_SZ) XY = aVars->csz; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_T) XY = aVars->ct; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_P) XY = aVars->cp; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_HDIV) XY = RAD2DEG*atan2(aVars->cvx,aVars->cvz); 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_VDIV) XY = RAD2DEG*atan2(aVars->cvy,aVars->cvz); 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_V) XY = sqrt(aVars->cvx*aVars->cvx+aVars->cvy*aVars->cvy+aVars->cvz*aVars->cvz); 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_RADIUS) XY = sqrt(aVars->cx*aVars->cx+aVars->cy*aVars->cy); 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_K) { XY = sqrt(aVars->cvx*aVars->cvx+aVars->cvy*aVars->cvy+aVars->cvz*aVars->cvz);  XY *= V2K; }
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_ENERGY) { XY = aVars->cvx*aVars->cvx+aVars->cvy*aVars->cvy+aVars->cvz*aVars->cvz;  XY *= VS2E; }
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_LAMBDA) { XY = sqrt(aVars->cvx*aVars->cvx+aVars->cvy*aVars->cvy+aVars->cvz*aVars->cvz);  XY *= V2K; if (XY != 0) XY = 2*PI/XY; }
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_NCOUNT) XY = Coord[i]+1; 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_ANGLE) 
            {  XY = sqrt(aVars->cvx*aVars->cvx+aVars->cvy*aVars->cvy+aVars->cvz*aVars->cvz);
               if (aVars->cvz != 0) 
               {
                 XY= RAD2DEG*atan2(XY,aVars->cvz);
               } else XY = 0;
            }
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_THETA)  { if (aVars->cz != 0) XY = RAD2DEG*atan2(aVars->cx,aVars->cz); } 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_PHI) { if (aVars->cz != 0) XY = RAD2DEG*atan2(aVars->cy,aVars->cz); } 
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_USER1) XY = aVars->UserVariable1;
                  else
            if (Set_aVars_Coord_Type == aDEFS->COORD_USER2) XY = aVars->UserVariable2;
                  else
            XY = 0;

            Coord[i] = XY;
            if (!aVars->Flag_Auto_Limits)
            {
              XY = (aVars->Coord_Max[i]-aVars->Coord_Min[i]);
              if (XY > 0) Coord_Index[i] = floor((Coord[i]-aVars->Coord_Min[i])*aVars->Coord_Bin[i]/XY);
              else Coord_Index[i] = 0;
              if (aVars->Flag_With_Borders)
              {
                if (Coord_Index[i] < 0) Coord_Index[i] = 0;
                if (Coord_Index[i] >= aVars->Coord_Bin[i]) Coord_Index[i] = aVars->Coord_Bin[i] - 1;
              }
            } /* else Auto_Limits will get Index later from Buffer */
          } /* end for i */
          While_End = 1;
        } /* end else if aVars->Flag_Auto_Limits == 2 */

        if (aVars->Flag_Auto_Limits != 2) /* not when reading auto limits Buffer */
        { /* now store Coord into Buffer (no index needed) if necessary */
          if ((aVars->Buffer_Counter < aVars->Buffer_Block) && ((aVars->Flag_List) || (aVars->Flag_Auto_Limits == 1)))
          {
            for (i = 0; i < aVars->Coord_Number; i++)
            {
              aVars->Mon2D_Buffer[i + aVars->Neutron_Counter*(aVars->Coord_Number+2)] = Coord[i+1];
            }
            aVars->Mon2D_Buffer[aVars->Coord_Number + aVars->Neutron_Counter*(aVars->Coord_Number+2)]   = pp;
            aVars->Mon2D_Buffer[(aVars->Coord_Number+1) + aVars->Neutron_Counter*(aVars->Coord_Number+2)] = pp*pp;
            aVars->Buffer_Counter++;
            if (aVars->Flag_Verbose && (aVars->Buffer_Counter >= aVars->Buffer_Block) && (aVars->Flag_List == 1)) printf("Monitor_nD: %s %li neutrons stored in List.\n", aVars->compcurname, aVars->Buffer_Counter);
          }
          aVars->Neutron_Counter++;
        } /* end (aVars->Flag_Auto_Limits != 2) */

        /* store n1d/2d section for Buffer or current neutron in while */

        if (aVars->Flag_Auto_Limits != 1) /* not when storing auto limits Buffer */
        {
        /* 1D and n1D case : aVars->Flag_Multiple */
          if (aVars->Flag_Multiple)
          { /* Dim : aVars->Coord_Number*aVars->Coord_Bin[i] vectors (intensity is not included) */ 
            for (i= 0; i < aVars->Coord_Number; i++)
            {
              j = Coord_Index[i+1];
              if (j >= 0 && j < aVars->Coord_Bin[i+1])
              {
                aVars->Mon2D_N[i][j]++;
                aVars->Mon2D_p[i][j] += pp;
                aVars->Mon2D_p2[i][j] += pp*pp;
              }
            }
          }
          else /* 2D case : aVars->Coord_Number==2 and !aVars->Flag_Multiple and !aVars->Flag_List */
          if ((aVars->Coord_Number == 2) && !aVars->Flag_Multiple)
          { /* Dim : aVars->Coord_Bin[1]*aVars->Coord_Bin[2] matrix */
            i = Coord_Index[1];
            j = Coord_Index[2];
            if (i >= 0 && i < aVars->Coord_Bin[1] && j >= 0 && j < aVars->Coord_Bin[2])
            {
              aVars->Mon2D_N[i][j]++;
              aVars->Mon2D_p[i][j] += pp;
              aVars->Mon2D_p2[i][j] += pp*pp;
            }
          }
        } /* end (aVars->Flag_Auto_Limits != 1) */
      } /* end while */
    } /* end if aVars->Coord_Number */
  }
  
#endif
  
  MonitornD_Defines_type DEFS;
  MonitornD_Variables_type Vars;

%}

INITIALIZE
%{

strcpy(Vars.compcurname, mccompcurname);
strcpy(Vars.option, options);

Monitor_nD_Init(&DEFS, &Vars, xwidth, yheigth, zthick, xmin,xmax,ymin,ymax,zmin,zmax);

%}
TRACE
%{    
  double  XY=0;
  double  t0 = 0;
  double  t1 = 0;
  long    i,j;
  double  pp;
  double  Coord[MONnD_COORD_NMAX];
  long    Coord_Index[MONnD_COORD_NMAX];
  char    While_End   =0;
  long    While_Buffer=0;
  int     intersect = 0;
  char    Set_Vars_Coord_Type = DEFS.COORD_NONE;
  double  cp  = p ;
  double  cx  = x ;
  double  cvx = vx;
  double  csx = sx;
  double  cy  = y ;
  double  cvy = vy;
  double  csy = sy;
  double  cz  = z ;
  double  cvz = vz;
  double  csz = sz;
  double  ct  = t ;    

  if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) /* square xy */
    intersect = (x>=Vars.mxmin && x<=Vars.mxmax && y>=Vars.mymin && y<=Vars.mymax);
  else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK)   /* disk xy */
    intersect = ((x*x + y*y) <= Vars.Sphere_Radius*Vars.Sphere_Radius);
  else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) /* sphere */
  {
    intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius);
  /*      intersect = (intersect && t0 > 0); */
  }
  else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) /* cylinder */
  {
    intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, Vars.Sphere_Radius, Vars.Cylinder_Height);
  }
  else if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX) /* box */
  {
    intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, fabs(Vars.mxmax-Vars.mxmin), fabs(Vars.mymax-Vars.mymin), fabs(Vars.mzmax-Vars.mzmin));
  }

  if (intersect)
  {
    if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX))
    {
      if (t0 < 0 && t1 > 0)
        t0 = t;  /* neutron was already inside ! */
      if (t1 < 0 && t0 > 0) /* neutron exit before entering !! */
        t1 = t;
      /* t0 is now time of incoming intersection with the sphere. */
      if ((Vars.Flag_Shape < 0) && (t1 > 0))
        PROP_DT(t1); /* t1 outgoing beam */
      else
        PROP_DT(t0); /* t0 incoming beam */
    }
    else
      PROP_Z0; 

    /* Now get the data to monitor: curent or keep from PreMonitor */
    if (Vars.Flag_UsePreMonitor != 1)
    {
      Vars.cp = p;
      Vars.cx = x;
      Vars.cvx = vx;
      Vars.csx = sx;
      Vars.cy = y;
      Vars.cvy = vy;
      Vars.csy = sy;
      Vars.cz = z;
      Vars.cvz = vz;
      Vars.csz = sz;
      Vars.ct = t;
    }

    if ((Vars.He3_pressure > 0) && (t1 != t0) && ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX)))
    {
      XY = exp(-7.417*Vars.He3_pressure*fabs(t1-t0)*2*PI/V2K);
      /* will monitor the absorbed part */
      Vars.cp *= 1-XY;
      /* and modify the neutron weight after monitor, only remains 1-p_detect */
      p *= XY;
    }

    if (Vars.Flag_per_cm2 && Vars.area != 0) Vars.cp /= Vars.area;
    Vars.Nsum++;
    Vars.psum += Vars.cp;
    Vars.p2sum += Vars.cp*Vars.cp;

    Monitor_nD_Trace(&DEFS, &Vars);
    
    /* now handles intermediate results saving */
    if ((Vars.Intermediate > 0) && (mcget_run_num() > Vars.IntermediateCnts))
    {
      Vars.IntermediateCnts += Vars.Intermediate*mcget_ncount();
      /* save results, but do not free pointers */
      Monitor_nD_OutPut(&DEFS, &Vars, 0);
    }
    if (Vars.Flag_parallel) /* back to neutron state before detection */
    {
      p  = cp ;
      x  = cx ;
      vx = cvx;
      sx = csx;
      y  = cy ;
      vy = cvy;
      sy = csy;
      z  = cz ;
      vz = cvz;
      sz = csz;
      t  = ct ;
    }
    
  } /* end if intersection */
  else if (Vars.Flag_Absorb) ABSORB;
  
 %}
 
FINALLY
 %{
    /* save results, and free pointers */
    Monitor_nD_OutPut(&DEFS, &Vars, 1);
%}

MCDISPLAY
%{
  double radius, h;
  double xmin;
  double xmax;
  double ymin;
  double ymax;
  double zmin;
  double zmax;
  
  radius = Vars.Sphere_Radius;
  h = Vars.Cylinder_Height;
  xmin = Vars.mxmin;
  xmax = Vars.mxmax;
  ymin = Vars.mymin;
  ymax = Vars.mymax;
  zmin = Vars.mzmin;
  zmax = Vars.mzmax;
  
  if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE)
  {
    magnify("");
    circle("xy",0,0,0,radius);
    circle("xz",0,0,0,radius);
    circle("yz",0,0,0,radius);
  }
  else
  if (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK)
  {
    magnify("");
    circle("xy",0,0,0,radius);
  }
  else
  if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE)
  {
    magnify("xy");
    multiline(5, (double)xmin, (double)ymin, 0.0,
           (double)xmax, (double)ymin, 0.0,
           (double)xmax, (double)ymax, 0.0,
           (double)xmin, (double)ymax, 0.0,
           (double)xmin, (double)ymin, 0.0);
  }
  else
  if (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND)
  {
    magnify("xyz");
    circle("xz", 0,  h/2.0, 0, radius);
    circle("xz", 0, -h/2.0, 0, radius);
    line(-radius, -h/2.0, 0, -radius, +h/2.0, 0);
    line(+radius, -h/2.0, 0, +radius, +h/2.0, 0);
    line(0, -h/2.0, -radius, 0, +h/2.0, -radius);
    line(0, -h/2.0, +radius, 0, +h/2.0, +radius);
  }
  else
  if (abs(Vars.Flag_Shape) == DEFS.SHAPE_BOX)
  {
    magnify("xyz");
    multiline(5, xmin, ymin, zmin,
                 xmax, ymin, zmin,
                 xmax, ymax, zmin,
                 xmin, ymax, zmin,
                 xmin, ymin, zmin);
    multiline(5, xmin, ymin, zmax,
                 xmax, ymin, zmax,
                 xmax, ymax, zmax,
                 xmin, ymax, zmax,
                 xmin, ymin, zmax);
    line(xmin, ymin, zmin, xmin, ymin, zmax);
    line(xmax, ymin, zmin, xmax, ymin, zmax);
    line(xmin, ymax, zmin, xmin, ymax, zmax);
    line(xmax, ymax, zmin, xmax, ymax, zmax);
  }
%}

END

From farhi at ill.fr  Mon Aug  6 12:09:00 2001
From: farhi at ill.fr (Emmanuel Farhi)
Date: Mon, 06 Aug 2001 12:09:00 +0200
Subject: Gravity_guide and Monitor_nD components for McStas 1.4
Message-ID: <3B6E6CBC.73F7DA08@ill.fr>

An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20010806/2535a50f/attachment.html>
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Gravity_guide.comp
*         Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: <a href="mailto:farhi at ill.fr">Emmanuel Farhi</a>
* Date: Aug 03 2001
* Version: $Revision: 1.4 $
* Origin: <a href="http://www.ill.fr">ILL (France)</a>. Aug 03 2001.
* Modified by: E. Farhi, from Gravity_guide by K. Lefmann (buggy).
*
* Neutron guide with gravity.
*
* %D
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane. Gravitation applies also when reaching the guide input
* window. The guide can be channeled (k,d parameters). The guide coating
* specifications may be entered via different ways (global, opposite side 
* pars, each wall m-value).
* For details on the geometry calculation see the description in the McStas
* reference manual.
*
* %P
* INPUT PARAMETERS:
*
* w1:      (m)    Width at the guide entry
* h1:      (m)    Height at the guide entry
* w2:      (m)    Width at the guide exit
* h2:      (m)    Height at the guide exit
* l:       (m)    length of guide
* R0:      (1)    Low-angle reflectivity
* Qc:      (AA-1) Critical scattering vector
* alpha:   (AA)   Slope of reflectivity
* m:       (1)    m-value of material. Zero means completely absorbing.
* W:       (AA-1) Width of supermirror cut-off
* d:       (m)    Thickness of subdividing walls [0]
* k:       (1)    Number of channels in the guide (>= 1) [1]
*
* Optional input parameters: (different ways for m-specifications)
*
* G:       (m/s2) Gravitation acceleration along y axis [-9.81]
* Gx:      (m/s2) Gravitation acceleration along x axis [0]
* Gy:      (m/s2) Gravitation acceleration along y axis [-9.81]
* Gz:      (m/s2) Gravitation acceleration along z axis [0]
* mh:      (1)    m-value of material for left/right vert. mirrors
* mv:      (1)    m-value of material for top/bottom horz. mirrors
* mx:      (1)    m-value of material for left/right vert. mirrors
* my:      (1)    m-value of material for top/bottom horz. mirrors
* mleft:   (1)    m-value of material for left.   vert. mirror
* mright:  (1)    m-value of material for right.  vert. mirror
* mtop:    (1)    m-value of material for top.    horz. mirror
* mbottom: (1)    m-value of material for bottom. horz. mirror
* 
*
* %D
* Example values: m=4 Qc=0.02 W=1/300 alpha=6.49 R0=1
*
* %E
*******************************************************************************/

DEFINE COMPONENT Gravity_guide
DEFINITION PARAMETERS ()
SETTING PARAMETERS (w1, h1, w2, h2, l, 
  R0=0.99, Qc=0.021, alpha=6.07, m=1.0, W=0.003, k=1, d=0.0005,
  Gx=0,Gy=-9.81,Gz=0, G=0,
  mh=-1,mv=-1,mx=-1, my=-1, 
  mleft=-1, mright=-1, mtop=-1, mbottom=-1)
OUTPUT PARAMETERS (gx,gy,gz,nx,ny,nz,wx,wy,wz,A,norm_n,norm_n2,N_reflection,w1c,w2c,M)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)

DECLARE 
%{
#ifndef Gravity_guide_here
#define Gravity_guide_here
 /* this function calculates the intersection between a neutron trajectory
  * and a plane with acceleration gx,gy,gz. The neutron starts at point x,y,z
  * with velocity vx, vy, vz. The plane has a normal vector nx,ny,nz and 
  * contains the point wx,wy,wz
  * The function returns 0 if no intersection occured after the neutron started
  * and 1 if there is an intersection. Then *Idt is the elapsed time until 
  * the neutron hits the roof.
  */
  /* Let n=(nx,ny,nz) be the normal plane vector (one of the six sides) 
   * Let W=(wx,wy,wz) be Any point on this plane (for instance at z=0)
   * The problem consists in solving the 2nd order equation:
   *      1/2.n.g.t^2 + n.v.t + n.(r-W) = 0 (1)
   * Without acceleration, t=-n.(r-W)/n.v
   */
  
  int plane_intersect_Gfast(double *Idt, 
                  double A,  double B,  double C)
  {
    /* plane_intersect_Gfast(&dt, A, B, C)
     * A = 0.5 n.g; B = n.v; C = n.(r-W);
     * no cceleration when A=0
     */
    double D, sD;
    double dt1, dt2;
    
    *Idt = -1;
    
    if (A == 0) /* this plane is parallel to the acceleration */
    {
      if (B == 0)  /* the speed is parallel to the plane, no intersection */
        return (0);
      else  /* no acceleration case */
        { *Idt = -C/B; 
          if (*Idt >= 0) return (2);
          else return (0); }
    }
    else
    {
      /* Delta > 0: neutron trajectory hits the mirror */
      D = B*B - 4*A*C;
      if (D >= 0)
      {
        sD = sqrt(D);
        dt1 = (-B + sD)/2/A;
        dt2 = (-B - sD)/2/A;
        if (dt1 <0 && dt2 >=0) *Idt = dt2;
        else
        if (dt2 <0 && dt1 >=0) *Idt = dt1;
        else
        if (dt1 <0 && dt2 < 0) return (0);
        else
        if (dt1 < dt2) *Idt = dt1;
        else
          *Idt = dt2;
        return (1);
      }
      else  /* Delta <0: no intersection */
        return (0);
    }     
  }
  
  #define PROP_GRAV_DT(dt, Ax, Ay, Az) \
  do { \
    mcnlx  += mcnlvx*dt + Ax*dt*dt/2; \
    mcnly  += mcnlvy*dt + Ay*dt*dt/2; \
    mcnlz  += mcnlvz*dt + Az*dt*dt/2; \
    mcnlvx += Ax*dt; \
    mcnlvy += Ay*dt; \
    mcnlvz += Az*dt; \
    mcnlt  += dt; \
  } while(0)
  
  /* The rotation of axes is a variable named "mcrota" ## mccompcurname */
  /* The position of comp is a variable named "mcposa" ## mccompcurname */
  
#endif
    
  double gx;
  double gy;
  double gz;
  double nx[6], ny[6], nz[6];
  double wx[6], wy[6], wz[6];
  double A[6], norm_n2[6], norm_n[6];
  long   N_reflection[7]={0,0,0,0,0,0,0};
  double w1c;
  double w2c;
  double M[5]={0,0,0,0,0};
%}

INITIALIZE
%{
  int i;
  
  gx = Gx; /* The gravitation vector in the current component axis system */
  if (G) gy = G; else gy = Gy;
  gz = Gz;
  
  if (k == 0) { fprintf(stderr,"%s: Fatal: no channel in this guide (k=0).\n", mccompcurname); exit(-1); }
  if (d < 0) { fprintf(stderr,"%s: Fatal: subdividing walls have negative thickness in this guide (d<0).\n", mccompcurname); exit(-1); }
  w1c = (w1 + d)/(double)k;
  w2c = (w2 + d)/(double)k;
  
  for (i=0; i <= 4; M[i++]=m);
  if (mx >= 0) { M[1] = mx; M[2] = mx; }
  if (mh >= 0) { M[1] = mh; M[2] = mh; }
  if (my >= 0) { M[1] = my; M[2] = my; }
  if (mv >= 0) { M[1] = mv; M[2] = mv; }
  if (mleft   >= 0) M[1] = mleft  ;
  if (mright  >= 0) M[2] = mright ;
  if (mtop    >= 0) M[3] = mtop   ;
  if (mbottom >= 0) M[4] = mbottom;
  
  /* This is now the downward gravitation vector */
                                                     
  nx[1] =  l; ny[1] =  0; nz[1] = -0.5*(w2c-w1c); /* 1:+X left       */
  nx[2] = -l; ny[2] =  0; nz[2] =  nz[1];         /* 2:-X right      */
  nx[3] =  0; ny[3] =  l; nz[3] = -0.5*(h2-h1);   /* 3:+Y top        */
  nx[4] =  0; ny[4] = -l; nz[4] =  nz[3];         /* 4:-Y bottom     */
  nx[5] =  0; ny[5] =  0; nz[5] =  1;             /* 5:+Z exit       */
  nx[0] =  0; ny[0] =  0; nz[0] = -1;             /* 0:Z0 input      */
  
  wx[1] = +(w1c-d)/2; wy[1] =  0;    wz[1] = 0;   /* 1:+X left       */
  wx[2] = -(w1c-d)/2; wy[2] =  0;    wz[2] = 0;   /* 2:-X right      */
  wx[3] =  0;         wy[3] = +h1/2; wz[3] = 0;   /* 3:+Y top        */
  wx[4] =  0;         wy[4] = -h1/2; wz[4] = 0;   /* 4:-Y bottom     */
  wx[5] =  0;         wy[5] =  0;    wz[5] = l;   /* 5:+Z exit       */
  wx[0] =  0;         wy[0] =  0;    wz[0] = 0;   /* 0:Z0 input      */
  
  for (i=0; i <= 5; i++)
  {
    A[i] = scalar_prod(nx[i], ny[i], nz[i], gx, gy, gz)/2;
    norm_n2[i] = nx[i]*nx[i] + ny[i]*ny[i] + nz[i]*nz[i];
    if (norm_n2[i] <= 0)
      { fprintf(stderr,"%s: Fatal: normal vector norm %i is null/negative ! Check guide dimensions.\n", mccompcurname, i); exit(-1); } /* should never occur */
    else
      norm_n[i] = sqrt(norm_n2[i]);
  }
  
  
%}

TRACE
%{
  double n_dot_v[6];
  double B, C, dt0, dt;
  double q;
  int    ret, side;
  double edge;
  double hadj;                                  /* Channel displacement */

  
  dt = -1; dt0 = -1;
  /* propagate to box input (with gravitation) in comp local coords */
  /* 0=Z0 side: n=(0, 0, 1) ; W = (0, 0, 0) (at z=0, guide input)*/
  B = -vz; C = -z;
  ret = plane_intersect_Gfast(&dt0, A[0], B, C);
  if (ret && dt0>0)
  { 
    dt = dt0; 
    PROP_GRAV_DT(dt, gx, gy, gz);
    N_reflection[6]++;
  }
  /* check if we are in the box input, else absorb */
  if(dt < 0 || x <= -w1/2 || x >= w1/2 || y <= -h1/2 || y >= h2/2)
    { ABSORB; }
    
  /* Shift origin to center of channel hit (absorb if hit dividing walls) */
  x += w1/2.0;
  edge = floor(x/w1c)*w1c;
  if(x - edge > w1c - d)
  {
    x -= w1/2.0; /* Re-adjust origin */
    ABSORB;
  }
  x -= (edge + (w1c - d)/2.0);
  hadj = edge + (w1c - d)/2.0 - w1/2.0;

  /* neutron is now in the input window of the guide */
  
  /* do loops on reflections in the box */
  for(;;)
  {
    /* get intersections for all box sides */
    /* A = 0.5 n.g; B = n.v; C = n.(r-W); 
    
       A = scalar_prod(nx,ny,nz,gx,gy,gz)/2;
       B = scalar_prod(nx,ny,nz,vx,vy,vz);
       C = scalar_prod(nx,ny,nz,x-wx,y-wy,z-wz);
    */
    side = 0;
    /* starts with the exit side intersection (the last one !)*/
    
   /* 5=+Z side: n=(0, 0, 1) ; W = (0, 0, l) (at z=l, guide exit)*/
    B = vz; C = z - wz[5];
    ret = plane_intersect_Gfast(&dt0, A[5], B, C);
    if (ret && dt0>0)     
    { dt = dt0; side=5;  
      n_dot_v[5] = B; }
    else
    { fprintf(stderr,"%s: warning: neutron trajectory is parallel to guide exit, and thus can not exit\n", mccompcurname); x += hadj; ABSORB; }
      
    /* now look if there is a previous intersection with guide sides */
    
    /* 1=+X side: n=(l, 0, -0.5*(w2-w1)) ; W = (+w1/2, 0, 0) (left)*/ 
    B = nx[1]*vx + nz[1]*vz; C = nx[1]*(x-wx[1]) + nz[1]*z; /* ny=wz=0 */
    ret = plane_intersect_Gfast(&dt0, A[1], B, C);
    if (ret && dt0>10e-10 && dt0<dt)
    { dt = dt0; side=1; n_dot_v[1] = B; }
    
    /* 2=-X side: n=(l, 0, +0.5*(w2-w1)) ; W = (-w1/2, 0, 0) (right) */
    B = nx[2]*vx + nz[2]*vz; C = nx[2]*(x-wx[2]) + nz[2]*z; /* ny=wz=0 */
    ret = plane_intersect_Gfast(&dt0, A[2], B, C);
    if (ret && dt0>10e-10 && dt0<dt)
    { dt = dt0; side=2; n_dot_v[2] = B; }
    
    /* 3=+Y side: n=(0, l, -0.5*(h2-h1)) ; W = (0, +h1/2, 0) (up) */
    B = ny[3]*vy + nz[3]*vz; C = ny[3]*(y-wy[3]) + nz[3]*z; /* nx=wz=0 */
    ret = plane_intersect_Gfast(&dt0, A[3], B, C);
    if (ret && dt0>10e-10 && dt0<dt)
    { dt = dt0; side=3; n_dot_v[3] = B; }
    
    /* 4=-Y side: n=(0, l, +0.5*(h2-h1)) ; W = (0, -h1/2, 0) (down) */
    B = ny[4]*vy + nz[4]*vz; C = ny[4]*(y-wy[4]) + nz[4]*z; /* nx=wz=0 */
    ret = plane_intersect_Gfast(&dt0, A[4], B, C);
    if (ret && dt0>10e-10 && dt0<dt)
    { dt = dt0; side=4; n_dot_v[4] = B; }
    
    /* only positive dt are valid */
    /* exit reflection loops if no intersection (neutron is after box) */
    if (side == 0 || dt < 0)
      { fprintf(stderr,"%s: warning: neutron has entered guide, but can not exit !\n", mccompcurname); ABSORB; } /* should never occur */   
     /* 
    if (side < 5 && (x < -w1 || x > w1 || y < -h1 || y > h2 ))  
      ABSORB;  */ /* neutron has left guide through wall */
      
    /* propagate to dt */
    PROP_GRAV_DT(dt, gx, gy, gz);
    
    /* do reflection on speed for l/r/u/d sides */
    if (side == 5) /* neutron reaches end of guide: end loop and exit comp */
      { N_reflection[side]++; break; }
    /* else reflection on a guide wall */
    if(M[side] == 0 || Qc == 0)  /* walls are absorbing */
      { x += hadj; ABSORB; }
    /* change/mirror velocity: v_f = v - n.2*n.v/|n|^2 */  
    N_reflection[side]++; /* norm_n2 > 0 was checked at INIT */
    dt0 = 2*n_dot_v[side]/norm_n2[side]; /* 2*n.v/|n|^2 */
    vx -= nx[side]*dt0;
    vy -= ny[side]*dt0;
    vz -= nz[side]*dt0;
    
    /* compute q and modify neutron weight */
    /* scattering q=|k_i-k_f| = V2Q*|vf - v| = V2Q*2*n.v/|n| */
    q = 2*V2Q*fabs(n_dot_v[side])/norm_n[side];
    if(q > Qc)
    {
      double arg;
      if (W>0)
        arg = (q-M[side]*Qc)/W;
      else
        arg = (q-M[side]*Qc)*10000; /* W = 0.00001 */

      if(arg < 10)
        p *= .5*(1-tanh(arg))*(1-alpha*(q-Qc));
      else
        { x += hadj; ABSORB; };                               /* Cutoff ~ 1E-10 */
     }
     p *= R0;
     x += hadj; SCATTER; x -= hadj;
     N_reflection[0]++;
     /* go to the next reflection */
   }
    x += hadj; /* Re-adjust origin after SCATTER */
%}

MCDISPLAY
%{
  double x;
  int i;

  magnify("xy");
  for(i = 0; i < k; i++)
  {
    multiline(5,
              i*w1c - w1/2.0, -h1/2.0, 0.0,
              i*w2c - w2/2.0, -h2/2.0, (double)l,
              i*w2c - w2/2.0,  h2/2.0, (double)l,
              i*w1c - w1/2.0,  h1/2.0, 0.0,
              i*w1c - w1/2.0, -h1/2.0, 0.0);
    multiline(5,
              (i+1)*w1c - d - w1/2.0, -h1/2.0, 0.0,
              (i+1)*w2c - d - w2/2.0, -h2/2.0, (double)l,
              (i+1)*w2c - d - w2/2.0,  h2/2.0, (double)l,
              (i+1)*w1c - d - w1/2.0,  h1/2.0, 0.0,
              (i+1)*w1c - d - w1/2.0, -h1/2.0, 0.0);
  }
  line(-w1/2.0, -h1/2.0, 0.0, w1/2.0, -h1/2.0, 0.0);
  line(-w2/2.0, -h2/2.0, (double)l, w2/2.0, -h2/2.0, (double)l);
%}

END

From per-olof.aastrand at risoe.dk  Wed Aug  8 18:20:26 2001
From: per-olof.aastrand at risoe.dk (Per-Olof Astrand)
Date: Wed, 08 Aug 2001 18:20:26 +0200 (CEST)
Subject: McStas: next version
Message-ID: <Pine.LNX.3.96.1010808171759.25095A-100000@fys-lin-1.risoe.dk>


Dear McStas user,

Version 1.5 of McStas including an updated manual is planned to be
released around the 10th of September. We think it would be very fruitful
for all of us if minor problems or extensions in the source code, official
components or the manual are reported to us so that they can be included
in the next official releas.

We will also separate the McStas manual into two documents: the "McStas
manual" (excluding the component descriptions) and a new "McStas component
manual" including component descriptions only. The latter document will
only be edited by the McStas administrators, and the actual author of the
component will hopefully also be the author of the contribution in
the McStas component manual. In this way, we hope that the authorship of a
McStas component is recognized appropriately. We therefore kindly ask the
McStas user community to send us new components to be included in the
next release. 

Hopefully, we can all benefit from a more extensive exchange of 
components. I have discussed this with some people and a regular user tend
to think that his/her components are simple and can be written by anyone.  
My main argument for submitting a component anyway is that it has already
been tested and validated extensively by the user and it is therefore a
very good chance for that it is correct and a lot of duplicate work is
thus saved.

Kind regards,

Per-Olof Astrand
mcstas at risoe.dk






From per-olof.aastrand at risoe.dk  Fri Aug 10 12:07:37 2001
From: per-olof.aastrand at risoe.dk (Per-Olof Astrand)
Date: Fri, 10 Aug 2001 12:07:37 +0200 (CEST)
Subject: McStas: representation of polarisation
Message-ID: <Pine.LNX.3.96.1010810114303.12915B-100000@fys-lin-1.risoe.dk>


Dear McStas user,

After some discussions with Kristian Nielsen, we have finally sorted out
how polarisation is supposed to work in McStas. Only components actually
modifying the spin vector explicitly (modifying the (sx,sy,sz) vector)
requires the POLARISATION PARAMETERS(sx,sy,sz) line. If the instrument
contains at least one component with this line it will handle polarization
correctly, i.e. the spin vector is correctly transformed to the local
coordinate system of the components and components containing this line
are allowed to modify (sx,sy,sz).

It is also of importance that non-polarising components do not contain
this line because a simulation of an instrument without polarisation is
more efficient if the transformation of the spin is excluded. This will be
cleaned up in the components of McStas 1.5 (but it already works in the
McStas 1.4 kernel) and it will also be properly documented in the manual.

You can follow the development of McStas 1.5 at
http://neutron.risoe.dk/mcstas/developments/developments.html . We have a
quite long list of things to fix, but further comments on problems and
extensions are welcome.

Best regards,

Per-Olof 


Per-Olof ?strand
Dept. of Chemistry, University of Copenhagen 
and Materials Research Department, Ris? National Laboratory 
per-olof.aastrand at theory.ki.ku.dk or per-olof.aastrand at risoe.dk
http://theochem.ki.ku.dk/~peo/ 




From kim.lefmann at risoe.dk  Tue Sep  4 14:59:07 2001
From: kim.lefmann at risoe.dk (Kim Lefmann)
Date: Tue, 04 Sep 2001 14:59:07 +0200 (CEST)
Subject: New monitor components
Message-ID: <Pine.LNX.3.96.1010904145208.23372I-100000@fys-lin-1.risoe.dk>

Dear McStas users,

I have cleaned up in my component directories before the 1.5 release. This
resulted in 5 "new" monitor components, that I have found useful during
the last two years. Note that there is a big overlap with Emmanuels
monitor_nD, which covers all but the last of my monitors.

PSDlin:     Rectangular, 1D, I vs. position
EPSD:       Rectangular, 2D, I vs. position and energy
PSDcyl:     Cylindrical, 2D, I vs. position (angle and height) 
TOF_cylPSD: Cylindrical, 2D, I vs. time and position (angle)
TOFlog:     Rectangular, 1D, I vs. time (logarithmically binned)

Yours,
          Kim

------------------------------
Kim Lefmann
Senior Scientist
Dept. Cond. Matt. Phys. & Chem.
Risoe National Laboratory

Phone: +45 46 77 47 26
Fax:   +45 46 77 47 90



From kim.lefmann at risoe.dk  Fri Sep  7 16:23:18 2001
From: kim.lefmann at risoe.dk (Kim Lefmann)
Date: Fri, 07 Sep 2001 16:23:18 +0200 (CEST)
Subject: New component: Source_Maxwell_3
Message-ID: <Pine.LNX.3.96.1010907162038.26667A-100000@fys-lin-1.risoe.dk>

Dear McStas users,

We have added a new component, which is a source described by the sum of 3
Maxwellian distributions. Source_Maxwell_3 . Flux parameters for the SINQ
cold source is given as example. 

Cheers,
          Kim

------------------------------
Kim Lefmann
Senior Scientist
Dept. Cond. Matt. Phys. & Chem.
Risoe National Laboratory

Phone: +45 46 77 47 26
Fax:   +45 46 77 47 90



From per-olof.aastrand at risoe.dk  Tue Sep 18 09:10:23 2001
From: per-olof.aastrand at risoe.dk (Per-Olof =?iso-8859-1?Q?=C5strand?=)
Date: Tue, 18 Sep 2001 09:10:23 +0200
Subject: McStas 1.5 delayed a couple of days
Message-ID: <3BA6F35F.6A274E74@risoe.dk>

Dear McStas user,

For those of you attending the ICNS meeting in Munich last week, you
perhaps remember the interesting comparison of different simulation
packages presented by Phil Seeger. Emmanuel Farhi has investigated the
differences for McStas and found an error in how the illumination is
calculated in the sources. This will be corrected in McStas version 1.5
and the release will therefore be delayed a couple of days.

Best regards,

Per-Olof ?strand

--
Per-Olof ?strand
Dept. of Chemistry, University of Copenhagen
and Materials Research Department, Ris? National Laboratory
per-olof.aastrand at theory.ki.ku.dk or per-olof.aastrand at risoe.dk
http://theochem.ki.ku.dk/~peo/





From per-olof.aastrand at risoe.dk  Tue Sep 18 12:24:47 2001
From: per-olof.aastrand at risoe.dk (Per-Olof =?iso-8859-1?Q?=C5strand?=)
Date: Tue, 18 Sep 2001 12:24:47 +0200
Subject: McStas: components on web-page
Message-ID: <3BA720EF.73AE7836@risoe.dk>

Dear McStas user,

Today, there is not a one-to-one correspondence between the components
in the official McStas release and the components generated by McDoc on
our web page. The web page also includes a number of unofficial
components and is identical to the ILL page. This has caused some
confusion and is changed when version 1.5 is released.

There will be a "one-to-one-to-one" correspondence between the
components included in the McStas release (the tar file that is
downloaded), the McStas web-page, and the forthcoming McStas component
manual, which means that the McStas web-page will only include official
components and components that will become official in the next release.
It is of course still possible for others (like ILL) to create a
web-page using McDoc which includes additional unofficial components .

It is still possible to submit components to be included in McStas 1.5,
and if you have a component that may be useful for the McStas community,
you should consider submitting it.  The first McStas component manual
will be released in late October/early November, and I will (soon)
return to component authors with details on how to include a
contribution to the manual.

Best regards,

Per-Olof ?strand

--
Per-Olof ?strand
Dept. of Chemistry, University of Copenhagen
and Materials Research Department, Ris? National Laboratory
per-olof.aastrand at theory.ki.ku.dk or per-olof.aastrand at risoe.dk
http://theochem.ki.ku.dk/~peo/





From jens.klenke at hmi.de  Thu Sep 27 14:08:36 2001
From: jens.klenke at hmi.de (Jens Klenke)
Date: Thu, 27 Sep 2001 14:08:36 +0200
Subject: PSD_monitor in McStas
Message-ID: <01092714204400.05298@nitrit>

Dear McStas user,

how can I get all three 2dim arrays (N(x,y),p(x,y),p2(x,y))
from the component psd_monitor. I only get a 2 dim array
with N(x,y). Kernel call DETEKTOR_OUT_2D should print out
three all 2dim arrays?
I use PSD_monitor revision 1.7, McStas V1.4 
on a Linux System.

Best regards
 
Jens
 -- 
Jens Klenke, Abt. SF2 Magnetismus
Hahn-Meitner-Institut
Glienickerstra?e 100, 14109 Berlin
Tel ++49-30-8062-3167
FAX: ++49-30-8062-2999