Still the 2e9...
Emmanuel Farhi
farhi at ill.fr
Mon Mar 19 20:28:53 CET 2001
Hy,
I just realized that the '2e9' bug is still there, but more vicious.
in fact, I passed ncount and run_num to long, because there may be
rounding errors when using float numbers, and comparisons. I got
unending simulations...
Please use the attached files ! (lib/mcstas)
Changing that requires to update Monitor_nD... (lib/mcstas/monitors)
Cheers.
--
What's up Doc ?
--------------------------------------------
Emmanuel FARHI, http://www.ill.fr/tas/people/Farhi.html \|/ ____ \|/
CS-Group ILL4/156, Institut Laue-Langevin (ILL) Grenoble ~@-/ oO \-@~
6 rue J. Horowitz, BP 156, 38042 Grenoble Cedex 9,France /_( \__/ )_\
Work :Tel (33/0) 4 76 20 71 35. Fax (33/0) 4 76 48 39 06 \__U_/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20010319/b5d7981f/attachment.html>
-------------- next part --------------
/*******************************************************************************
* Runtime system for McStas.
*
* Project: Monte Carlo Simulation of Triple Axis Spectrometers
* File name: mcstas-r.h
*
* Author: K.N. Aug 29, 1997
*
* Copyright (C) Risoe National Laboratory, 1997-1998, All rights reserved
*******************************************************************************/
#ifndef MCSTAS_R_H
#define MCSTAS_R_H
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
/* If the runtime is embedded in the simulation program, some definitions can
be made static. */
#ifdef MC_EMBEDDED_RUNTIME
#define mcstatic static
#else
#define mcstatic
#endif
#ifdef __dest_os
#if (__dest_os == __mac_os)
#define MAC
#endif
#endif
#ifdef WIN32
#define MC_PATHSEP_C '\\'
#define MC_PATHSEP_S "\\"
#else /* !WIN32 */
#ifdef MAC
#define MC_PATHSEP_C ':'
#define MC_PATHSEP_S ":"
#else /* !MAC */
#define MC_PATHSEP_C '/'
#define MC_PATHSEP_S "/"
#endif /* !MAC */
#endif /* !WIN32 */
#ifndef MAC
#ifndef WIN32
#include <signal.h>
#endif /* !MAC */
#endif /* !WIN32 */
typedef double MCNUM;
typedef struct {MCNUM x, y, z;} Coords;
typedef MCNUM Rotation[3][3];
/* Note: the enum instr_formal_types definition MUST be kept
synchronized with the one in mcstas.h and with the
instr_formal_type_names array in cogen.c. */
enum instr_formal_types
{
instr_type_double, instr_type_int, instr_type_string
};
struct mcinputtable_struct {
char *name;
void *par;
enum instr_formal_types type;
};
extern struct mcinputtable_struct mcinputtable[];
extern int mcnumipar;
extern char mcinstrument_name[], mcinstrument_source[];
extern int mctraceenabled, mcdefaultmain;
void mcinit(void);
void mcraytrace(void);
void mcfinally(void);
void mcdisplay(void);
/* Adaptive search tree definitions. */
typedef double adapt_t;
/*******************************************************************************
* Structure of an adaptive search tree. The v array runs from 0 to N-1 (a
* total of N elements) and holds the values of each node. The sum of all v
* values is in total.
* The s array runs from 0 to N and is used to represents the cumulative sum
* of v[0] through v[i-1]. The number represented is the sum of s[i] and all
* its parents up to the root node.
*******************************************************************************/
struct adapt_tree
{
adapt_t *s, *v, total;
int N; /* < 1 << (depth+1) */
int depth;
int root; /* = (1 << depth) - 1 */
int initstep; /* = 1 << (depth-1) */
};
int adapt_tree_search(struct adapt_tree *t, adapt_t v);
void adapt_tree_add(struct adapt_tree *t, int i, adapt_t v);
struct adapt_tree * adapt_tree_init(int N);
void adapt_tree_free(struct adapt_tree *t);
#define SCATTER do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
mcnlt,mcnlsx,mcnlsy, mcnlp);} while(0)
#define ABSORB do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
mcnlt,mcnlsx,mcnlsy, mcnlp); mcDEBUG_ABSORB(); goto mcabsorb;} while(0)
/* Note: The two-stage approach to MC_GETPAR is NOT redundant; without it,
* after #define C sample, MC_GETPAR(C,x) would refer to component C, not to
* component sample. Such are the joys of ANSI C.
*/
#define MC_GETPAR2(comp, par) (mcc ## comp ## _ ## par)
#define MC_GETPAR(comp, par) MC_GETPAR2(comp,par)
#define DETECTOR_OUT(p0,p1,p2) mcdetector_out(mccompcurname,p0,p1,p2,NULL)
#define DETECTOR_OUT_0D(t,p0,p1,p2) mcdetector_out_0D(t,p0,p1,p2,mccompcurname)
#define DETECTOR_OUT_1D(t,xl,yl,xvar,x1,x2,n,p0,p1,p2,f) \
mcdetector_out_1D(t,xl,yl,xvar,x1,x2,n,p0,p1,p2,f,mccompcurname)
#define DETECTOR_OUT_2D(t,xl,yl,x1,x2,y1,y2,m,n,p0,p1,p2,f) \
mcdetector_out_2D(t,xl,yl,x1,x2,y1,y2,m,n,p0,p1,p2,f,mccompcurname)
#ifdef MC_TRACE_ENABLED
#define DEBUG
#endif
#ifdef DEBUG
#define mcDEBUG_INSTR() if(!mcdotrace); else printf("INSTRUMENT:\n");
#define mcDEBUG_COMPONENT(name,c,t) if(!mcdotrace); else \
printf("COMPONENT: \"%s\"\n" \
"POS: %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g\n", \
name, c.x, c.y, c.z, t[0][0], t[0][1], t[0][2], \
t[1][0], t[1][1], t[1][2], t[2][0], t[2][1], t[2][2]);
#define mcDEBUG_INSTR_END() if(!mcdotrace); else printf("INSTRUMENT END:\n");
#define mcDEBUG_ENTER() if(!mcdotrace); else printf("ENTER:\n");
#define mcDEBUG_COMP(c) if(!mcdotrace); else printf("COMP: \"%s\"\n", c);
#define mcDEBUG_STATE(x,y,z,vx,vy,vz,t,s1,s2,p) if(!mcdotrace); else \
printf("STATE: %g, %g, %g, %g, %g, %g, %g, %g, %g, %g\n", \
x,y,z,vx,vy,vz,t,s1,s2,p);
#define mcDEBUG_SCATTER(x,y,z,vx,vy,vz,t,s1,s2,p) if(!mcdotrace); else \
printf("SCATTER: %g, %g, %g, %g, %g, %g, %g, %g, %g, %g\n", \
x,y,z,vx,vy,vz,t,s1,s2,p);
#define mcDEBUG_LEAVE() if(!mcdotrace); else printf("LEAVE:\n");
#define mcDEBUG_ABSORB() if(!mcdotrace); else printf("ABSORB:\n");
#else
#define mcDEBUG_INSTR()
#define mcDEBUG_COMPONENT(name,c,t)
#define mcDEBUG_INSTR_END()
#define mcDEBUG_ENTER()
#define mcDEBUG_COMP(c)
#define mcDEBUG_STATE(x,y,z,vx,vy,vz,t,s1,s2,p)
#define mcDEBUG_SCATTER(x,y,z,vx,vy,vz,t,s1,s2,p)
#define mcDEBUG_LEAVE()
#define mcDEBUG_ABSORB()
#endif
#ifdef TEST
#define test_printf printf
#else
#define test_printf while(0) printf
#endif
void mcdis_magnify(char *);
void mcdis_line(double, double, double, double, double, double);
void mcdis_multiline(int, ...);
void mcdis_circle(char *, double, double, double, double);
#define RAD2MIN ((180*60)/PI)
#define MIN2RAD (PI/(180*60))
#define DEG2RAD (PI/180)
#define RAD2DEG (180/PI)
#define AA2MS 629.719 /* Convert k[1/AA] to v[m/s] */
#define MS2AA 1.58801E-3 /* Convert v[m/s] to k[1/AA] */
#define K2V AA2MS
#define V2K MS2AA
#define Q2V AA2MS
#define V2Q MS2AA
#define SE2V 437.3949 /* Convert sqrt(E)[meV] to v[m/s] */
#define VS2E 5.227e-6 /* Convert (v[m/s])**2 to E[meV] */
#define FWHM2RMS 0.424660900144 /* Convert between full-width-half-max and */
#define RMS2FWHM 2.35482004503 /* root-mean-square (standard deviation) */
#define HBAR 1.05459E-34
#define MNEUTRON 1.67492E-27
#ifndef PI
# ifdef M_PI
# define PI M_PI
# else
# define PI 3.14159265358979323846
# endif
#endif
typedef int mc_int32_t;
mc_int32_t mc_random(void);
void mc_srandom (unsigned int x);
unsigned long mt_random(void);
void mt_srandom (unsigned long x);
#ifndef MC_RAND_ALG
#define MC_RAND_ALG 1
#endif
#if MC_RAND_ALG == 0
/* Use system random() (not recommended). */
# define MC_RAND_MAX RAND_MAX
#elif MC_RAND_ALG == 1
/* "Mersenne Twister", by Makoto Matsumoto and Takuji Nishimura. */
# define MC_RAND_MAX ((unsigned long)0xffffffff)
# define random mt_random
# define srandom mt_srandom
#elif MC_RAND_ALG == 2
/* Algorithm used in McStas 1.1 and earlier (not recommended). */
# define MC_RAND_MAX 0x7fffffff
# define random mc_random
# define srandom mc_srandom
#else
# error "Bad value for random number generator choice."
#endif
#define rand01() ( ((double)random())/((double)MC_RAND_MAX+1) )
#define randpm1() ( ((double)random()) / (((double)MC_RAND_MAX+1)/2) - 1 )
#define rand0max(max) ( ((double)random()) / (((double)MC_RAND_MAX+1)/(max)) )
#define randminmax(min,max) ( rand0max((max)-(min)) + (min) )
#define PROP_X0 \
do { \
double mc_dt; \
if(mcnlvx == 0) ABSORB; \
mc_dt = -mcnlx/mcnlvx; \
if(mc_dt < 0) ABSORB; \
mcnly += mcnlvy*mc_dt; \
mcnlz += mcnlvz*mc_dt; \
mcnlt += mc_dt; \
mcnlx = 0; \
} while(0)
#define PROP_Y0 \
do { \
double mc_dt; \
if(mcnlvy == 0) ABSORB; \
mc_dt = -mcnly/mcnlvy; \
if(mc_dt < 0) ABSORB; \
mcnlx += mcnlvx*mc_dt; \
mcnlz += mcnlvz*mc_dt; \
mcnlt += mc_dt; \
mcnly = 0; \
} while(0)
#define PROP_Z0 \
do { \
double mc_dt; \
if(mcnlvz == 0) ABSORB; \
mc_dt = -mcnlz/mcnlvz; \
if(mc_dt < 0) ABSORB; \
mcnlx += mcnlvx*mc_dt; \
mcnly += mcnlvy*mc_dt; \
mcnlt += mc_dt; \
mcnlz = 0; \
} while(0)
#define PROP_DT(dt) \
do { \
if(dt < 0) ABSORB; \
mcnlx += mcnlvx*(dt); \
mcnly += mcnlvy*(dt); \
mcnlz += mcnlvz*(dt); \
mcnlt += (dt); \
} while(0)
#define vec_prod(x, y, z, x1, y1, z1, x2, y2, z2) \
do { \
double mcvp_tmpx, mcvp_tmpy, mcvp_tmpz; \
mcvp_tmpx = (y1)*(z2) - (y2)*(z1); \
mcvp_tmpy = (z1)*(x2) - (z2)*(x1); \
mcvp_tmpz = (x1)*(y2) - (x2)*(y1); \
(x) = mcvp_tmpx; (y) = mcvp_tmpy; (z) = mcvp_tmpz; \
} while(0)
#define scalar_prod(x1, y1, z1, x2, y2, z2) \
((x1)*(x2) + (y1)*(y2) + (z1)*(z2))
#define NORM(x,y,z) \
do { \
double mcnm_tmp = sqrt((x)*(x) + (y)*(y) + (z)*(z)); \
if(mcnm_tmp != 0.0) \
{ \
(x) /= mcnm_tmp; \
(y) /= mcnm_tmp; \
(z) /= mcnm_tmp; \
} \
} while(0)
#define rotate(x, y, z, vx, vy, vz, phi, ax, ay, az) \
do { \
double mcrt_tmpx = (ax), mcrt_tmpy = (ay), mcrt_tmpz = (az); \
double mcrt_vp, mcrt_vpx, mcrt_vpy, mcrt_vpz; \
double mcrt_vnx, mcrt_vny, mcrt_vnz, mcrt_vn1x, mcrt_vn1y, mcrt_vn1z; \
double mcrt_bx, mcrt_by, mcrt_bz, mcrt_v1x, mcrt_v1y, mcrt_v1z; \
double mcrt_cos, mcrt_sin; \
NORM(mcrt_tmpx, mcrt_tmpy, mcrt_tmpz); \
mcrt_vp = scalar_prod((vx), (vy), (vz), mcrt_tmpx, mcrt_tmpy, mcrt_tmpz); \
mcrt_vpx = mcrt_vp*mcrt_tmpx; \
mcrt_vpy = mcrt_vp*mcrt_tmpy; \
mcrt_vpz = mcrt_vp*mcrt_tmpz; \
mcrt_vnx = (vx) - mcrt_vpx; \
mcrt_vny = (vy) - mcrt_vpy; \
mcrt_vnz = (vz) - mcrt_vpz; \
vec_prod(mcrt_bx, mcrt_by, mcrt_bz, \
mcrt_tmpx, mcrt_tmpy, mcrt_tmpz, mcrt_vnx, mcrt_vny, mcrt_vnz); \
mcrt_cos = cos((phi)); mcrt_sin = sin((phi)); \
mcrt_vn1x = mcrt_vnx*mcrt_cos + mcrt_bx*mcrt_sin; \
mcrt_vn1y = mcrt_vny*mcrt_cos + mcrt_by*mcrt_sin; \
mcrt_vn1z = mcrt_vnz*mcrt_cos + mcrt_bz*mcrt_sin; \
(x) = mcrt_vpx + mcrt_vn1x; \
(y) = mcrt_vpy + mcrt_vn1y; \
(z) = mcrt_vpz + mcrt_vn1z; \
} while(0)
Coords coords_set(MCNUM x, MCNUM y, MCNUM z);
Coords coords_add(Coords a, Coords b);
Coords coords_sub(Coords a, Coords b);
Coords coords_neg(Coords a);
void rot_set_rotation(Rotation t, double phx, double phy, double phz);
void rot_mul(Rotation t1, Rotation t2, Rotation t3);
void rot_copy(Rotation dest, Rotation src);
void rot_transpose(Rotation src, Rotation dst);
Coords rot_apply(Rotation t, Coords a);
void mccoordschange(Coords a, Rotation t, double *x, double *y, double *z,
double *vx, double *vy, double *vz, double *time,
double *s1, double *s2);
void mccoordschange_polarisation(Rotation t,
double *sx, double *sy, double *sz);
double mcestimate_error(int N, double p1, double p2);
void mcsiminfo_out(char *format, ...);
void mcdetector_out(char *cname, double p0, double p1, double p2,
char *filename);
void mcdetector_out_0D(char *t, double p0, double p1, double p2, char *cname);
void mcdetector_out_1D(char *t, char *xl, char *yl,
char *xvar, double x1, double x2, int n,
int *p0, double *p1, double *p2, char *f, char *c);
void mcdetector_out_2D(char *t, char *xl, char *yl, double x1, double x2,
double y1,double y2,int m, int n,
int *p0, double *p1, double *p2, char *f, char *c);
void mcreadparams(void);
void mcsetstate(double x, double y, double z, double vx, double vy, double vz,
double t, double sx, double sy, double sz, double p);
void mcgenstate(void);
double randnorm(void);
void normal_vec(double *nx, double *ny, double *nz, double x, double y, double z);
int box_intersect(double *dt_in, double *dt_out, double x, double y, double z,
double vx, double vy, double vz, double dx, double dy, double dz);
int cylinder_intersect(double *t0, double *t1, double x, double y, double z,
double vx, double vy, double vz, double r, double h);
int sphere_intersect(double *t0, double *t1, double x, double y, double z,
double vx, double vy, double vz, double r);
void randvec_target_sphere(double *xo, double *yo, double *zo, double *solid_angle,
double xi, double yi, double zi, double radius);
void extend_list(int count, void **list, int *size, size_t elemsize);
void mcset_ncount(long count);
long mcget_ncount(void);
long mcget_run_num(void);
int mcstas_main(int argc, char *argv[]);
#endif /* MCSTAS_R_H */
-------------- next part --------------
/*******************************************************************************
* Runtime system for McStas.
*
* Project: Monte Carlo Simulation of Triple Axis Spectrometers
* File name: mcstas-r.c
*
* Author: K.N. Aug 27, 1997
*
* Copyright (C) Risoe National Laboratory, 1997-1999, All rights reserved
*******************************************************************************/
#include <stdarg.h>
#include <limits.h>
#include <math.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#ifndef MCSTAS_R_H
#include "mcstas-r.h"
#endif
#ifdef MC_ANCIENT_COMPATIBILITY
int mctraceenabled = 0;
int mcdefaultmain = 0;
#endif
static long mcseed = 0;
mcstatic int mcdotrace = 0;
static int mcascii_only = 0;
static int mcdisable_output_files = 0;
static int mcsingle_file = 0;
static FILE *mcsiminfo_file = NULL;
static char *mcdirname = NULL;
static char *mcsiminfo_name = "mcstas.sim";
/* MCDISPLAY support. */
void mcdis_magnify(char *what){
printf("MCDISPLAY: magnify('%s')\n", what);
}
void mcdis_line(double x1, double y1, double z1,
double x2, double y2, double z2){
printf("MCDISPLAY: multiline(2,%g,%g,%g,%g,%g,%g)\n",
x1,y1,z1,x2,y2,z2);
}
void mcdis_multiline(int count, ...){
va_list ap;
double x,y,z;
printf("MCDISPLAY: multiline(%d", count);
va_start(ap, count);
while(count--)
{
x = va_arg(ap, double);
y = va_arg(ap, double);
z = va_arg(ap, double);
printf(",%g,%g,%g", x, y, z);
}
va_end(ap);
printf(")\n");
}
void mcdis_circle(char *plane, double x, double y, double z, double r){
printf("MCDISPLAY: circle('%s',%g,%g,%g,%g)\n", plane, x, y, z, r);
}
/* Assign coordinates. */
Coords
coords_set(MCNUM x, MCNUM y, MCNUM z)
{
Coords a;
a.x = x;
a.y = y;
a.z = z;
return a;
}
/* Add two coordinates. */
Coords
coords_add(Coords a, Coords b)
{
Coords c;
c.x = a.x + b.x;
c.y = a.y + b.y;
c.z = a.z + b.z;
return c;
}
/* Subtract two coordinates. */
Coords
coords_sub(Coords a, Coords b)
{
Coords c;
c.x = a.x - b.x;
c.y = a.y - b.y;
c.z = a.z - b.z;
return c;
}
/* Negate coordinates. */
Coords
coords_neg(Coords a)
{
Coords b;
b.x = -a.x;
b.y = -a.y;
b.z = -a.z;
return b;
}
/*******************************************************************************
* Get transformation for rotation first phx around x axis, then phy around y,
* then phz around z.
*******************************************************************************/
void
rot_set_rotation(Rotation t, double phx, double phy, double phz)
{
double cx = cos(phx);
double sx = sin(phx);
double cy = cos(phy);
double sy = sin(phy);
double cz = cos(phz);
double sz = sin(phz);
t[0][0] = cy*cz;
t[0][1] = sx*sy*cz + cx*sz;
t[0][2] = sx*sz - cx*sy*cz;
t[1][0] = -cy*sz;
t[1][1] = cx*cz - sx*sy*sz;
t[1][2] = sx*cz + cx*sy*sz;
t[2][0] = sy;
t[2][1] = -sx*cy;
t[2][2] = cx*cy;
}
/*******************************************************************************
* Matrix multiplication of transformations (this corresponds to combining
* transformations). After rot_mul(T1, T2, T3), doing T3 is equal to doing
* first T2, then T1.
* Note that T3 must not alias (use the same array as) T1 or T2.
*******************************************************************************/
void
rot_mul(Rotation t1, Rotation t2, Rotation t3)
{
int i,j, k;
for(i = 0; i < 3; i++)
for(j = 0; j < 3; j++)
t3[i][j] = t1[i][0]*t2[0][j] + t1[i][1]*t2[1][j] + t1[i][2]*t2[2][j];
}
/*******************************************************************************
* Copy a rotation transformation (needed since arrays cannot be assigned in C).
*******************************************************************************/
void
rot_copy(Rotation dest, Rotation src)
{
dest[0][0] = src[0][0];
dest[0][1] = src[0][1];
dest[0][2] = src[0][2];
dest[1][0] = src[1][0];
dest[1][1] = src[1][1];
dest[1][2] = src[1][2];
dest[2][0] = src[2][0];
dest[2][1] = src[2][1];
dest[2][2] = src[2][2];
}
void
rot_transpose(Rotation src, Rotation dst)
{
dst[0][0] = src[0][0];
dst[0][1] = src[1][0];
dst[0][2] = src[2][0];
dst[1][0] = src[0][1];
dst[1][1] = src[1][1];
dst[1][2] = src[2][1];
dst[2][0] = src[0][2];
dst[2][1] = src[1][2];
dst[2][2] = src[2][2];
}
Coords
rot_apply(Rotation t, Coords a)
{
Coords b;
b.x = t[0][0]*a.x + t[0][1]*a.y + t[0][2]*a.z;
b.y = t[1][0]*a.x + t[1][1]*a.y + t[1][2]*a.z;
b.z = t[2][0]*a.x + t[2][1]*a.y + t[2][2]*a.z;
return b;
}
void
mccoordschange(Coords a, Rotation t, double *x, double *y, double *z,
double *vx, double *vy, double *vz, double *time,
double *s1, double *s2)
{
Coords b, c;
b.x = *x;
b.y = *y;
b.z = *z;
c = rot_apply(t, b);
b = coords_add(c, a);
*x = b.x;
*y = b.y;
*z = b.z;
b.x = *vx;
b.y = *vy;
b.z = *vz;
c = rot_apply(t, b);
*vx = c.x;
*vy = c.y;
*vz = c.z;
/* ToDo: What to do about the spin? */
}
void
mccoordschange_polarisation(Rotation t, double *sx, double *sy, double *sz)
{
Coords b, c;
b.x = *sx;
b.y = *sy;
b.z = *sz;
c = rot_apply(t, b);
*sx = c.x;
*sy = c.y;
*sz = c.z;
}
/*******************************************************************************
* Find i in adaptive search tree t s.t. v(i) <= v < v(i+1).
*******************************************************************************/
int
adapt_tree_search(struct adapt_tree *t, adapt_t v)
{
adapt_t F = 0; /* Current value. */
int i = 0; /* Current candidate. */
int step = t->initstep;
adapt_t *s = t->s;
int j;
for(j = t->root; step > 0; step >>= 1)
{
F += s[j]; /* Cumulative value in current node */
if(v < F)
j -= step; /* Value is to the left or above. */
else
i = j, j += step; /* Value is current or to the right. */
}
/* Now j is at the bottom of a tree (a leaf node). */
if(v < F + s[j])
return i;
else
return j;
}
/*******************************************************************************
* Add v to v[i], updating the cumulative sums appropriately.
*******************************************************************************/
void
adapt_tree_add(struct adapt_tree *t, int i, adapt_t v)
{
int j = t->root;
int step = t->initstep;
adapt_t *s = t->s;
t->total += v;
t->v[i++] += v;
for(;;)
{
while(j < i)
j += step, step >>= 1;
s[j] += v;
while(j > i)
j -= step, step >>= 1;
if(j == i)
break;
s[j] -= v;
}
if(step)
s[j - step] -= v;
}
/*******************************************************************************
* Initialise an adaptive search tree. The tree has N nodes, and all nodes are
* initialized to zero. Any N > 0 is allowed, but is rounded up to the nearest
* value of the form N = 2**k - 2.
*******************************************************************************/
struct adapt_tree *
adapt_tree_init(int N)
{
struct adapt_tree *t;
int i;
int depth;
/* Round up to nearest 2**k - 2 */
for(depth = 0; ((1 << (depth + 1)) - 2) < N; depth++);
N = (1 << (depth + 1)) - 2;
t = malloc(sizeof(*t));
if(t)
{
t->s = malloc((N + 1) * sizeof(*(t->s)));
t->v = malloc(N * sizeof(*(t->v)));
}
if(!(t && t->s && t->v))
{
fprintf(stderr, "Fatal error: out of memory\n");
exit(1);
}
t->N = N;
t->depth = depth;
t->root = (1 << t->depth) - 1;
t->initstep = (1 << (t->depth - 1));
for(i = 0; i < t->N; i++)
{
t->s[i] = 0.0;
t->v[i] = 0.0;
}
t->s[i] = 0.0;
t->total = 0.0;
return t;
}
/*******************************************************************************
* Free memory allocated to an adaptive search tree.
*******************************************************************************/
void
adapt_tree_free(struct adapt_tree *t)
{
free(t->v);
free(t->s);
free(t);
}
double
mcestimate_error(int N, double p1, double p2)
{
double pmean, n1;
if(N <= 1)
return p1;
pmean = p1 / N;
n1 = N - 1;
/* Note: underflow may cause p2 to become zero; the fabs() below guards
against this. */
return sqrt((N/n1)*fabs(p2 - pmean*pmean));
}
/* Instrument input parameter type handling. */
static int
mcparm_double(char *s, void *vptr)
{
char *p;
double *v = vptr;
*v = strtod(s, &p);
if(*s == '\0' || (p != NULL && *p != '\0') || errno == ERANGE)
return 0; /* Failed */
else
return 1; /* Success */
}
static char *
mcparminfo_double(char *parmname)
{
return "double";
}
static void
mcparmerror_double(char *parm, char *val)
{
fprintf(stderr, "Error: Invalid value '%s' for floating point parameter %s\n",
val, parm);
}
static void
mcparmprinter_double(FILE *f, void *vptr)
{
double *v = vptr;
fprintf(f, "%g", *v);
}
static int
mcparm_int(char *s, void *vptr)
{
char *p;
int *v = vptr;
long x;
*v = 0;
x = strtol(s, &p, 10);
if(x < INT_MIN || x > INT_MAX)
return 0; /* Under/overflow */
*v = x;
if(*s == '\0' || (p != NULL && *p != '\0') || errno == ERANGE)
return 0; /* Failed */
else
return 1; /* Success */
}
static char *
mcparminfo_int(char *parmname)
{
return "int";
}
static void
mcparmerror_int(char *parm, char *val)
{
fprintf(stderr, "Error: Invalid value '%s' for integer parameter %s\n",
val, parm);
}
static void
mcparmprinter_int(FILE *f, void *vptr)
{
int *v = vptr;
fprintf(f, "%d", *v);
}
static int
mcparm_string(char *s, void *vptr)
{
char **v = vptr;
*v = malloc(strlen(s) + 1);
if(*v == NULL)
{
fprintf(stderr, "mcparm_string: Out of memory.\n");
exit(1);
}
strcpy(*v, s);
return 1; /* Success */
}
static char *
mcparminfo_string(char *parmname)
{
return "string";
}
static void
mcparmerror_string(char *parm, char *val)
{
fprintf(stderr, "Error: Invalid value '%s' for string parameter %s\n",
val, parm);
}
static void
mcparmprinter_string(FILE *f, void *vptr)
{
char **v = vptr;
char *p;
fprintf(f, "\"");
for(p = *v; *p != '\0'; p++)
{
switch(*p)
{
case '\n':
fprintf(f, "\\n");
break;
case '\r':
fprintf(f, "\\r");
break;
case '"':
fprintf(f, "\\\"");
break;
case '\\':
fprintf(f, "\\\\");
break;
default:
fprintf(f, "%c", *p);
}
}
fprintf(f, "\"");
}
static struct
{
int (*getparm)(char *, void *);
char * (*parminfo)(char *);
void (*error)(char *, char *);
void (*printer)(FILE *, void *);
} mcinputtypes[] =
{
mcparm_double, mcparminfo_double, mcparmerror_double,
mcparmprinter_double,
mcparm_int, mcparminfo_int, mcparmerror_int,
mcparmprinter_int,
mcparm_string, mcparminfo_string, mcparmerror_string,
mcparmprinter_string
};
FILE *
mcnew_file(char *name)
{
int dirlen;
char *mem;
FILE *file;
dirlen = mcdirname ? strlen(mcdirname) : 0;
mem = malloc(dirlen + 1 + strlen(name) + 1);
if(!mem)
{
fprintf(stderr, "Error: Out of memory\n");
exit(1);
}
strcpy(mem, "");
if(dirlen)
{
strcat(mem, mcdirname);
if(mcdirname[dirlen - 1] != MC_PATHSEP_C &&
name[0] != MC_PATHSEP_C)
strcat(mem, MC_PATHSEP_S);
}
strcat(mem, name);
file = fopen(mem, "w");
if(!file)
fprintf(stderr, "Warning: could not open output file '%s'\n", mem);
free(mem);
return file;
}
static void
mcinfo_out(char *pre, FILE *f)
{
int i;
fprintf(f, "%sName: %s\n", pre, mcinstrument_name);
fprintf(f, "%sParameters:", pre);
for(i = 0; i < mcnumipar; i++)
fprintf(f, " %s(%s)", mcinputtable[i].name,
(*mcinputtypes[mcinputtable[i].type].parminfo)
(mcinputtable[i].name));
fprintf(f, "\n");
fprintf(f, "%sInstrument-source: %s\n", pre, mcinstrument_source);
fprintf(f, "%sTrace-enabled: %s\n", pre, mctraceenabled ? "yes" : "no");
fprintf(f, "%sDefault-main: %s\n", pre, mcdefaultmain ? "yes" : "no");
fprintf(f, "%sEmbedded-runtime: %s\n", pre,
#ifdef MC_EMBEDDED_RUNTIME
"yes"
#else
"no"
#endif
);
}
static void
mcruninfo_out(char *pre, FILE *f)
{
int i;
time_t t;
if(!f)
return;
time(&t);
fprintf(f, "%sDate: %s", pre, ctime(&t)); /* Note: ctime adds '\n' ! */
fprintf(f, "%sNcount: %i\n", pre, mcget_ncount());
fprintf(f, "%sTrace: %s\n", pre, mcdotrace ? "yes" : "no");
if(mcseed)
fprintf(f, "%sSeed: %ld\n", pre, mcseed);
for(i = 0; i < mcnumipar; i++)
{
fprintf(f, "%sParam: %s=", pre, mcinputtable[i].name);
(*mcinputtypes[mcinputtable[i].type].printer)(f, mcinputtable[i].par);
fprintf(f, "\n");
}
}
void
mcsiminfo_init(void)
{
if(mcdisable_output_files)
return;
mcsiminfo_file = mcnew_file(mcsiminfo_name);
if(!mcsiminfo_file)
fprintf(stderr,
"Warning: could not open simulation description file '%s'\n",
mcsiminfo_name);
else
{
mcsiminfo_out("begin instrument\n");
mcinfo_out(" ", mcsiminfo_file);
mcsiminfo_out("end instrument\n");
mcsiminfo_out("\nbegin simulation\n");
mcruninfo_out(" ", mcsiminfo_file);
mcsiminfo_out("end simulation\n");
}
}
void
mcsiminfo_out(char *format, ...)
{
va_list ap;
double x,y,z;
if(mcsiminfo_file)
{
va_start(ap, format);
vfprintf(mcsiminfo_file, format, ap);
va_end(ap);
}
}
void
mcsiminfo_close()
{
if(mcsiminfo_file)
fclose(mcsiminfo_file);
}
void
mcdetector_out(char *cname, double p0, double p1, double p2, char *filename)
{
printf("Detector: %s_I=%g %s_ERR=%g %s_N=%g",
cname, p1, cname, mcestimate_error(p0,p1,p2), cname, p0);
if(filename)
printf(" \"%s\"", filename);
printf("\n");
}
static void
mcdatainfo_out_0D(char *pre, FILE *f, char *t,
double I, double I_err, double N, char *c)
{
if(!f)
return;
fprintf(f, "%stype: array_0d\n", pre);
fprintf(f, "%scomponent: %s\n", pre, c);
fprintf(f, "%stitle: %s\n", pre, t);
fprintf(f, "%svariables: I I_err N\n", pre);
fprintf(f, "%svalues: %g %g %g\n", pre, I, I_err, N);
}
static void
mcdatainfo_out_1D(char *pre, FILE *f, char *t, char *xl, char *yl, char *xvar,
double x1, double x2, int n, char *fname, char *c, int do_errb)
{
if(!f)
return;
fprintf(f, "%stype: array_1d(%d)\n", pre, n);
fprintf(f, "%scomponent: %s\n", pre, c);
fprintf(f, "%stitle: %s\n", pre, t);
if(fname)
fprintf(f, "%sfilename: '%s'\n", pre, fname);
fprintf(f, "%svariables: %s I%s\n", pre, xvar, do_errb ? " I_err N" : "");
fprintf(f, "%sxvar: %s\n", pre, xvar);
fprintf(f, "%syvar: %s\n", pre, do_errb ? "(I,I_err)" : "I");
fprintf(f, "%sxlabel: '%s'\n%sylabel: '%s'\n", pre, xl, pre, yl);
fprintf(f, "%sxlimits: %g %g\n", pre, x1, x2);
}
static void
mcdatainfo_out_2D(char *pre, FILE *f, char *t, char *xl, char *yl,
double x1, double x2, double y1, double y2,
int m, int n, char *fname, char *c)
{
if(!f)
return;
fprintf(f, "%stype: array_2d(%d,%d)\n", pre, m, n);
fprintf(f, "%scomponent: %s\n", pre, c);
fprintf(f, "%stitle: %s\n", pre, t);
if(fname)
fprintf(f, "%sfilename: '%s'\n", pre, fname);
fprintf(f, "%sxlabel: '%s'\n%sylabel: '%s'\n", pre, xl, pre, yl);
fprintf(f, "%sxylimits: %g %g %g %g\n", pre, x1, x2, y1, y2);
}
/*******************************************************************************
* Output single detector/monitor data (p0, p1, p2).
* Title is t, component name is c.
*******************************************************************************/
void
mcdetector_out_0D(char *t, double p0, double p1, double p2, char *c)
{
/* Write data set information to simulation description file. */
mcsiminfo_out("\nbegin data\n");
mcdatainfo_out_0D(" ", mcsiminfo_file, t,
p1, mcestimate_error(p0, p1, p2), p0, c);
mcsiminfo_out("end data\n");
/* Finally give 0D detector output. */
mcdetector_out(c, p0, p1, p2, NULL);
}
/*******************************************************************************
* Output 1d detector data (p0, p1, p2) for n bins linearly
* distributed across the range x1..x2 (x1 is lower limit of first
* bin, x2 is upper limit of last bin). Title is t, axis labels are xl
* and yl. File name is f, component name is c.
*******************************************************************************/
void
mcdetector_out_1D(char *t, char *xl, char *yl,
char *xvar, double x1, double x2, int n,
int *p0, double *p1, double *p2, char *f, char *c)
{
int i;
FILE *outfile = NULL;
int Nsum;
double Psum, P2sum;
char *pre;
int do_errb = p0 && p2;
/* Write data set information to simulation description file. */
mcsiminfo_out("\nbegin data\n");
mcdatainfo_out_1D(" ", mcsiminfo_file, t, xl, yl, xvar, x1, x2,
n, f, c, do_errb);
/* Loop over array elements, computing total sums and writing to file. */
Nsum = Psum = P2sum = 0;
pre = "";
if(mcsingle_file)
{
outfile = mcsiminfo_file;
f = NULL;
mcsiminfo_out(" begin array2D (%d,%d)\n", do_errb ? 4 : 2, n);
pre = " ";
}
else if(f && !mcdisable_output_files) /* Don't write if filename is NULL */
outfile = mcnew_file(f);
if(outfile && !mcascii_only && !mcsingle_file)
{
fprintf(outfile, "# Instrument-source: %s\n", mcinstrument_source);
mcruninfo_out("# ", outfile);
mcdatainfo_out_1D("# ", outfile, t, xl, yl, xvar, x1, x2,
n, f, c, do_errb);
}
for(i = 0; i < n; i++)
{
if(outfile)
{
fprintf(outfile, "%s%g %g", pre, x1 + (i + 0.5)/n*(x2 - x1), p1[i]);
if(do_errb)
fprintf(outfile, " %g %d", mcestimate_error(p0[i],p1[i],p2[i]), p0[i]);
fprintf(outfile, "\n");
}
Nsum += p0 ? p0[i] : 1;
Psum += p1[i];
P2sum += p2 ? p2[i] : p1[i]*p1[i];
}
if(mcsingle_file)
mcsiminfo_out(" end array2D\n");
else if(outfile)
fclose(outfile);
mcsiminfo_out("end data\n");
/* Finally give 0D detector output. */
mcdetector_out(c, Nsum, Psum, P2sum, f);
}
void
mcdetector_out_2D(char *t, char *xl, char *yl,
double x1, double x2, double y1, double y2, int m,
int n, int *p0, double *p1, double *p2, char *f, char *c)
{
int i, j;
FILE *outfile = NULL;
int Nsum;
double Psum, P2sum;
char *pre;
/* Write data set information to simulation description file. */
mcsiminfo_out("\nbegin data\n");
mcdatainfo_out_2D(" ", mcsiminfo_file,
t, xl, yl, x1, x2, y1, y2, m, n, f, c);
/* Loop over array elements, computing total sums and writing to file. */
Nsum = Psum = P2sum = 0;
pre = "";
if(mcsingle_file)
{
outfile = mcsiminfo_file;
f = NULL;
mcsiminfo_out(" begin array2D (%d,%d)\n", m, n);
pre = " ";
}
else if(f && !mcdisable_output_files) /* Don't write if filename is NULL */
outfile = mcnew_file(f);
if(outfile && !mcascii_only && !mcsingle_file)
{
fprintf(outfile, "# Instrument-source: %s\n", mcinstrument_source);
mcruninfo_out("# ", outfile);
mcdatainfo_out_2D("# ", outfile,
t, xl, yl, x1, x2, y1, y2, m, n, f, c);
}
for(j = 0; j < n; j++)
{
if(outfile)
fprintf(outfile,"%s", pre);
for(i = 0; i < m; i++)
{
if(outfile)
fprintf(outfile, "%g ", p1[i*n + j]);
Nsum += p0 ? p0[i*n + j] : 1;
Psum += p1[i*n + j];
P2sum += p2 ? p2[i*n + j] : p1[i*n + j]*p1[i*n + j];
}
if(outfile)
fprintf(outfile,"\n");
}
if(mcsingle_file)
mcsiminfo_out(" end array2D\n");
else if(outfile)
fclose(outfile);
mcsiminfo_out("end data\n");
/* Finally give 0D detector output. */
mcdetector_out(c, Nsum, Psum, P2sum, f);
}
void
mcreadparams(void)
{
int i,j,status;
char buf[1024];
char *p;
int len;
for(i = 0; mcinputtable[i].name != 0; i++)
{
do
{
printf("Set value of instrument parameter %s (%s):\n",
mcinputtable[i].name,
(*mcinputtypes[mcinputtable[i].type].parminfo)
(mcinputtable[i].name));
fflush(stdout);
p = fgets(buf, 1024, stdin);
if(p == NULL)
{
fprintf(stderr, "Error: empty input\n");
exit(1);
}
len = strlen(buf);
for(j = 0; j < 2; j++)
{
if(len > 0 && (buf[len - 1] == '\n' || buf[len - 1] == '\r'))
{
len--;
buf[len] = '\0';
}
}
status = (*mcinputtypes[mcinputtable[i].type].getparm)
(buf, mcinputtable[i].par);
if(!status)
{
(*mcinputtypes[mcinputtable[i].type].error)(mcinputtable[i].name, buf);
}
} while(!status);
}
}
void
mcsetstate(double x, double y, double z, double vx, double vy, double vz,
double t, double sx, double sy, double sz, double p)
{
extern double mcnx, mcny, mcnz, mcnvx, mcnvy, mcnvz;
extern double mcnt, mcnsx, mcnsy, mcnsz, mcnp;
mcnx = x;
mcny = y;
mcnz = z;
mcnvx = vx;
mcnvy = vy;
mcnvz = vz;
mcnt = t;
mcnsx = sx;
mcnsy = sy;
mcnsz = sz;
mcnp = p;
}
void
mcgenstate(void)
{
mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1);
}
/* McStas random number routine. */
/*
* Copyright (c) 1983 Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms are permitted
* provided that the above copyright notice and this paragraph are
* duplicated in all such forms and that any documentation,
* advertising materials, and other materials related to such
* distribution and use acknowledge that the software was developed
* by the University of California, Berkeley. The name of the
* University may not be used to endorse or promote products derived
* from this software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*/
/*
* This is derived from the Berkeley source:
* @(#)random.c 5.5 (Berkeley) 7/6/88
* It was reworked for the GNU C Library by Roland McGrath.
* Rewritten to use reentrant functions by Ulrich Drepper, 1995.
*/
/*******************************************************************************
* Modified for McStas from glibc 2.0.7pre1 stdlib/random.c and
* stdlib/random_r.c.
*
* This way random() is more than four times faster compared to calling
* standard glibc random() on ix86 Linux, probably due to multithread support,
* ELF shared library overhead, etc. It also makes McStas generated
* simulations more portable (more likely to behave identically across
* platforms, important for parrallel computations).
*******************************************************************************/
#define TYPE_3 3
#define BREAK_3 128
#define DEG_3 31
#define SEP_3 3
static mc_int32_t randtbl[DEG_3 + 1] =
{
TYPE_3,
-1726662223, 379960547, 1735697613, 1040273694, 1313901226,
1627687941, -179304937, -2073333483, 1780058412, -1989503057,
-615974602, 344556628, 939512070, -1249116260, 1507946756,
-812545463, 154635395, 1388815473, -1926676823, 525320961,
-1009028674, 968117788, -123449607, 1284210865, 435012392,
-2017506339, -911064859, -370259173, 1132637927, 1398500161,
-205601318,
};
static mc_int32_t *fptr = &randtbl[SEP_3 + 1];
static mc_int32_t *rptr = &randtbl[1];
static mc_int32_t *state = &randtbl[1];
#define rand_deg DEG_3
#define rand_sep SEP_3
static mc_int32_t *end_ptr = &randtbl[sizeof (randtbl) / sizeof (randtbl[0])];
mc_int32_t
mc_random (void)
{
mc_int32_t result;
*fptr += *rptr;
/* Chucking least random bit. */
result = (*fptr >> 1) & 0x7fffffff;
++fptr;
if (fptr >= end_ptr)
{
fptr = state;
++rptr;
}
else
{
++rptr;
if (rptr >= end_ptr)
rptr = state;
}
return result;
}
void
mc_srandom (unsigned int x)
{
/* We must make sure the seed is not 0. Take arbitrarily 1 in this case. */
state[0] = x ? x : 1;
{
long int i;
for (i = 1; i < rand_deg; ++i)
{
/* This does:
state[i] = (16807 * state[i - 1]) % 2147483647;
but avoids overflowing 31 bits. */
long int hi = state[i - 1] / 127773;
long int lo = state[i - 1] % 127773;
long int test = 16807 * lo - 2836 * hi;
state[i] = test + (test < 0 ? 2147483647 : 0);
}
fptr = &state[rand_sep];
rptr = &state[0];
for (i = 0; i < 10 * rand_deg; ++i)
random ();
}
}
/* "Mersenne Twister", by Makoto Matsumoto and Takuji Nishimura. */
/* See http://www.math.keio.ac.jp/~matumoto/emt.html for original source. */
/* Coded by Takuji Nishimura, considering the suggestions by */
/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Library General Public */
/* License as published by the Free Software Foundation; either */
/* version 2 of the License, or (at your option) any later */
/* version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */
/* See the GNU Library General Public License for more details. */
/* You should have received a copy of the GNU Library General */
/* Public License along with this library; if not, write to the */
/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */
/* 02111-1307 USA */
/* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. */
/* When you use this, send an email to: matumoto at math.keio.ac.jp */
/* with an appropriate reference to your work. */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y) (y >> 11)
#define TEMPERING_SHIFT_S(y) (y << 7)
#define TEMPERING_SHIFT_T(y) (y << 15)
#define TEMPERING_SHIFT_L(y) (y >> 18)
static unsigned long mt[N]; /* the array for the state vector */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
/* initializing the array with a NONZERO seed */
void
mt_srandom(unsigned long seed)
{
/* setting initial seeds to mt[N] using */
/* the generator Line 25 of Table 1 in */
/* [KNUTH 1981, The Art of Computer Programming */
/* Vol. 2 (2nd Ed.), pp102] */
mt[0]= seed & 0xffffffff;
for (mti=1; mti<N; mti++)
mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
}
unsigned long
mt_random(void)
{
unsigned long y;
static unsigned long mag01[2]={0x0, MATRIX_A};
/* mag01[x] = x * MATRIX_A for x=0,1 */
if (mti >= N) { /* generate N words at one time */
int kk;
if (mti == N+1) /* if sgenrand() has not been called, */
mt_srandom(4357); /* a default initial seed is used */
for (kk=0;kk<N-M;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
}
for (;kk<N-1;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
}
y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
mti = 0;
}
y = mt[mti++];
y ^= TEMPERING_SHIFT_U(y);
y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
y ^= TEMPERING_SHIFT_L(y);
return y;
}
#undef N
#undef M
#undef MATRIX_A
#undef UPPER_MASK
#undef LOWER_MASK
#undef TEMPERING_MASK_B
#undef TEMPERING_MASK_C
#undef TEMPERING_SHIFT_U
#undef TEMPERING_SHIFT_S
#undef TEMPERING_SHIFT_T
#undef TEMPERING_SHIFT_L
/* End of "Mersenne Twister". */
/* End of McStas random number routine. */
double
randnorm(void)
{
static double v1, v2, s;
static int phase = 0;
double X, u1, u2;
if(phase == 0)
{
do
{
u1 = rand01();
u2 = rand01();
v1 = 2*u1 - 1;
v2 = 2*u2 - 1;
s = v1*v1 + v2*v2;
} while(s >= 1 || s == 0);
X = v1*sqrt(-2*log(s)/s);
}
else
{
X = v2*sqrt(-2*log(s)/s);
}
phase = 1 - phase;
return X;
}
/* Compute normal vector to (x,y,z). */
void normal_vec(double *nx, double *ny, double *nz,
double x, double y, double z)
{
double ax = fabs(x);
double ay = fabs(y);
double az = fabs(z);
double l;
if(x == 0 && y == 0 && z == 0)
{
*nx = 0;
*ny = 0;
*nz = 0;
return;
}
if(ax < ay)
{
if(ax < az)
{ /* Use X axis */
l = sqrt(z*z + y*y);
*nx = 0;
*ny = z/l;
*nz = -y/l;
return;
}
}
else
{
if(ay < az)
{ /* Use Y axis */
l = sqrt(z*z + x*x);
*nx = z/l;
*ny = 0;
*nz = -x/l;
return;
}
}
/* Use Z axis */
l = sqrt(y*y + x*x);
*nx = y/l;
*ny = -x/l;
*nz = 0;
}
/* If intersection with box dt_in and dt_out is returned */
/* This function written by Stine Nyborg, 1999. */
int box_intersect(double *dt_in, double *dt_out,
double x, double y, double z,
double vx, double vy, double vz,
double dx, double dy, double dz)
{
double x_in, y_in, z_in, tt, t[6], a, b;
int i, count, s;
/* Calculate intersection time for each of the six box surface planes
* If the box surface plane is not hit, the result is zero.*/
if(vx != 0)
{
tt = -(dx/2 + x)/vx;
y_in = y + tt*vy;
z_in = z + tt*vz;
if( y_in > -dy/2 && y_in < dy/2 && z_in > -dz/2 && z_in < dz/2)
t[0] = tt;
else
t[0] = 0;
tt = (dx/2 - x)/vx;
y_in = y + tt*vy;
z_in = z + tt*vz;
if( y_in > -dy/2 && y_in < dy/2 && z_in > -dz/2 && z_in < dz/2)
t[1] = tt;
else
t[1] = 0;
}
else
t[0] = t[1] = 0;
if(vy != 0)
{
tt = -(dy/2 + y)/vy;
x_in = x + tt*vx;
z_in = z + tt*vz;
if( x_in > -dx/2 && x_in < dx/2 && z_in > -dz/2 && z_in < dz/2)
t[2] = tt;
else
t[2] = 0;
tt = (dy/2 - y)/vy;
x_in = x + tt*vx;
z_in = z + tt*vz;
if( x_in > -dx/2 && x_in < dx/2 && z_in > -dz/2 && z_in < dz/2)
t[3] = tt;
else
t[3] = 0;
}
else
t[2] = t[3] = 0;
if(vz != 0)
{
tt = -(dz/2 + z)/vz;
x_in = x + tt*vx;
y_in = y + tt*vy;
if( x_in > -dx/2 && x_in < dx/2 && y_in > -dy/2 && y_in < dy/2)
t[4] = tt;
else
t[4] = 0;
tt = (dz/2 - z)/vz;
x_in = x + tt*vx;
y_in = y + tt*vy;
if( x_in > -dx/2 && x_in < dx/2 && y_in > -dy/2 && y_in < dy/2)
t[5] = tt;
else
t[5] = 0;
}
else
t[4] = t[5] = 0;
/* The intersection is evaluated and *dt_in and *dt_out are assigned */
a = b = s = 0;
count = 0;
for( i = 0; i < 6; i = i + 1 )
if( t[i] == 0 )
s = s+1;
else if( count == 0 )
{
a = t[i];
count = 1;
}
else
{
b = t[i];
count = 2;
}
if ( a == 0 && b == 0 )
return 0;
else if( a < b )
{
*dt_in = a;
*dt_out = b;
return 1;
}
else
{
*dt_in = b;
*dt_out = a;
return 1;
}
}
/* Written by: EM,NB,ABA 4.2.98 */
int
cylinder_intersect(double *t0, double *t1, double x, double y, double z,
double vx, double vy, double vz, double r, double h)
{
double v, D, t_in, t_out, y_in, y_out;
v = sqrt(vx*vx+vy*vy+vz*vz);
D = (2*vx*x + 2*vz*z)*(2*vx*x + 2*vz*z)
- 4*(vx*vx + vz*vz)*(x*x + z*z - r*r);
if (D>=0)
{
t_in = (-(2*vz*z + 2*vx*x) - sqrt(D))/(2*(vz*vz + vx*vx));
t_out = (-(2*vz*z + 2*vx*x) + sqrt(D))/(2*(vz*vz + vx*vx));
y_in = vy*t_in + y;
y_out =vy*t_out + y;
if ( (y_in > h/2 && y_out > h/2) || (y_in < -h/2 && y_out < -h/2) )
return 0;
else
{
if (y_in > h/2)
t_in = ((h/2)-y)/vy;
if (y_in < -h/2)
t_in = ((-h/2)-y)/vy;
if (y_out > h/2)
t_out = ((h/2)-y)/vy;
if (y_out < -h/2)
t_out = ((-h/2)-y)/vy;
}
*t0 = t_in;
*t1 = t_out;
return 1;
}
else
{
*t0 = *t1 = 0;
return 0;
}
}
/* Calculate intersection between line and sphere. */
int
sphere_intersect(double *t0, double *t1, double x, double y, double z,
double vx, double vy, double vz, double r)
{
double A, B, C, D, v;
v = sqrt(vx*vx + vy*vy + vz*vz);
A = v*v;
B = 2*(x*vx + y*vy + z*vz);
C = x*x + y*y + z*z - r*r;
D = B*B - 4*A*C;
if(D < 0)
return 0;
D = sqrt(D);
*t0 = (-B - D) / (2*A);
*t1 = (-B + D) / (2*A);
return 1;
}
/* Choose random direction towards target at (x,y,z) with given radius. */
/* If radius is zero, choose random direction in full 4PI, no target. */
/* ToDo: It should be possible to optimize this to avoid computing angles. */
void
randvec_target_sphere(double *xo, double *yo, double *zo, double *solid_angle,
double xi, double yi, double zi, double radius)
{
double l, theta0, phi, theta, nx, ny, nz, xt, yt, zt;
if(radius == 0.0)
{
/* No target, choose uniformly a direction in full 4PI solid angle. */
theta = acos (1 - rand0max(2));
phi = rand0max(2 * PI);
if(solid_angle)
*solid_angle = 4*PI;
nx = 1;
ny = 0;
nz = 0;
xi = 0;
yi = 1;
zi = 0;
}
else
{
l = sqrt(xi*xi + yi*yi + zi*zi); /* Distance to target. */
theta0 = asin(radius / l); /* Target size as seen from origin */
if(solid_angle)
{
/* Compute solid angle of target as seen from origin. */
*solid_angle = 2*PI*(1 - cos(theta0));
}
/* Now choose point uniformly on sphere surface within angle theta0 */
theta = acos (1 - rand0max(1 - cos(theta0)));
phi = rand0max(2 * PI);
/* Now, to obtain the desired vector rotate (x,y,z) angle theta around a
perpendicular axis (nx,ny,nz) and then angle phi around (x,y,z). */
if(xi == 0 && yi == 0)
{
nx = 1;
ny = 0;
nz = 0;
}
else
{
nx = -yi;
ny = xi;
nz = 0;
}
}
rotate(xt, yt, zt, xi, yi, zi, theta, nx, ny, nz);
rotate(*xo, *yo, *zo, xt, yt, zt, phi, xi, yi, zi);
}
/* Make sure a list is big enough to hold element COUNT.
*
* The list is an array, and the argument 'list' is a pointer to a pointer to
* the array start. The argument 'size' is a pointer to the number of elements
* in the array. The argument 'elemsize' is the sizeof() an element. The
* argument 'count' is the minimum number of elements needed in the list.
*
* If the old array is to small (or if *list is NULL or *size is 0), a
* sufficuently big new array is allocated, and *list and *size are updated.
*/
void extend_list(int count, void **list, int *size, size_t elemsize)
{
if(count >= *size)
{
void *oldlist = *list;
if(*size > 0)
*size *= 2;
else
*size = 32;
*list = malloc(*size*elemsize);
if(!*list)
{
fprintf(stderr, "\nFatal error: Out of memory.\n");
exit(1);
}
if(oldlist)
{
memcpy(*list, oldlist, count*elemsize);
free(oldlist);
}
}
}
/* Number of neutron histories to simulate. */
static long mcncount = 1e6;
long mcrun_num = 0;
void
mcset_ncount(long count)
{
mcncount = count;
}
long
mcget_ncount(void)
{
return mcncount;
}
long
mcget_run_num(void)
{
return mcrun_num;
}
static void
mcsetn_arg(char *arg)
{
mcset_ncount(strtod(arg, NULL));
}
static void
mcsetseed(char *arg)
{
mcseed = atol(arg);
if(mcseed)
srandom(mcseed);
else
{
fprintf(stderr, "Error seed most not be zero.\n");
exit(1);
}
}
static void
mchelp(char *pgmname)
{
int i;
fprintf(stderr, "Usage: %s [options] [parm=value ...]\n", pgmname);
fprintf(stderr,
"Options are:\n"
" -s SEED --seed=SEED Set random seed (must be != 0)\n"
" -n COUNT --ncount=COUNT Set number of neutrons to simulate.\n"
" -d DIR --dir=DIR Put all data files in directory DIR.\n"
" -f FILE --file=FILE Put all data in a single file.\n"
" -t --trace Enable trace of neutron through instrument.\n"
" -a --ascii-only Do not put any headers in the data files.\n"
" --no-output-files Do not write any data files.\n"
" -h --help Show help message.\n"
" -i --info Detailed instrument information.\n"
);
if(mcnumipar > 0)
{
fprintf(stderr, "Instrument parameters are:\n");
for(i = 0; i < mcnumipar; i++)
fprintf(stderr, " %-16s(%s)\n", mcinputtable[i].name,
(*mcinputtypes[mcinputtable[i].type].parminfo)(mcinputtable[i].name));
}
}
static void
mcshowhelp(char *pgmname)
{
mchelp(pgmname);
exit(0);
}
static void
mcusage(char *pgmname)
{
fprintf(stderr, "Error: incorrect command line arguments\n");
mchelp(pgmname);
exit(1);
}
static void
mcinfo(void)
{
mcinfo_out("", stdout);
exit(0);
}
static void
mcenabletrace(void)
{
if(mctraceenabled)
mcdotrace = 1;
else
{
fprintf(stderr,
"Error: trace not enabled.\n"
"Please re-run the McStas compiler"
"with the --trace option, or rerun the\n"
"C compiler with the MC_TRACE_ENABLED macro defined.\n");
exit(1);
}
}
static void
mcuse_dir(char *dir)
{
#ifdef MC_PORTABLE
fprintf(stderr, "Error: "
"Directory output cannot be used with portable simulation.\n");
exit(1);
#else /* !MC_PORTABLE */
if(mkdir(dir, 0777))
{
fprintf(stderr, "Error: unable to create directory '%s'.\n", dir);
fprintf(stderr, "(Maybe the directory already exists?)\n");
exit(1);
}
mcdirname = dir;
#endif /* !MC_PORTABLE */
}
static void
mcuse_file(char *file)
{
mcsiminfo_name = file;
mcsingle_file = 1;
}
void
mcparseoptions(int argc, char *argv[])
{
int i, j, pos;
char *p;
int paramset = 0, *paramsetarray;
/* Add one to mcnumipar to aviod allocating zero size memory block. */
paramsetarray = malloc((mcnumipar + 1)*sizeof(*paramsetarray));
if(paramsetarray == NULL)
{
fprintf(stderr, "Error: insufficient memory\n");
exit(1);
}
for(j = 0; j < mcnumipar; j++)
paramsetarray[j] = 0;
for(i = 1; i < argc; i++)
{
if(!strcmp("-s", argv[i]) && (i + 1) < argc)
mcsetseed(argv[++i]);
else if(!strncmp("-s", argv[i], 2))
mcsetseed(&argv[i][2]);
else if(!strcmp("--seed", argv[i]) && (i + 1) < argc)
mcsetseed(argv[++i]);
else if(!strncmp("--seed=", argv[i], 7))
mcsetseed(&argv[i][7]);
else if(!strcmp("-n", argv[i]) && (i + 1) < argc)
mcsetn_arg(argv[++i]);
else if(!strncmp("-n", argv[i], 2))
mcsetn_arg(&argv[i][2]);
else if(!strcmp("--ncount", argv[i]) && (i + 1) < argc)
mcsetn_arg(argv[++i]);
else if(!strncmp("--ncount=", argv[i], 9))
mcsetn_arg(&argv[i][9]);
else if(!strcmp("-d", argv[i]) && (i + 1) < argc)
mcuse_dir(argv[++i]);
else if(!strncmp("-d", argv[i], 2))
mcuse_dir(&argv[i][2]);
else if(!strcmp("--dir", argv[i]) && (i + 1) < argc)
mcuse_dir(argv[++i]);
else if(!strncmp("--dir=", argv[i], 6))
mcuse_dir(&argv[i][6]);
else if(!strcmp("-f", argv[i]) && (i + 1) < argc)
mcuse_file(argv[++i]);
else if(!strncmp("-f", argv[i], 2))
mcuse_file(&argv[i][2]);
else if(!strcmp("--file", argv[i]) && (i + 1) < argc)
mcuse_file(argv[++i]);
else if(!strncmp("--file=", argv[i], 7))
mcuse_file(&argv[i][7]);
else if(!strcmp("-h", argv[i]))
mcshowhelp(argv[0]);
else if(!strcmp("--help", argv[i]))
mcshowhelp(argv[0]);
else if(!strcmp("-i", argv[i]))
mcinfo();
else if(!strcmp("--info", argv[i]))
mcinfo();
else if(!strcmp("-t", argv[i]))
mcenabletrace();
else if(!strcmp("--trace", argv[i]))
mcenabletrace();
else if(!strcmp("-a", argv[i]))
mcascii_only = 1;
else if(!strcmp("--ascii-only", argv[i]))
mcascii_only = 1;
else if(!strcmp("--no-output-files", argv[i]))
mcdisable_output_files = 1;
else if(argv[i][0] != '-' && (p = strchr(argv[i], '=')) != NULL)
{
*p++ = '\0';
for(j = 0; j < mcnumipar; j++)
if(!strcmp(mcinputtable[j].name, argv[i]))
{
int status;
status = (*mcinputtypes[mcinputtable[j].type].getparm)(p,
mcinputtable[j].par);
if(!status)
{
(*mcinputtypes[mcinputtable[j].type].error)
(mcinputtable[j].name, p);
exit(1);
}
paramsetarray[j] = 1;
paramset = 1;
break;
}
if(j == mcnumipar)
{ /* Unrecognized parameter name */
fprintf(stderr, "Error: unrecognized parameter %s\n", argv[i]);
exit(1);
}
}
else
mcusage(argv[0]);
}
if(!paramset)
mcreadparams(); /* Prompt for parameters if not specified. */
else
{
for(j = 0; j < mcnumipar; j++)
if(!paramsetarray[j])
{
fprintf(stderr, "Error: Instrument parameter %s left unset\n",
mcinputtable[j].name);
exit(1);
}
}
free(paramsetarray);
}
#ifndef MAC
#ifndef WIN32
/* This is the signal handler that makes simulation stop, and save results */
void sighandler(int sig)
{
char *tmp;
char *atfile;
FILE *fnum;
time_t t1;
t1 = time(NULL);
printf("\n# McStas: Signal %i detected in %s (%s): ", sig, mcinstrument_name, mcinstrument_source);
switch (sig) {
case SIGINT : printf(" SIGINT "); break; /* Ctrl-C */
case SIGQUIT : printf(" SIGQUIT "); break;
case SIGABRT : printf(" SIGABRT "); break;
case SIGTRAP : printf(" SIGTRAP "); break;
case SIGTERM : printf(" SIGTERM "); break;
case SIGPIPE : printf(" SIGPIPE "); break;
case SIGPWR : printf(" SIGPWR "); break;
case SIGUSR1 : printf(" SIGUSR1 "); break;
case SIGUSR2 : printf(" SIGUSR2 "); break;
case SIGILL : printf(" SIGILL "); break;
case SIGFPE : printf(" SIGFPE "); break;
case SIGBUS : printf(" SIGBUS "); break;
case SIGSEGV : printf(" SIGSEGV "); break;
case SIGURG : printf(" SIGURG "); break;
default : break;
}
printf("\n");
printf("# Date %s",ctime(&t1));
printf("# McStas: simulation now at %.2f %% (%i/%i)\n", 100*(double)(mcget_run_num())/mcget_ncount(), mcget_run_num(), mcget_ncount());
fflush(stdout);
if (sig == SIGUSR1)
{
return;
}
else
if ((sig == SIGUSR2) || (sig == SIGQUIT) || (sig == SIGTERM) || (sig == SIGABRT) || (sig == SIGALRM) || (sig == SIGINT))
{
printf("# McStas: finishing simulation\n");
mcset_ncount(mcget_run_num());
return;
}
else
{
printf("McStas: SYSTEM stop\n");
exit(-1);
}
}
#endif /* !MAC */
#endif /* !WIN32 */
/* McStas main() function. */
int
mcstas_main(int argc, char *argv[])
{
/* int run_num = 0; */
#ifdef MAC
argc = ccommand(&argv);
#endif
#ifndef MAC
#ifndef WIN32
/* install sig handler, but only once !! */
signal( SIGINT ,sighandler); /* interrupt (rubout) */
signal( SIGQUIT ,sighandler); /* quit (ASCII FS) */
signal( SIGABRT ,sighandler); /* used by abort, replace SIGIOT in the future */
signal( SIGTRAP ,sighandler); /* trace trap (not reset when caught) */
signal( SIGTERM ,sighandler); /* software termination signal from kill */
signal( SIGPIPE ,sighandler); /* write on a pipe with no one to read it */
signal( SIGPWR ,sighandler);
signal( SIGUSR1 ,sighandler); /* display simulation status */
signal( SIGUSR2 ,sighandler);
signal( SIGILL ,sighandler); /* illegal instruction (not reset when caught) */
signal( SIGFPE ,sighandler); /* floating point exception */
signal( SIGBUS ,sighandler); /* bus error */
signal( SIGSEGV ,sighandler); /* segmentation violation */
#ifdef SIGSYS
signal( SIGSYS ,sighandler); /* bad argument to system call */
#endif
signal( SIGURG ,sighandler); /* urgent socket condition */
#endif /* !MAC */
#endif /* !WIN32 */
srandom(time(NULL));
mcparseoptions(argc, argv);
mcsiminfo_init();
mcinit();
while(mcrun_num < mcncount)
{
mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1);
mcraytrace();
mcrun_num++;
}
mcfinally();
mcsiminfo_close();
return 0;
}
-------------- 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.13.9
* 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 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, 2tht Feb 2001 : 3He gas absorption handling (0.13.8)
* Modified by: EF, 5tht Mar 2001 : intermediate savings (0.13.9)
*
* 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.
* 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 (deg) vertical divergence (y)
* hdiv divergence (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)
* cylinder To get a cylindrical monitor (e.g. banana)
* (radius along x, height along y).
* disk Disk flat xy monitor
* radius is max abs value of xmin xmax ymin ymax
* file=string Detector image file name. default is component
* name, plus date and variable extension.
* square Square flat xy monitor
* limits=[min max] Lower/Upper limits for axes
* (see below 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
* sphrere To get a spherical monitor (e.g. a 4PI)
* radius is max abs value of xmin xmax ymin ymax
* 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
*
* <b>EXAMPLES:</b>
* MyMon = Monitor_nD(
* xmin = -0.1, xmax = 0.1,
* ymin = -0.1, ymax = 0.1,
* 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(
* xmin = -0.05, xmax = 0.05,
* ymin = -0.1, ymax = 0.1,
* 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:
*
* 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
* options:(str) String that specifies the configuration of the monitor</br>
* The general syntax is "[x] options..." (see <b>Descr.</b>)
*
* OUTPUT PARAMETERS:
*
* Nsum: Number of neutron hits (counts)
* psum: Sum of neutron weights (flux)
* p2sum: 2nd moment of neutron weights (error)
* Mon2D_N: output array containing detected signals counts (counts)
* Mon2D_p: (flux)
* Mon2D_p2: (error)
* Mon2D_Buffer: used for long lists and auto limits
* 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 (xmin, xmax, ymin, ymax)
/* these are protected C variables */
OUTPUT PARAMETERS (Nsum, psum, p2sum, Mon2D_N, Mon2D_p, Mon2D_p2, Mon2D_Buffer, 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_Here
#define MONnD_COORD_NMAX 25 /* 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;
} 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 */
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];
} MonitornD_Variables_type;
/* this routine is used to save results at end of simulation, but also
* during simulation (SIGUSR... intermediate steps) */
void Monitor_nD_OutPut(aNsum, apsum, ap2sum, aMon2D_N, aMon2D_p, aMon2D_p2, aMon2D_Buffer, aDEFS, aVars, dofree, amccompcurname)
int aNsum;
double apsum, ap2sum;
int **aMon2D_N;
double **aMon2D_p;
double **aMon2D_p2;
double *aMon2D_Buffer;
MonitornD_Defines_type aDEFS;
MonitornD_Variables_type aVars;
int dofree;
char *amccompcurname;
/* dofree = 0 : no free, =1 : free variables (last call) */
#ifdef mccompcurname
#define mccompcurname_sav mccompcurname
#undef mccompcurname
#define mccompcurname amccompcurname
#endif
{
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", amccompcurname, 100*(double)(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", amccompcurname, 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 = aMon2D_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) */
/* 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] = aMon2D_Buffer[While_Buffer*(aVars.Coord_Number+2) + (i-1)];
Coord[0] = aMon2D_Buffer[While_Buffer*(aVars.Coord_Number+2) + (aVars.Coord_Number)];
pp = Coord[0];
if (XY > 0) Coord_Index[i] = floor((aMon2D_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(aMon2D_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])
{
aMon2D_N[i][j]++;
aMon2D_p[i][j] += pp;
aMon2D_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])
{
aMon2D_N[i][j]++;
aMon2D_p[i][j] += pp;
aMon2D_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", amccompcurname, aVars.Monitor_Label);
printf("Monitor_nD: with options=%s\n", 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, aNsum, apsum, ap2sum);
}
else
if (strlen(aVars.Mon_File) > 0)
{
fname = (char*)malloc(strlen(aVars.Mon_File)+10*aVars.Coord_Number);
if (aVars.Flag_List) /* List */
{
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", amccompcurname, 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*(double)(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 transpose. Skipping.\n", amccompcurname);
if (p0m != NULL) free(p0m);
DETECTOR_OUT_2D(
label,
aVars.Coord_Label[0],
Coord_X_Label,
min2d, max2d,
min1d, max1d,
bin2d,
bin1d,
(int *)aMon2D_Buffer,aMon2D_Buffer,aMon2D_Buffer,
fname);
}
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] = aMon2D_Buffer[j*(aVars.Coord_Number+2) + i];
p0m[i*aVars.Buffer_Size+j] = 1;
}
DETECTOR_OUT_2D(
label,
Coord_X_Label,
aVars.Coord_Label[0],
min1d, max1d,
min2d, max2d,
bin1d,
bin2d,
p0m,p1m,aMon2D_Buffer,
fname);
free(p0m);
free(p1m);
}
}
if (aVars.Flag_Multiple) /* n1D */
{
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*(double)(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", amccompcurname, 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;
DETECTOR_OUT_1D(
label,
aVars.Coord_Label[i+1],
aVars.Coord_Label[0],
aVars.Coord_Var[i+1],
min1d, max1d,
aVars.Coord_Bin[i+1],
aMon2D_N[i],aMon2D_p[i],aMon2D_p2[i],
fname);
}
}
else
if (aVars.Coord_Number == 2) /* 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", amccompcurname, 3*aVars.Coord_Bin[1]*aVars.Coord_Bin[2]*sizeof(double));
if (p0m != NULL) free(p0m);
if (p1m != NULL) free(p1m);
}
else
{
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]] = aMon2D_N[i][j];
p1m[j + i*aVars.Coord_Bin[2]] = aMon2D_p[i][j];
p2m[j + i*aVars.Coord_Bin[2]] = aMon2D_p2[i][j];
}
}
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", amccompcurname, 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*(double)(mcget_run_num())/mcget_ncount());
}
else
strcpy(label, aVars.Monitor_Label);
DETECTOR_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);
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 (aMon2D_Buffer != NULL && dofree) free(aMon2D_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(aMon2D_N[i]);
free(aMon2D_p[i]);
free(aMon2D_p2[i]);
}
free(aMon2D_N);
free(aMon2D_p);
free(aMon2D_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(aMon2D_N[i]);
free(aMon2D_p[i]);
free(aMon2D_p2[i]);
}
free(aMon2D_N);
free(aMon2D_p);
free(aMon2D_p2);
}
}
#ifdef mccompcurname_sav
#undef mccompcurname
#define mccompcurname mccompcurname_sav
#undef mccompcurname_sav
#endif
#define Monitor_nD_Here
#endif
int Nsum;
double psum, p2sum;
int **Mon2D_N;
double **Mon2D_p;
double **Mon2D_p2;
double *Mon2D_Buffer;
MonitornD_Defines_type DEFS;
MonitornD_Variables_type Vars;
%}
INITIALIZE
%{
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_Vars_Coord_Type;
char Set_Coord_Flag = 0;
char Set_Vars_Coord_Label[30];
char Set_Vars_Coord_Var[30];
char Short_Label[MONnD_COORD_NMAX][30];
char Set_Coord_Mode;
long i=0;
double lmin, lmax;
long t;
t = (long)time(NULL);
DEFS.COORD_NONE =0;
DEFS.COORD_X =1;
DEFS.COORD_Y =2;
DEFS.COORD_Z =3;
DEFS.COORD_VX =4;
DEFS.COORD_VY =5;
DEFS.COORD_VZ =6;
DEFS.COORD_T =7;
DEFS.COORD_P =8;
DEFS.COORD_SX =9;
DEFS.COORD_SY =10;
DEFS.COORD_SZ =11;
DEFS.COORD_KX =12;
DEFS.COORD_KY =13;
DEFS.COORD_KZ =14;
DEFS.COORD_K =15;
DEFS.COORD_V =16;
DEFS.COORD_ENERGY =17;
DEFS.COORD_LAMBDA =18;
DEFS.COORD_RADIUS =19;
DEFS.COORD_HDIV =20;
DEFS.COORD_VDIV =21;
DEFS.COORD_ANGLE =22;
DEFS.COORD_NCOUNT =23;
DEFS.COORD_THETA =24;
DEFS.COORD_PHI =25;
DEFS.COORD_USER1 =26;
DEFS.COORD_USER2 =26;
/* token modifiers */
DEFS.COORD_VAR =0; /* next token should be a variable or normal option */
DEFS.COORD_MIN =1; /* next token is a min value */
DEFS.COORD_MAX =2; /* next token is a max value */
DEFS.COORD_DIM =3; /* next token is a bin value */
DEFS.COORD_FIL =4; /* next token is a filename */
DEFS.COORD_EVNT =5; /* next token is a buffer size value */
DEFS.COORD_3HE =6; /* next token is a 3He pressure value */
DEFS.COORD_INTERM =7; /* next token is an intermediate save value (%) */
strcpy(DEFS.TOKEN_DEL, " =,;[](){}:"); /* token separators */
DEFS.SHAPE_SQUARE =0; /* shape of the monitor */
DEFS.SHAPE_DISK =1;
DEFS.SHAPE_SPHERE =2;
DEFS.SHAPE_CYLIND =3;
Vars.Sphere_Radius = 0;
Vars.Cylinder_Height = 0;
Vars.Flag_With_Borders = 0; /* 2 means xy borders too */
Vars.Flag_List = 0; /* 1 store 1 buffer, 2 is list all */
Vars.Flag_Multiple = 0; /* 1 when n1D, 0 for 2D */
Vars.Flag_Verbose = 0;
Vars.Flag_Shape = DEFS.SHAPE_SQUARE;
Vars.Flag_Auto_Limits = 0; /* get limits from first Buffer */
Vars.Flag_Absorb = 0; /* monitor is also a slit */
Vars.Flag_per_cm2 = 0; /* flux is per cm2 */
Vars.Coord_Number = 0; /* total number of variables to monitor, plus intensity (0) */
Vars.Buffer_Block = 1000; /* Buffer size for list or auto limits */
Vars.Neutron_Counter = 0; /* event counter, simulation total counts is mcget_ncount() */
Vars.Buffer_Counter = 0; /* index in Buffer size (for realloc) */
Vars.Buffer_Size = 0;
Vars.UserVariable1 = 0;
Vars.UserVariable2 = 0;
Vars.He3_pressure = 0;
Vars.IntermediateCnts = 0;
Set_Vars_Coord_Type = DEFS.COORD_NONE;
Set_Coord_Mode = DEFS.COORD_VAR;
/* parse option string */
strcpy(Vars.option, options);
option_copy = (char*)malloc(strlen(options));
if (option_copy == NULL)
{
printf("Monitor_nD: %s cannot allocate option_copy (%i). Fatal.\n", mccompcurname, strlen(options));
exit(-1);
}
if (strlen(options))
{
Flag_End = 0;
strcpy(option_copy, options);
}
if (strstr(options, "cm2") || strstr(options, "cm^2")) Vars.Flag_per_cm2 = 1;
if (Vars.Flag_per_cm2) strcpy(Vars.Coord_Label[0],"Intensity [n/cm^2/s]");
else strcpy(Vars.Coord_Label[0],"Intensity [n/s]");
strcpy(Vars.Coord_Var[0],"p");
Vars.Coord_Type[0] = DEFS.COORD_P;
Vars.Coord_Bin[0] = 1;
Vars.Coord_Min[0] = 0;
Vars.Coord_Max[0] = FLT_MAX;
/* default file name is comp name+dateID */
sprintf(Vars.Mon_File, "%s_%i", mccompcurname, 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,DEFS.TOKEN_DEL);
else token=(char *)strtok(NULL,DEFS.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 == DEFS.COORD_MAX)
{
if (!Flag_All)
Vars.Coord_Max[Vars.Coord_Number] = atof(token);
else
for (i = 0; i <= Vars.Coord_Number; Vars.Coord_Max[i++] = atof(token));
Set_Coord_Mode = DEFS.COORD_VAR; Flag_All = 0;
}
if (Set_Coord_Mode == DEFS.COORD_MIN)
{
if (!Flag_All)
Vars.Coord_Min[Vars.Coord_Number] = atof(token);
else
for (i = 0; i <= Vars.Coord_Number; Vars.Coord_Min[i++] = atof(token));
Set_Coord_Mode = DEFS.COORD_MAX;
}
if (Set_Coord_Mode == DEFS.COORD_DIM)
{
if (!Flag_All)
Vars.Coord_Bin[Vars.Coord_Number] = atoi(token);
else
for (i = 0; i <= Vars.Coord_Number; Vars.Coord_Bin[i++] = atoi(token));
Set_Coord_Mode = DEFS.COORD_VAR; Flag_All = 0;
}
if (Set_Coord_Mode == DEFS.COORD_FIL)
{
if (!Flag_No) strcpy(Vars.Mon_File,token);
else { strcpy(Vars.Mon_File,""); Vars.Coord_Number = 0; Flag_End = 1;}
Set_Coord_Mode = DEFS.COORD_VAR;
}
if (Set_Coord_Mode == DEFS.COORD_EVNT)
{
if (!strcmp(token, "all") || Flag_All) Vars.Flag_List = 2;
else { i = atoi(token); if (i) Vars.Buffer_Counter = i; }
Set_Coord_Mode = DEFS.COORD_VAR; Flag_All = 0;
}
if (Set_Coord_Mode == DEFS.COORD_3HE)
{
Vars.He3_pressure = atof(token);
Set_Coord_Mode = DEFS.COORD_VAR; Flag_All = 0;
}
if (Set_Coord_Mode == DEFS.COORD_INTERM)
{
Vars.Intermediate = atof(token);
Set_Coord_Mode = DEFS.COORD_VAR; Flag_All = 0;
}
/* now look for general option keywords */
if (!strcmp(token, "borders"))
{ if (Flag_No) { Vars.Flag_With_Borders = 0; Flag_No = 0; }
else Vars.Flag_With_Borders = 1; }
if (!strcmp(token, "verbose"))
{ if (Flag_No) { Vars.Flag_Verbose = 0; Flag_No = 0; }
else Vars.Flag_Verbose = 1; }
if (!strcmp(token, "multiple"))
{ if (Flag_No) { Vars.Flag_Multiple = 0; Flag_No = 0; }
else Vars.Flag_Multiple = 1; }
if (!strcmp(token, "list"))
{ if (Flag_No) { Vars.Flag_List = 0; Flag_No = 0; }
else Vars.Flag_List = 1;
Set_Coord_Mode = DEFS.COORD_EVNT; }
if (!strcmp(token, "limits") || !strcmp(token, "min")) Set_Coord_Mode = DEFS.COORD_MIN;
if (!strcmp(token, "slit") || !strcmp(token, "absorb"))
{ if (Flag_No) { Vars.Flag_Absorb = 0; Flag_No = 0; }
else Vars.Flag_Absorb = 1; }
if (!strcmp(token, "max")) Set_Coord_Mode = DEFS.COORD_MAX;
if (!strcmp(token, "bins")) Set_Coord_Mode = DEFS.COORD_DIM;
if (!strcmp(token, "file"))
{ Set_Coord_Mode = DEFS.COORD_FIL;
if (Flag_No) { strcpy(Vars.Mon_File,""); Vars.Coord_Number = 0; Flag_End = 1;}}
if (!strcmp(token, "unactivate")) { Flag_End = 1; Vars.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) { Vars.Flag_Shape = DEFS.SHAPE_SQUARE; Flag_No = 0; }
else Vars.Flag_Shape = DEFS.SHAPE_SPHERE; }
if (!strcmp(token, "cylinder"))
{ if (Flag_No) { Vars.Flag_Shape = DEFS.SHAPE_SQUARE; Flag_No = 0; }
else Vars.Flag_Shape = DEFS.SHAPE_CYLIND; }
if (!strcmp(token, "square")) Vars.Flag_Shape = DEFS.SHAPE_SQUARE;
if (!strcmp(token, "disk")) Vars.Flag_Shape = DEFS.SHAPE_DISK;
if (!strcmp(token, "auto"))
{ if (Flag_No) { Vars.Flag_Auto_Limits = 0; Flag_No = 0; }
else Vars.Flag_Auto_Limits = 1; }
if (!strcmp(token, "premonitor"))
{ if (Flag_No) { Vars.Flag_UsePreMonitor = 0; Flag_No = 0; }
else Vars.Flag_UsePreMonitor == 1; }
if (!strcmp(token, "3He_pressure"))
{ if (!Flag_No) Set_Coord_Mode = DEFS.COORD_3HE;
Vars.He3_pressure = 3; }
if (!strcmp(token, "intermediate"))
{ if (!Flag_No) Set_Coord_Mode = DEFS.COORD_INTERM;
Vars.Intermediate = 5; }
if (!strcmp(token, "no") || !strcmp(token, "not")) Flag_No = 1;
/* now look for variable names to monitor */
Set_Vars_Coord_Type = DEFS.COORD_NONE; lmin = 0; lmax = 0;
if (!strcmp(token, "x"))
{ Set_Vars_Coord_Type = DEFS.COORD_X; strcpy(Set_Vars_Coord_Label,"x [m]"); strcpy(Set_Vars_Coord_Var,"x"); lmin = xmin; lmax = xmax; }
if (!strcmp(token, "y"))
{ Set_Vars_Coord_Type = DEFS.COORD_Y; strcpy(Set_Vars_Coord_Label,"y [m]"); strcpy(Set_Vars_Coord_Var,"y"); lmin = ymin; lmax = ymax; }
if (!strcmp(token, "z"))
{ Set_Vars_Coord_Type = DEFS.COORD_Z; strcpy(Set_Vars_Coord_Label,"z [m]"); strcpy(Set_Vars_Coord_Var,"z"); lmin = 0; lmax = 100; }
if (!strcmp(token, "k") || !strcmp(token, "wavevector"))
{ Set_Vars_Coord_Type = DEFS.COORD_K; strcpy(Set_Vars_Coord_Label,"|k| [Angs-1]"); strcpy(Set_Vars_Coord_Var,"k"); lmin = 0; lmax = 10; }
if (!strcmp(token, "v"))
{ Set_Vars_Coord_Type = DEFS.COORD_V; strcpy(Set_Vars_Coord_Label,"Velocity [m/s]"); strcpy(Set_Vars_Coord_Var,"v"); lmin = 0; lmax = 10000; }
if (!strcmp(token, "t") || !strcmp(token, "time"))
{ Set_Vars_Coord_Type = DEFS.COORD_T; strcpy(Set_Vars_Coord_Label,"TOF [s]"); strcpy(Set_Vars_Coord_Var,"t"); lmin = 0; lmax = .1; }
if ((Vars.Coord_Number > 0) && (!strcmp(token, "p") || !strcmp(token, "intensity") || !strcmp(token, "flux")))
{ Set_Vars_Coord_Type = DEFS.COORD_P;
if (Vars.Flag_per_cm2) strcpy(Set_Vars_Coord_Label,"Intensity [n/cm^2/s]");
else strcpy(Set_Vars_Coord_Label,"Intensity [n/s]");
strcpy(Set_Vars_Coord_Var,"I");
lmin = 0; lmax = FLT_MAX; }
if (!strcmp(token, "vx"))
{ Set_Vars_Coord_Type = DEFS.COORD_VX; strcpy(Set_Vars_Coord_Label,"vx [m/s]"); strcpy(Set_Vars_Coord_Var,"vx"); lmin = -1000; lmax = 1000; }
if (!strcmp(token, "vy"))
{ Set_Vars_Coord_Type = DEFS.COORD_VY; strcpy(Set_Vars_Coord_Label,"vy [m/s]"); strcpy(Set_Vars_Coord_Var,"vy"); lmin = -1000; lmax = 1000; }
if (!strcmp(token, "vz"))
{ Set_Vars_Coord_Type = DEFS.COORD_VZ; strcpy(Set_Vars_Coord_Label,"vz [m/s]"); strcpy(Set_Vars_Coord_Var,"vz"); lmin = -10000; lmax = 10000; }
if (!strcmp(token, "kx"))
{ Set_Vars_Coord_Type = DEFS.COORD_KX; strcpy(Set_Vars_Coord_Label,"kx [Angs-1]"); strcpy(Set_Vars_Coord_Var,"kx"); lmin = -1; lmax = 1; }
if (!strcmp(token, "ky"))
{ Set_Vars_Coord_Type = DEFS.COORD_KY; strcpy(Set_Vars_Coord_Label,"ky [Angs-1]"); strcpy(Set_Vars_Coord_Var,"ky"); lmin = -1; lmax = 1; }
if (!strcmp(token, "kz"))
{ Set_Vars_Coord_Type = DEFS.COORD_KZ; strcpy(Set_Vars_Coord_Label,"kz [Angs-1]"); strcpy(Set_Vars_Coord_Var,"kz"); lmin = -10; lmax = 10; }
if (!strcmp(token, "sx"))
{ Set_Vars_Coord_Type = DEFS.COORD_SX; strcpy(Set_Vars_Coord_Label,"sx [1]"); strcpy(Set_Vars_Coord_Var,"sx"); lmin = -1; lmax = 1; }
if (!strcmp(token, "sy"))
{ Set_Vars_Coord_Type = DEFS.COORD_SY; strcpy(Set_Vars_Coord_Label,"sy [1]"); strcpy(Set_Vars_Coord_Var,"sy"); lmin = -1; lmax = 1; }
if (!strcmp(token, "sz"))
{ Set_Vars_Coord_Type = DEFS.COORD_SZ; strcpy(Set_Vars_Coord_Label,"sz [1]"); strcpy(Set_Vars_Coord_Var,"sz"); lmin = -1; lmax = 1; }
if (!strcmp(token, "energy") || !strcmp(token, "omega"))
{ Set_Vars_Coord_Type = DEFS.COORD_ENERGY; strcpy(Set_Vars_Coord_Label,"Energy [meV]"); strcpy(Set_Vars_Coord_Var,"E"); lmin = 0; lmax = 100; }
if (!strcmp(token, "lambda") || !strcmp(token, "wavelength"))
{ Set_Vars_Coord_Type = DEFS.COORD_LAMBDA; strcpy(Set_Vars_Coord_Label,"Wavelength [Angs]"); strcpy(Set_Vars_Coord_Var,"L"); lmin = 0; lmax = 100; }
if (!strcmp(token, "radius"))
{ Set_Vars_Coord_Type = DEFS.COORD_RADIUS; strcpy(Set_Vars_Coord_Label,"Radius [m]"); strcpy(Set_Vars_Coord_Var,"R"); lmin = 0; lmax = xmax; }
if (!strcmp(token, "angle"))
{ Set_Vars_Coord_Type = DEFS.COORD_ANGLE; strcpy(Set_Vars_Coord_Label,"Angle [deg]"); strcpy(Set_Vars_Coord_Var,"A"); lmin = -5; lmax = 5; }
if (!strcmp(token, "hdiv")|| !strcmp(token, "divergence"))
{ Set_Vars_Coord_Type = DEFS.COORD_HDIV; strcpy(Set_Vars_Coord_Label,"Hor. Divergence [deg]"); strcpy(Set_Vars_Coord_Var,"HD"); lmin = -5; lmax = 5; }
if (!strcmp(token, "vdiv"))
{ Set_Vars_Coord_Type = DEFS.COORD_VDIV; strcpy(Set_Vars_Coord_Label,"Vert. Divergence [deg]"); strcpy(Set_Vars_Coord_Var,"VD"); lmin = -5; lmax = 5; }
if (!strcmp(token, "theta") || !strcmp(token, "longitude"))
{ Set_Vars_Coord_Type = DEFS.COORD_THETA; strcpy(Set_Vars_Coord_Label,"Longitude [deg]"); strcpy(Set_Vars_Coord_Var,"th"); lmin = -180; lmax = 180; }
if (!strcmp(token, "phi") || !strcmp(token, "lattitude"))
{ Set_Vars_Coord_Type = DEFS.COORD_PHI; strcpy(Set_Vars_Coord_Label,"Lattitude [deg]"); strcpy(Set_Vars_Coord_Var,"ph"); lmin = -180; lmax = 180; }
if (!strcmp(token, "ncounts"))
{ Set_Vars_Coord_Type = DEFS.COORD_NCOUNT; strcpy(Set_Vars_Coord_Label,"Neutrons [1]"); strcpy(Set_Vars_Coord_Var,"N"); lmin = 0; lmax = 1e10; }
if (!strcmp(token, "user") || !strcmp(token, "user1"))
{ Set_Vars_Coord_Type = DEFS.COORD_USER1; strncpy(Set_Vars_Coord_Label,Vars.UserName1,32); strcpy(Set_Vars_Coord_Var,"U1"); lmin = -1e10; lmax = 1e10; }
if (!strcmp(token, "user2"))
{ Set_Vars_Coord_Type = DEFS.COORD_USER2; strncpy(Set_Vars_Coord_Label,Vars.UserName2,32); strcpy(Set_Vars_Coord_Var,"U2"); lmin = -1e10; lmax = 1e10; }
/* now stores variable keywords detected, if any */
if (Set_Vars_Coord_Type != DEFS.COORD_NONE)
{
if (Vars.Coord_Number < MONnD_COORD_NMAX) Vars.Coord_Number++;
else if (Vars.Flag_Verbose) printf("Monitor_nD: %s reached max number of variables (%i).\n", mccompcurname, MONnD_COORD_NMAX);
Vars.Coord_Type[Vars.Coord_Number] = Set_Vars_Coord_Type;
strcpy(Vars.Coord_Label[Vars.Coord_Number], Set_Vars_Coord_Label);
strcpy(Vars.Coord_Var[Vars.Coord_Number], Set_Vars_Coord_Var);
Vars.Coord_Min[Vars.Coord_Number] = lmin;
Vars.Coord_Max[Vars.Coord_Number] = lmax;
Vars.Coord_Bin[Vars.Coord_Number] = 20;
Set_Coord_Mode = DEFS.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", mccompcurname, 128);
if (strstr(options,"unactivate") && Vars.Flag_Verbose) printf("Monitor_nD: %s is unactivated (0D)\n", mccompcurname);
/* now setting Monitor Name from variable labels */
strcpy(Vars.Monitor_Label,"");
for (i = 0; i <= Vars.Coord_Number; i++)
{
Set_Vars_Coord_Type = Vars.Coord_Type[i];
if ((Set_Vars_Coord_Type == DEFS.COORD_THETA)
|| (Set_Vars_Coord_Type == DEFS.COORD_PHI)
|| (Set_Vars_Coord_Type == DEFS.COORD_X)
|| (Set_Vars_Coord_Type == DEFS.COORD_Y)
|| (Set_Vars_Coord_Type == DEFS.COORD_Z)
|| (Set_Vars_Coord_Type == DEFS.COORD_RADIUS))
strcpy(Short_Label[i],"Position");
else
if ((Set_Vars_Coord_Type == DEFS.COORD_VX)
|| (Set_Vars_Coord_Type == DEFS.COORD_VY)
|| (Set_Vars_Coord_Type == DEFS.COORD_VZ)
|| (Set_Vars_Coord_Type == DEFS.COORD_V))
strcpy(Short_Label[i],"Velocity");
else
if ((Set_Vars_Coord_Type == DEFS.COORD_KX)
|| (Set_Vars_Coord_Type == DEFS.COORD_KY)
|| (Set_Vars_Coord_Type == DEFS.COORD_KZ)
|| (Set_Vars_Coord_Type == DEFS.COORD_K))
strcpy(Short_Label[i],"Wavevector");
else
if ((Set_Vars_Coord_Type == DEFS.COORD_SX)
|| (Set_Vars_Coord_Type == DEFS.COORD_SY)
|| (Set_Vars_Coord_Type == DEFS.COORD_SZ))
strcpy(Short_Label[i],"Spin");
else
if ((Set_Vars_Coord_Type == DEFS.COORD_HDIV)
|| (Set_Vars_Coord_Type == DEFS.COORD_VDIV)
|| (Set_Vars_Coord_Type == DEFS.COORD_ANGLE))
strcpy(Short_Label[i],"Divergence");
else
if (Set_Vars_Coord_Type == DEFS.COORD_ENERGY)
strcpy(Short_Label[i],"Energy");
else
if (Set_Vars_Coord_Type == DEFS.COORD_LAMBDA)
strcpy(Short_Label[i],"Wavelength");
else
if (Set_Vars_Coord_Type == DEFS.COORD_NCOUNT)
strcpy(Short_Label[i],"Neutron counts");
else
if (Set_Vars_Coord_Type == DEFS.COORD_T)
strcpy(Short_Label[i],"Time Of Flight");
else
if (Set_Vars_Coord_Type == DEFS.COORD_P)
strcpy(Short_Label[i],"Intensity");
else
if (Set_Vars_Coord_Type == DEFS.COORD_USER1)
strncpy(Short_Label[i],Vars.UserName1,32);
else
if (Set_Vars_Coord_Type == DEFS.COORD_USER2)
strncpy(Short_Label[i],Vars.UserName2,32);
else
strcpy(Short_Label[i],"Unknown");
strcat(Vars.Monitor_Label, " ");
strcat(Vars.Monitor_Label, Short_Label[i]);
} /* end for Short_Label */
strcat(Vars.Monitor_Label, " Monitor");
if (Vars.Flag_Shape == DEFS.SHAPE_SQUARE) strcat(Vars.Monitor_Label, " (Square)");
if (Vars.Flag_Shape == DEFS.SHAPE_DISK) strcat(Vars.Monitor_Label, " (Disk)");
if (Vars.Flag_Shape == DEFS.SHAPE_SPHERE) strcat(Vars.Monitor_Label, " (Sphere)");
if (Vars.Flag_Shape == DEFS.SHAPE_CYLIND) strcat(Vars.Monitor_Label, " (Cylinder)");
if (((Vars.Flag_Shape == DEFS.SHAPE_CYLIND) || (Vars.Flag_Shape == DEFS.SHAPE_SPHERE))
&& strstr(options, "outgoing"))
{
Vars.Flag_Shape *= -1;
strcat(Vars.Monitor_Label, " [out]");
}
if (Vars.Flag_UsePreMonitor == 1)
{
strcat(Vars.Monitor_Label, " at ");
strncat(Vars.Monitor_Label, Vars.UserName1,32);
}
/* Vars.Coord_Number 0 : intensity
* Vars.Coord_Number 1:n : detector variables */
/* now allocate memory to store variables in TRACE */
if ((Vars.Coord_Number != 2) && !Vars.Flag_Multiple && !Vars.Flag_List)
{ Vars.Flag_Multiple = 1; Vars.Flag_List = 0; } /* default is n1D */
/* list and auto limits case : Vars.Flag_List or Vars.Flag_Auto_Limits
* -> Buffer to flush and suppress after Vars.Flag_Auto_Limits
*/
if ((Vars.Flag_Auto_Limits || Vars.Flag_List) && Vars.Coord_Number)
{ /* Dim : (Vars.Coord_Number+2)*Vars.Buffer_Block matrix (for p, dp) */
Mon2D_Buffer = (double *)malloc((Vars.Coord_Number+2)*Vars.Buffer_Block*sizeof(double));
if (Mon2D_Buffer == NULL)
{ printf("Monitor_nD: %s cannot allocate Mon2D_Buffer (%li). No list and auto limits.\n", mccompcurname, Vars.Buffer_Block*(Vars.Coord_Number+2)*sizeof(double)); Vars.Flag_List = 0; Vars.Flag_Auto_Limits = 0; }
Vars.Buffer_Size = Vars.Buffer_Block;
}
/* 1D and n1D case : Vars.Flag_Multiple */
if (Vars.Flag_Multiple && Vars.Coord_Number)
{ /* Dim : Vars.Coord_Number*Vars.Coord_Bin[i] vectors */
Mon2D_N = (int **)malloc((Vars.Coord_Number)*sizeof(int *));
Mon2D_p = (double **)malloc((Vars.Coord_Number)*sizeof(double *));
Mon2D_p2 = (double **)malloc((Vars.Coord_Number)*sizeof(double *));
if ((Mon2D_N == NULL) || (Mon2D_p == NULL) || (Mon2D_p2 == NULL))
{ printf("Monitor_nD: %s n1D cannot allocate Mon2D_N/p/2p (%i). Fatal.\n", mccompcurname, (Vars.Coord_Number)*sizeof(double *)); exit(-1); }
for (i= 1; i <= Vars.Coord_Number; i++)
{
Mon2D_N[i-1] = (int *)malloc(Vars.Coord_Bin[i]*sizeof(int));
Mon2D_p[i-1] = (double *)malloc(Vars.Coord_Bin[i]*sizeof(double));
Mon2D_p2[i-1] = (double *)malloc(Vars.Coord_Bin[i]*sizeof(double));
if ((Mon2D_N == NULL) || (Mon2D_p == NULL) || (Mon2D_p2 == NULL))
{ printf("Monitor_nD: %s n1D cannot allocate %s Mon2D_N/p/2p[%li] (%i). Fatal.\n", mccompcurname, Vars.Coord_Var[i], i, (Vars.Coord_Bin[i])*sizeof(double *)); exit(-1); }
}
}
else /* 2D case : Vars.Coord_Number==2 and !Vars.Flag_Multiple and !Vars.Flag_List */
if ((Vars.Coord_Number == 2) && !Vars.Flag_Multiple)
{ /* Dim : Vars.Coord_Bin[1]*Vars.Coord_Bin[2] matrix */
Mon2D_N = (int **)malloc((Vars.Coord_Bin[1])*sizeof(int *));
Mon2D_p = (double **)malloc((Vars.Coord_Bin[1])*sizeof(double *));
Mon2D_p2 = (double **)malloc((Vars.Coord_Bin[1])*sizeof(double *));
if ((Mon2D_N == NULL) || (Mon2D_p == NULL) || (Mon2D_p2 == NULL))
{ printf("Monitor_nD: %s 2D cannot allocate %s Mon2D_N/p/2p (%i). Fatal.\n", mccompcurname, Vars.Coord_Var[1], (Vars.Coord_Bin[1])*sizeof(double *)); exit(-1); }
for (i= 0; i < Vars.Coord_Bin[1]; i++)
{
Mon2D_N[i] = (int *)malloc(Vars.Coord_Bin[2]*sizeof(int));
Mon2D_p[i] = (double *)malloc(Vars.Coord_Bin[2]*sizeof(double));
Mon2D_p2[i] = (double *)malloc(Vars.Coord_Bin[2]*sizeof(double));
if ((Mon2D_N == NULL) || (Mon2D_p == NULL) || (Mon2D_p2 == NULL))
{ printf("Monitor_nD: %s 2D cannot allocate %s Mon2D_N/p/2p[%li] (%i). Fatal.\n", mccompcurname, Vars.Coord_Var[1], i, (Vars.Coord_Bin[2])*sizeof(double *)); exit(-1); }
}
}
/* no Mon2D allocated for
* (Vars.Coord_Number != 2) && !Vars.Flag_Multiple && Vars.Flag_List */
psum = 0;
p2sum = 0;
Nsum = 0;
Vars.area = (xmax - xmin)*(ymax - ymin)*1E4; /* in cm**2 for square shapes */
Vars.Sphere_Radius = 0;
if (fabs(xmin) > Vars.Sphere_Radius) Vars.Sphere_Radius = fabs(xmin);
if (fabs(xmax) > Vars.Sphere_Radius) Vars.Sphere_Radius = fabs(xmax);
if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE))
{
if (fabs(ymin) > Vars.Sphere_Radius) Vars.Sphere_Radius = fabs(ymin);
if (fabs(ymax) > Vars.Sphere_Radius) Vars.Sphere_Radius = fabs(ymax);
Vars.area = PI*Vars.Sphere_Radius*Vars.Sphere_Radius; /* disk shapes */
}
Vars.Cylinder_Height = fabs(ymax-ymin);
if (Vars.Intermediate < 0) Vars.Intermediate = 0;
if (Vars.Intermediate > 1) Vars.Intermediate /= 100;
Vars.IntermediateCnts = Vars.Intermediate*mcget_ncount();
%}
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;
if (abs(Vars.Flag_Shape) == DEFS.SHAPE_SQUARE) /* square xy */
intersect = (x>=xmin && x<=xmax && y>=ymin && y<=ymax);
if (abs(Vars.Flag_Shape) == DEFS.SHAPE_DISK) /* disk xy */
intersect = ((x*x + y*y) <= Vars.Sphere_Radius*Vars.Sphere_Radius);
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); */
}
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);
}
if (intersect)
{
if ((abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND))
{
if (t0 < 0 && t1 > 0)
t0 = 0; /* neutron was already inside ! */
if (t1 < 0 && t0 > 0)
t1 = 0;
/* 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) && (abs(Vars.Flag_Shape) == DEFS.SHAPE_SPHERE) || (abs(Vars.Flag_Shape) == DEFS.SHAPE_CYLIND))
{
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;
}
pp = Vars.cp;
if (Vars.Flag_per_cm2 && Vars.area != 0) pp /= Vars.area;
Nsum++;
psum += pp;
p2sum += pp*pp;
if (Vars.Coord_Number > 0)
{
/* Vars.Flag_Auto_Limits */
if ((Vars.Buffer_Counter >= Vars.Buffer_Block) && (Vars.Flag_Auto_Limits == 1))
{
/* auto limits case : get limits in Buffer for each variable */
/* Dim : (Vars.Coord_Number+2)*Vars.Buffer_Block matrix (for p, dp) */
if (Vars.Flag_Verbose) printf("Monitor_nD: %s getting %i Auto Limits from List (%li).\n", mccompcurname, Vars.Coord_Number, Vars.Buffer_Counter);
for (i = 1; i <= Vars.Coord_Number; i++)
{
Vars.Coord_Min[i] = FLT_MAX;
Vars.Coord_Max[i] = -FLT_MAX;
for (j = 0; j < Vars.Buffer_Block; j++)
{
XY = Mon2D_Buffer[j*(Vars.Coord_Number+2) + (i-1)]; /* scanning variables in Buffer */
if (XY < Vars.Coord_Min[i]) Vars.Coord_Min[i] = XY;
if (XY > Vars.Coord_Max[i]) Vars.Coord_Max[i] = XY;
}
}
Vars.Flag_Auto_Limits = 2; /* pass to 2nd auto limits step */
}
/* manage realloc for list all if Buffer size exceeded */
if ((Vars.Buffer_Counter >= Vars.Buffer_Block) && (Vars.Flag_List == 2))
{
Mon2D_Buffer = (double *)realloc(Mon2D_Buffer, (Vars.Coord_Number+2)*(Vars.Neutron_Counter+Vars.Buffer_Block)*sizeof(double));
if (Mon2D_Buffer == NULL)
{ printf("Monitor_nD: %s cannot reallocate Mon2D_Buffer[%li] (%li). Skipping.\n", mccompcurname, i, (Vars.Neutron_Counter+Vars.Buffer_Block)*sizeof(double)); Vars.Flag_List = 1; }
else { Vars.Buffer_Counter = 0; Vars.Buffer_Size = Vars.Neutron_Counter+Vars.Buffer_Block; }
}
while (!While_End)
{ /* we generate Coord[] and Coord_index[] from Buffer (auto limits) or passing neutron */
if (Vars.Flag_Auto_Limits == 2)
{
if (While_Buffer < Vars.Buffer_Block)
{
/* first while loops (While_Buffer) */
/* auto limits case : scan Buffer within limits and store in Mon2D */
for (i = 1; i <= Vars.Coord_Number; i++)
{
/* scanning variables in Buffer */
XY = (Vars.Coord_Max[i]-Vars.Coord_Min[i]);
Coord[i] = Mon2D_Buffer[While_Buffer*(Vars.Coord_Number+2) + (i-1)];
Coord[0] = Mon2D_Buffer[While_Buffer*(Vars.Coord_Number+2) + (Vars.Coord_Number)];
pp = Coord[0];
if (XY > 0) Coord_Index[i] = floor((Mon2D_Buffer[(i-1) + While_Buffer*(Vars.Coord_Number+2)]-Vars.Coord_Min[i])*Vars.Coord_Bin[i]/XY);
else Coord_Index[i] = 0;
if (Vars.Flag_With_Borders)
{
if (Coord_Index[i] < 0) Coord_Index[i] = 0;
if (Coord_Index[i] >= Vars.Coord_Bin[i]) Coord_Index[i] = Vars.Coord_Bin[i] - 1;
}
} /* end for */
While_Buffer++;
} /* end if in Buffer */
else /* (While_Buffer >= Vars.Buffer_Block) && (Vars.Flag_Auto_Limits == 2) */
{
Vars.Flag_Auto_Limits = 0;
if (!Vars.Flag_List) /* free Buffer not needed (no list to output) */
{ /* Dim : (Vars.Coord_Number+2)*Vars.Buffer_Block matrix (for p, dp) */
free(Mon2D_Buffer);
}
}
}
else /* Vars.Flag_Auto_Limits == 0 or 1 */
{
for (i = 0; i <= Vars.Coord_Number; i++)
{ /* handle current neutron : last while */
pp = Vars.cp;
XY = 0;
Set_Vars_Coord_Type = Vars.Coord_Type[i];
if (Set_Vars_Coord_Type == DEFS.COORD_X) XY = Vars.cx;
else
if (Set_Vars_Coord_Type == DEFS.COORD_Y) XY = Vars.cy;
else
if (Set_Vars_Coord_Type == DEFS.COORD_Z) XY = Vars.cz;
else
if (Set_Vars_Coord_Type == DEFS.COORD_VX) XY = Vars.cvx;
else
if (Set_Vars_Coord_Type == DEFS.COORD_VY) XY = Vars.cvy;
else
if (Set_Vars_Coord_Type == DEFS.COORD_VZ) XY = Vars.cvz;
else
if (Set_Vars_Coord_Type == DEFS.COORD_KX) XY = V2K*Vars.cvx;
else
if (Set_Vars_Coord_Type == DEFS.COORD_KY) XY = V2K*Vars.cvy;
else
if (Set_Vars_Coord_Type == DEFS.COORD_KZ) XY = V2K*Vars.cvz;
else
if (Set_Vars_Coord_Type == DEFS.COORD_SX) XY = Vars.csx;
else
if (Set_Vars_Coord_Type == DEFS.COORD_SY) XY = Vars.csy;
else
if (Set_Vars_Coord_Type == DEFS.COORD_SZ) XY = Vars.csz;
else
if (Set_Vars_Coord_Type == DEFS.COORD_T) XY = Vars.ct;
else
if (Set_Vars_Coord_Type == DEFS.COORD_P) XY = pp;
else
if (Set_Vars_Coord_Type == DEFS.COORD_HDIV) XY = RAD2DEG*atan2(Vars.cvx,Vars.cvz);
else
if (Set_Vars_Coord_Type == DEFS.COORD_VDIV) XY = RAD2DEG*atan2(Vars.cvy,Vars.cvz);
else
if (Set_Vars_Coord_Type == DEFS.COORD_V) XY = sqrt(Vars.cvx*Vars.cvx+Vars.cvy*Vars.cvy+Vars.cvz*Vars.cvz);
else
if (Set_Vars_Coord_Type == DEFS.COORD_RADIUS) XY = sqrt(Vars.cx*Vars.cx+Vars.cy*Vars.cy);
else
if (Set_Vars_Coord_Type == DEFS.COORD_K) { XY = sqrt(Vars.cvx*Vars.cvx+Vars.cvy*Vars.cvy+Vars.cvz*Vars.cvz); XY *= V2K; }
else
if (Set_Vars_Coord_Type == DEFS.COORD_ENERGY) { XY = Vars.cvx*Vars.cvx+Vars.cvy*Vars.cvy+Vars.cvz*Vars.cvz; XY *= VS2E; }
else
if (Set_Vars_Coord_Type == DEFS.COORD_LAMBDA) { XY = sqrt(Vars.cvx*Vars.cvx+Vars.cvy*Vars.cvy+Vars.cvz*Vars.cvz); XY *= V2K; if (XY != 0) XY = 2*PI/XY; }
else
if (Set_Vars_Coord_Type == DEFS.COORD_NCOUNT) XY = Coord[i]+1;
else
if (Set_Vars_Coord_Type == DEFS.COORD_ANGLE)
{ XY = sqrt(Vars.cvx*Vars.cvx+Vars.cvy*Vars.cvy+Vars.cvz*Vars.cvz);
if (Vars.cvz != 0)
{
XY= RAD2DEG*atan2(XY,Vars.cvz);
} else XY = 0;
}
else
if (Set_Vars_Coord_Type == DEFS.COORD_THETA) { if (Vars.cz != 0) XY = RAD2DEG*atan2(Vars.cx,Vars.cz); }
else
if (Set_Vars_Coord_Type == DEFS.COORD_PHI) { if (Vars.cz != 0) XY = RAD2DEG*atan2(Vars.cy,Vars.cz); }
else
if (Set_Vars_Coord_Type == DEFS.COORD_USER1) XY = Vars.UserVariable1;
else
if (Set_Vars_Coord_Type == DEFS.COORD_USER2) XY = Vars.UserVariable2;
else
XY = 0;
Coord[i] = XY;
if (!Vars.Flag_Auto_Limits)
{
XY = (Vars.Coord_Max[i]-Vars.Coord_Min[i]);
if (XY > 0) Coord_Index[i] = floor((Coord[i]-Vars.Coord_Min[i])*Vars.Coord_Bin[i]/XY);
else Coord_Index[i] = 0;
if (Vars.Flag_With_Borders)
{
if (Coord_Index[i] < 0) Coord_Index[i] = 0;
if (Coord_Index[i] >= Vars.Coord_Bin[i]) Coord_Index[i] = Vars.Coord_Bin[i] - 1;
}
} /* else Auto_Limits will get Index later from Buffer */
} /* end for i */
While_End = 1;
} /* end else if Vars.Flag_Auto_Limits == 2 */
if (Vars.Flag_Auto_Limits != 2) /* not when reading auto limits Buffer */
{ /* now store Coord into Buffer (no index needed) if necessary */
if ((Vars.Buffer_Counter < Vars.Buffer_Block) && ((Vars.Flag_List) || (Vars.Flag_Auto_Limits == 1)))
{
for (i = 0; i < Vars.Coord_Number; i++)
{
Mon2D_Buffer[i + Vars.Neutron_Counter*(Vars.Coord_Number+2)] = Coord[i+1];
}
Mon2D_Buffer[Vars.Coord_Number + Vars.Neutron_Counter*(Vars.Coord_Number+2)] = pp;
Mon2D_Buffer[(Vars.Coord_Number+1) + Vars.Neutron_Counter*(Vars.Coord_Number+2)] = pp*pp;
Vars.Buffer_Counter++;
if (Vars.Flag_Verbose && (Vars.Buffer_Counter >= Vars.Buffer_Block) && (Vars.Flag_List == 1)) printf("Monitor_nD: %s %li neutrons stored in List.\n", mccompcurname, Vars.Buffer_Counter);
}
Vars.Neutron_Counter++;
} /* end (Vars.Flag_Auto_Limits != 2) */
/* store n1d/2d section for Buffer or current neutron in while */
if (Vars.Flag_Auto_Limits != 1) /* not when storing auto limits Buffer */
{
/* 1D and n1D case : Vars.Flag_Multiple */
if (Vars.Flag_Multiple)
{ /* Dim : Vars.Coord_Number*Vars.Coord_Bin[i] vectors (intensity is not included) */
for (i= 0; i < Vars.Coord_Number; i++)
{
j = Coord_Index[i+1];
if (j >= 0 && j < Vars.Coord_Bin[i+1])
{
Mon2D_N[i][j]++;
Mon2D_p[i][j] += pp;
Mon2D_p2[i][j] += pp*pp;
}
}
}
else /* 2D case : Vars.Coord_Number==2 and !Vars.Flag_Multiple and !Vars.Flag_List */
if ((Vars.Coord_Number == 2) && !Vars.Flag_Multiple)
{ /* Dim : Vars.Coord_Bin[1]*Vars.Coord_Bin[2] matrix */
i = Coord_Index[1];
j = Coord_Index[2];
if (i >= 0 && i < Vars.Coord_Bin[1] && j >= 0 && j < Vars.Coord_Bin[2])
{
Mon2D_N[i][j]++;
Mon2D_p[i][j] += pp;
Mon2D_p2[i][j] += pp*pp;
}
}
} /* end (Vars.Flag_Auto_Limits != 1) */
} /* end while */
} /* end if Vars.Coord_Number */
/* 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(Nsum, psum, p2sum, Mon2D_N, Mon2D_p, Mon2D_p2, Mon2D_Buffer, DEFS, Vars, 0, mccompcurname);
}
} /* end if intersection */
else if (Vars.Flag_Absorb) ABSORB;
%}
FINALLY
%{
/* save results, and free pointers */
Monitor_nD_OutPut(Nsum, psum, p2sum, Mon2D_N, Mon2D_p, Mon2D_p2, Mon2D_Buffer, DEFS, Vars, 1, mccompcurname);
%}
MCDISPLAY
%{
double radius, h;
radius = Vars.Sphere_Radius;
h = Vars.Cylinder_Height;
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);
}
%}
END
More information about the mcstas-users
mailing list