Mcstas example programs
stuart
stuart at studsvik.uu.se
Fri Feb 12 14:07:59 CET 1999
Hi
Thankyou for your reply. The files are attached.
Stuart Rycroft
Kristian Nielsen wrote:
> We have some alpha machines at Risø, so if you could send me the
> prisma2.c file generated by McStas and the prisma2 executable produced
> by the C compiler, I could try to repeat your problems and find the
> cause.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: prisma2.out
Type: application/x-unknown-content-type-out_auto_file
Size: 82124 bytes
Desc: not available
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/19990212/0ae9d47c/attachment.bin>
-------------- next part --------------
/* Automatically generated file. Do not edit. */
#define MC_USE_DEFAULT_MAIN
#define MC_EMBEDDED_RUNTIME
#line 1 "mcstas-r.h"
/*******************************************************************************
* Runtime system for McStas.
*
* Project: Monte Carlo Simulation of Triple Axis Spectrometers
* File name: mcstas-r.h
*
* Author: K.N. Aug 29, 1997
*
* $Id: mcstas-r.h,v 1.20 1998/10/09 07:53:48 kn Exp kn $
*
* $Log: mcstas-r.h,v $
* Revision 1.20 1998/10/09 07:53:48 kn
* Added some unit conversion constants.
*
* Revision 1.19 1998/10/02 08:38:36 kn
* Added DETECTOR_OUT support.
* Fixed header comment.
*
* Revision 1.18 1998/10/01 08:12:42 kn
* Support for embedding the file in the output from McStas.
* Added mcstas_main() function.
* Added support for command line arguments.
*
* Revision 1.17 1998/09/24 13:01:39 kn
* Minor conversion factor additions.
*
* Revision 1.16 1998/09/23 13:52:08 kn
* Added conversion factors.
* McStas now uses its own random() implementation (unless
* USE_SYSTEM_RANDOM is defined).
*
* Revision 1.15 1998/05/19 07:59:45 kn
* Hack to make random number generation work with HP's CC C compiler.
*
* Revision 1.14 1998/04/17 11:50:31 kn
* Added sphere_intersect.
*
* Revision 1.13 1998/04/17 10:53:08 kn
* Added randvec_target_sphere.
*
* Revision 1.12 1998/03/25 07:23:24 kn
* Fixed RAND_MAX #define for HPUX.
*
* Revision 1.11 1998/03/24 13:59:26 lefmann
* Added #define for RAND_MAX, needed on HPUX.
*
* Revision 1.10 1998/03/24 13:24:40 lefmann
* Added HBAR, MNEUTRON.
*
* Revision 1.9 1998/03/24 07:42:35 kn
* Added definition for PI.
*
* Revision 1.8 1998/03/24 07:36:20 kn
* Make ABSORB macro work better in control structures.
* Add test_printf().
* Add rand01(), randpm1(), rand0max(), and randminmax().
* Add PROP_X0, PROP_Y0, PROP_DT, vec_prod(), scalar_prod(), NORM(), and
* rotate().
* Fix typos.
*
* Revision 1.7 1998/03/20 14:20:10 lefmann
* Added a few definitions.
*
* Revision 1.6 1998/03/18 13:21:48 elu_krni
* Added definition of PROP_Z0 macro.
*
* Revision 1.5 1998/03/16 08:04:16 kn
* Added normal distributed random number function randnorm().
*
* Revision 1.4 1997/12/03 13:34:19 kn
* Added definition of ABSORB macro.
*
* Revision 1.3 1997/10/16 14:27:28 kn
* Added debugging output used by the "display" graphical visualization
* tool.
*
* Revision 1.2 1997/09/08 11:31:27 kn
* Added mcsetstate() function.
*
* Revision 1.1 1997/09/08 10:39:44 kn
* Initial revision
*
*
* 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
typedef double MCNUM;
typedef struct {MCNUM x, y, z;} Coords;
typedef MCNUM Rotation[3][3];
struct mcinputtable_struct {
char *name;
MCNUM *par;
};
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);
#define ABSORB do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \
mcnlt,mcnls1,mcnls2, mcnlp); mcDEBUG_ABSORB(); goto mcabsorb;} while(0)
#define MC_GETPAR(comp, par) mcc ## comp ## _ ## par
#define DETECTOR_OUT(p0,p1,p2) printf("Detector: %s_I=%g %s_ERR=%g\n", \
mccompcurname, p1, mccompcurname, mcestimate_error(p0,p1,p2))
#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_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_LEAVE()
#define mcDEBUG_ABSORB()
#endif
#ifdef TEST
#define test_printf printf
#else
#define test_printf while(0) printf
#endif
#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] to sqrt(E)[meV] */
#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);
#ifndef USE_SYSTEM_RANDOM
#ifdef RAND_MAX
# undef RAND_MAX
#endif
#define RAND_MAX 0x7fffffff
#define random mc_random
#define srandom mc_srandom
#endif /* !USE_SYSTEM_RANDOM */
#define rand01() ( ((double)random())/((double)RAND_MAX+1) )
#define randpm1() ( ((double)random()) / (((double)RAND_MAX+1)/2) - 1 )
#define rand0max(max) ( ((double)random()) / (((double)RAND_MAX+1)/(max)) )
#define randminmax(min,max) ( rand0max((max)-(min)) - (min) )
#define PROP_X0 \
{ \
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; \
}
#define PROP_Y0 \
{ \
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; \
}
#define PROP_Z0 \
{ \
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; \
}
#define PROP_DT(dt) \
{ \
mcnlx += mcnlvx*(dt); \
mcnly += mcnlvy*(dt); \
mcnlz += mcnlvz*(dt); \
mcnlt += (dt); \
}
#define vec_prod(x, y, z, x1, y1, z1, x2, y2, z2) \
{ \
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; \
}
#define scalar_prod(x1, y1, z1, x2, y2, z2) \
((x1)*(x2) + (y1)*(y2) + (z1)*(z2))
#define NORM(x,y,z) \
{ \
double mcnm_tmp = sqrt((x)*(x) + (y)*(y) + (z)*(z)); \
if(mcnm_tmp != 0.0) \
{ \
(x) /= mcnm_tmp; \
(y) /= mcnm_tmp; \
(z) /= mcnm_tmp; \
} \
}
#define rotate(x, y, z, vx, vy, vz, phi, ax, ay, az) \
{ \
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; \
}
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);
double mcestimate_error(int N, double p1, double p2);
void mcreadparams(void);
void mcsetstate(double x, double y, double z, double vx, double vy, double vz,
double t, double s1, double s2, double p);
void mcgenstate(void);
double randnorm(void);
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 mcset_ncount(double count);
int mcstas_main(int argc, char *argv[]);
#endif /* MCSTAS_R_H */
/* End of file "mcstas-r.h". */
#line 1 "mcstas-r.c"
/*******************************************************************************
* Runtime system for McStas.
*
* Project: Monte Carlo Simulation of Triple Axis Spectrometers
* File name: mcstas-r.c
*
* Author: K.N. Aug 27, 1997
*
* $Id: mcstas-r.c,v 1.15 1998/10/02 08:38:27 kn Exp kn $
*
* $Log: mcstas-r.c,v $
* Revision 1.15 1998/10/02 08:38:27 kn
* Added DETECTOR_OUT support.
* Fixed header comment.
*
* Revision 1.14 1998/10/01 08:12:26 kn
* Support for embedding the file in the output from McStas.
* Added mcstas_main() function.
* Added support for command line arguments.
*
* Revision 1.13 1998/09/23 13:51:35 kn
* McStas now uses its own random() implementation (unless
* USE_SYSTEM_RANDOM is defined).
*
* Revision 1.12 1998/05/19 07:54:05 kn
* In randvec_target_sphere, accept radius=0 to mean that no focusing is to
* be done (choose direction uniformly in full 4PI solid angle).
*
* Revision 1.11 1998/05/11 12:08:49 kn
* Fix bug in cylinder_intersect that caused an infinite cylinder height in
* some cases.
*
* Revision 1.10 1998/04/17 11:50:08 kn
* Added sphere_intersect.
*
* Revision 1.9 1998/04/17 10:52:27 kn
* Better names in randvec_target_sphere.
*
* Revision 1.8 1998/04/16 14:21:49 kn
* Added randvec_target() function.
*
* Revision 1.7 1998/03/24 07:34:03 kn
* Use rand01() in randnorm(). Fix typos.
*
* Revision 1.6 1998/03/20 14:19:52 lefmann
* Added cylinder_intersect().
*
* Revision 1.5 1998/03/16 08:03:41 kn
* Added normal distributed random number function randnorm().
*
* Revision 1.4 1997/10/16 14:27:05 kn
* Add missing #include. Change in mcreadparams() to fit better with the
* "display" visualization tool.
*
* Revision 1.3 1997/09/08 11:31:22 kn
* Added mcsetstate() function.
*
* Revision 1.2 1997/09/08 11:16:43 kn
* Bug fix in mccoordschange().
*
* Revision 1.1 1997/09/08 10:40:03 kn
* Initial revision
*
*
* Copyright (C) Risoe National Laboratory, 1997-1998, All rights reserved
*******************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#ifndef MCSTAS_R_H
#include "mcstas-r.h"
#endif
#ifdef MC_ANCIENT_COMPATIBILITY
int mctraceenabled = 0;
int mcdefaultmain = 0;
#endif
/* 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? */
}
double mcestimate_error(int N, double p1, double p2)
{
double pmean, n1;
if(N <= 1)
return 0;
pmean = p1 / N;
n1 = N - 1;
return sqrt((N/n1)*(p2 - pmean*pmean));
}
void
mcreadparams(void)
{
extern struct mcinputtable_struct mcinputtable[];
int i;
for(i = 0; mcinputtable[i].name != 0; i++)
{
printf("Set value of instrument parameter %s:\n", mcinputtable[i].name);
fflush(stdout);
scanf("%lf", mcinputtable[i].par);
}
}
void
mcsetstate(double x, double y, double z, double vx, double vy, double vz,
double t, double s1, double s2, double p)
{
extern double mcnx, mcny, mcnz, mcnvx, mcnvy, mcnvz;
extern double mcnt, mcns1, mcns2, mcnp;
mcnx = x;
mcny = y;
mcnz = z;
mcnvx = vx;
mcnvy = vy;
mcnvz = vz;
mcnt = t;
mcns1 = s1;
mcns2 = s2;
mcnp = p;
}
void
mcgenstate(void)
{
mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 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 ();
}
}
/* 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;
}
/* 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);
}
/* Number of neutron histories to simulate. */
static double mcncount = 1e6;
void
mcset_ncount(double count)
{
mcncount = count;
}
mcstatic int mcdotrace = 0;
static void
mcsetn_arg(char *arg)
{
mcset_ncount(strtod(arg, NULL));
}
static void
mcsetseed(char *arg)
{
srandom(atol(arg));
}
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\n"
" -n count --ncount=count Set number of neutrons to simulate.\n"
" -t --trace Enable trace of neutron through instrument\n"
" -h --help Show help message\n"
" -i --info Detailed instrument information\n"
"Instrument parameters are:\n");
for(i = 0; i < mcnumipar; i++)
fprintf(stderr, " %s\n", mcinputtable[i].name);
}
static void
mcshowhelp(char *pgmname)
{
mchelp(pgmname);
exit(0);
}
static void
mcusage(char *pgmname)
{
fprintf(stderr, "Error: incorrect commad line arguments\n");
mchelp(pgmname);
exit(1);
}
static void
mcinfo(void)
{
int i;
printf("Name: %s\n", mcinstrument_name);
printf("Parameters:");
for(i = 0; i < mcnumipar; i++)
printf(" %s", mcinputtable[i].name);
printf("\n");
printf("Instrument-source: %s\n", mcinstrument_source);
printf("Trace-enabled: %s\n", mctraceenabled ? "yes" : "no");
printf("Default-main: %s\n", mcdefaultmain ? "yes" : "no");
printf("Embedded-runtime: %s\n",
#ifdef MC_EMBEDDED_RUNTIME
"yes"
#else
"no"
#endif
);
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\n");
exit(1);
}
}
void
mcparseoptions(int argc, char *argv[])
{
int i, j, pos;
char *p;
int paramset = 0, *paramsetarray;
paramsetarray = malloc(mcnumipar*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("-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(argv[i][0] != '-' && (p = strchr(argv[i], '=')) != NULL)
{
*p++ = '\0';
for(j = 0; j < mcnumipar; j++)
if(!strcmp(mcinputtable[j].name, argv[i]))
{
*(mcinputtable[j].par) = strtod(p, NULL);
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);
}
}
}
/* McStas main() function. */
int
mcstas_main(int argc, char *argv[])
{
int run_num = 0;
srandom(time(NULL));
mcparseoptions(argc, argv);
mcinit();
while(run_num < mcncount)
{
mcsetstate(0, 0, 0, 0, 0, 1, 0, 0, 0, 1);
mcraytrace();
run_num++;
}
mcfinally();
return 0;
}
/* End of file "mcstas-r.c". */
int mctraceenabled = 0;
int mcdefaultmain = 1;
char mcinstrument_name[] = "prisma2";
char mcinstrument_source[] = "examples/prisma2.instr";
int main(int argc, char *argv[]){return mcstas_main(argc, argv);}
void mcinit(void);
void mcraytrace(void);
void mcfinally(void);
/* Instrument parameters. */
MCNUM mcipTT;
MCNUM mcipPHA;
MCNUM mcipPHA1;
MCNUM mcipPHA2;
MCNUM mcipPHA3;
MCNUM mcipPHA4;
MCNUM mcipPHA5;
MCNUM mcipPHA6;
MCNUM mcipPHA7;
MCNUM mcipTTA;
#define mcNUMIPAR 10
int mcnumipar = 10;
struct mcinputtable_struct mcinputtable[mcNUMIPAR+1] = {
"TT", &mcipTT,
"PHA", &mcipPHA,
"PHA1", &mcipPHA1,
"PHA2", &mcipPHA2,
"PHA3", &mcipPHA3,
"PHA4", &mcipPHA4,
"PHA5", &mcipPHA5,
"PHA6", &mcipPHA6,
"PHA7", &mcipPHA7,
"TTA", &mcipTTA,
NULL, NULL
};
/* User declarations from instrument definition. */
#define TT mcipTT
#define PHA mcipPHA
#define PHA1 mcipPHA1
#define PHA2 mcipPHA2
#define PHA3 mcipPHA3
#define PHA4 mcipPHA4
#define PHA5 mcipPHA5
#define PHA6 mcipPHA6
#define PHA7 mcipPHA7
#define TTA mcipTTA
#line 171 "examples/prisma2.instr"
int neu_color; /* "Color" of current neutron */
/* 30' mosaicity used on analysator */
double prisma_ana_mosaic = 30;
/* Q vector for bragg scattering with monochromator and analysator */
double prisma_ana_q = 1.87325;
double prisma_ana_r0 = 0.6;
double focus_x,focus_z;
double apos1, apos2, apos3, apos4, apos5, apos6, apos7;
#undef TTA
#undef PHA7
#undef PHA6
#undef PHA5
#undef PHA4
#undef PHA3
#undef PHA2
#undef PHA1
#undef PHA
#undef TT
/* Declarations of component definition and setting parameters. */
/* Definition parameters for component 'mod'. */
#define mccmod_radius 0.0707
#define mccmod_E0 10
#define mccmod_E1 15
#define mccmod_dist 9.035
#define mccmod_xw 0.021
#define mccmod_yh 0.021
#define mccmod_t0 37.15
#define mccmod_Ec 9
#define mccmod_gam 39.1
/* Definition parameters for component 'modslit'. */
#define mccmodslit_xmin -0.05
#define mccmodslit_xmax 0.05
#define mccmodslit_ymin -0.05
#define mccmodslit_ymax 0.05
/* Definition parameters for component 'tof_test'. */
#define mcctof_test_xmin -0.05
#define mcctof_test_xmax 0.05
#define mcctof_test_ymin -0.05
#define mcctof_test_ymax 0.05
#define mcctof_test_nchan 10000
#define mcctof_test_dt 10
#define mcctof_test_filename "sim/prisma2.mon"
/* Definition parameters for component 'mon1'. */
#define mccmon1_xmin -0.1
#define mccmon1_xmax 0.1
#define mccmon1_ymin -0.1
#define mccmon1_ymax 0.1
/* Definition parameters for component 'slit1'. */
#define mccslit1_xmin -0.05
#define mccslit1_xmax 0.05
#define mccslit1_ymin -0.05
#define mccslit1_ymax 0.05
/* Definition parameters for component 'slit2'. */
#define mccslit2_xmin -0.02
#define mccslit2_xmax 0.02
#define mccslit2_ymin -0.03
#define mccslit2_ymax 0.03
/* Definition parameters for component 'mon2'. */
#define mccmon2_xmin -0.1
#define mccmon2_xmax 0.1
#define mccmon2_ymin -0.1
#define mccmon2_ymax 0.1
/* Definition parameters for component 'sample'. */
#define mccsample_radius_i 1e-05
#define mccsample_radius_o 0.01
#define mccsample_h 0.02
#define mccsample_pack 1
#define mccsample_focus_r 0.03
/* Setting parameters for component 'sample'. */
MCNUM mccsample_target_x;
MCNUM mccsample_target_y;
MCNUM mccsample_target_z;
/* Definition parameters for component 'mon3'. */
#define mccmon3_xmin -0.1
#define mccmon3_xmax 0.1
#define mccmon3_ymin -0.1
#define mccmon3_ymax 0.1
/* Definition parameters for component 'coll2'. */
#define mcccoll2_xmin -0.015
#define mcccoll2_xmax 0.015
#define mcccoll2_ymin -0.025
#define mcccoll2_ymax 0.025
#define mcccoll2_len 0.12
#define mcccoll2_divergence 120
/* Definition parameters for component 'mon4'. */
#define mccmon4_xmin -0.1
#define mccmon4_xmax 0.1
#define mccmon4_ymin -0.1
#define mccmon4_ymax 0.1
/* Definition parameters for component 'ana1'. */
#define mccana1_zmin -0.006
#define mccana1_zmax 0.006
#define mccana1_ymin -0.0375
#define mccana1_ymax 0.0375
#define mccana1_mosaich prisma_ana_mosaic
#define mccana1_mosaicv prisma_ana_mosaic
#define mccana1_r0 prisma_ana_r0
#define mccana1_Q prisma_ana_q
#define mccana1_color 0
/* Definition parameters for component 'ana2'. */
#define mccana2_zmin -0.006
#define mccana2_zmax 0.006
#define mccana2_ymin -0.0375
#define mccana2_ymax 0.0375
#define mccana2_mosaich prisma_ana_mosaic
#define mccana2_mosaicv prisma_ana_mosaic
#define mccana2_r0 prisma_ana_r0
#define mccana2_Q prisma_ana_q
#define mccana2_color 1
/* Definition parameters for component 'ana3'. */
#define mccana3_zmin -0.006
#define mccana3_zmax 0.006
#define mccana3_ymin -0.0375
#define mccana3_ymax 0.0375
#define mccana3_mosaich prisma_ana_mosaic
#define mccana3_mosaicv prisma_ana_mosaic
#define mccana3_r0 prisma_ana_r0
#define mccana3_Q prisma_ana_q
#define mccana3_color 2
/* Definition parameters for component 'ana4'. */
#define mccana4_zmin -0.006
#define mccana4_zmax 0.006
#define mccana4_ymin -0.0375
#define mccana4_ymax 0.0375
#define mccana4_mosaich prisma_ana_mosaic
#define mccana4_mosaicv prisma_ana_mosaic
#define mccana4_r0 prisma_ana_r0
#define mccana4_Q prisma_ana_q
#define mccana4_color 3
/* Definition parameters for component 'ana5'. */
#define mccana5_zmin -0.006
#define mccana5_zmax 0.006
#define mccana5_ymin -0.0375
#define mccana5_ymax 0.0375
#define mccana5_mosaich prisma_ana_mosaic
#define mccana5_mosaicv prisma_ana_mosaic
#define mccana5_r0 prisma_ana_r0
#define mccana5_Q prisma_ana_q
#define mccana5_color 4
/* Definition parameters for component 'ana6'. */
#define mccana6_zmin -0.006
#define mccana6_zmax 0.006
#define mccana6_ymin -0.0375
#define mccana6_ymax 0.0375
#define mccana6_mosaich prisma_ana_mosaic
#define mccana6_mosaicv prisma_ana_mosaic
#define mccana6_r0 prisma_ana_r0
#define mccana6_Q prisma_ana_q
#define mccana6_color 5
/* Definition parameters for component 'ana7'. */
#define mccana7_zmin -0.006
#define mccana7_zmax 0.006
#define mccana7_ymin -0.0375
#define mccana7_ymax 0.0375
#define mccana7_mosaich prisma_ana_mosaic
#define mccana7_mosaicv prisma_ana_mosaic
#define mccana7_r0 prisma_ana_r0
#define mccana7_Q prisma_ana_q
#define mccana7_color 6
/* Definition parameters for component 'mon5'. */
#define mccmon5_xmin -0.05
#define mccmon5_xmax 0.05
#define mccmon5_ymin -0.05
#define mccmon5_ymax 0.05
/* Definition parameters for component 'mon6'. */
#define mccmon6_xmin -0.1
#define mccmon6_xmax 0.1
#define mccmon6_ymin -0.1
#define mccmon6_ymax 0.1
/* Definition parameters for component 'psd'. */
#define mccpsd_xmin -0.05
#define mccpsd_xmax 0.05
#define mccpsd_ymin -0.05
#define mccpsd_ymax 0.05
#define mccpsd_nx 100
#define mccpsd_ny 100
#define mccpsd_filename "sim/prisma2.psd"
/* Definition parameters for component 'detector'. */
#define mccdetector_xmin -0.05
#define mccdetector_xmax 0.05
#define mccdetector_ymin -0.05
#define mccdetector_ymax 0.05
#define mccdetector_nchan 10000
#define mccdetector_dt 10
#define mccdetector_filename "sim/prisma2.tof"
#define mccdetector_maxcolor 6
/* Definition parameters for component 'mon9'. */
#define mccmon9_xmin -0.1
#define mccmon9_xmax 0.1
#define mccmon9_ymin -0.1
#define mccmon9_ymax 0.1
/* User component declarations. */
/* User declarations for component 'mod'. */
#define mccompcurname "mod"
#define radius mccmod_radius
#define E0 mccmod_E0
#define E1 mccmod_E1
#define dist mccmod_dist
#define xw mccmod_xw
#define yh mccmod_yh
#define t0 mccmod_t0
#define Ec mccmod_Ec
#define gam mccmod_gam
#line 33 "/usr/users/batman/mc01/lib/mcstas/Moderator.comp"
double hdiv,vdiv; /* ToDo: Should be component local */
double moder_v0, moder_dv;
double p_in;
#undef gam
#undef Ec
#undef t0
#undef yh
#undef xw
#undef dist
#undef E1
#undef E0
#undef radius
#undef mccompcurname
/* User declarations for component 'tof_test'. */
#define mccompcurname "tof_test"
#define xmin mcctof_test_xmin
#define xmax mcctof_test_xmax
#define ymin mcctof_test_ymin
#define ymax mcctof_test_ymax
#define nchan mcctof_test_nchan
#define dt mcctof_test_dt
#define filename mcctof_test_filename
#define TOF_N mcctof_test_TOF_N
#define TOF_p mcctof_test_TOF_p
#define TOF_p2 mcctof_test_TOF_p2
#line 39 "/usr/users/batman/mc01/lib/mcstas/TOF_monitor.comp"
int TOF_N[nchan];
double TOF_p[nchan];
double TOF_p2[nchan];
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'mon1'. */
#define mccompcurname "mon1"
#define xmin mccmon1_xmin
#define xmax mccmon1_xmax
#define ymin mccmon1_ymin
#define ymax mccmon1_ymax
#define Nsum mccmon1_Nsum
#define psum mccmon1_psum
#define p2sum mccmon1_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'mon2'. */
#define mccompcurname "mon2"
#define xmin mccmon2_xmin
#define xmax mccmon2_xmax
#define ymin mccmon2_ymin
#define ymax mccmon2_ymax
#define Nsum mccmon2_Nsum
#define psum mccmon2_psum
#define p2sum mccmon2_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'sample'. */
#define mccompcurname "sample"
#define radius_i mccsample_radius_i
#define radius_o mccsample_radius_o
#define h mccsample_h
#define pack mccsample_pack
#define focus_r mccsample_focus_r
#define target_x mccsample_target_x
#define target_y mccsample_target_y
#define target_z mccsample_target_z
#line 38 "/usr/users/batman/mc01/lib/mcstas/V_sample.comp"
/* ToDo: Should be component local names. */
#define V_sigma_a 5.08 /* Absorption cross section per atom (barns) */
#define V_sigma_i 4.935 /* Incoherent scattering cross section per atom (barns) */
#define V_rho (2*pack/(3.024*3.024*3.024)) /* Density of atoms (AA-3) */
#define V_my_s (V_rho * 100 * V_sigma_i)
#define V_my_a_v (V_rho * 100 * V_sigma_a * 2200)
#undef target_z
#undef target_y
#undef target_x
#undef focus_r
#undef pack
#undef h
#undef radius_o
#undef radius_i
#undef mccompcurname
/* User declarations for component 'mon3'. */
#define mccompcurname "mon3"
#define xmin mccmon3_xmin
#define xmax mccmon3_xmax
#define ymin mccmon3_ymin
#define ymax mccmon3_ymax
#define Nsum mccmon3_Nsum
#define psum mccmon3_psum
#define p2sum mccmon3_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'coll2'. */
#define mccompcurname "coll2"
#define xmin mcccoll2_xmin
#define xmax mcccoll2_xmax
#define ymin mcccoll2_ymin
#define ymax mcccoll2_ymax
#define len mcccoll2_len
#define divergence mcccoll2_divergence
#define slope mcccoll2_slope
#line 36 "/usr/users/batman/mc01/lib/mcstas/Soller.comp"
double slope;
#undef slope
#undef divergence
#undef len
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'mon4'. */
#define mccompcurname "mon4"
#define xmin mccmon4_xmin
#define xmax mccmon4_xmax
#define ymin mccmon4_ymin
#define ymax mccmon4_ymax
#define Nsum mccmon4_Nsum
#define psum mccmon4_psum
#define p2sum mccmon4_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'ana1'. */
#define mccompcurname "ana1"
#define zmin mccana1_zmin
#define zmax mccana1_zmax
#define ymin mccana1_ymin
#define ymax mccana1_ymax
#define mosaich mccana1_mosaich
#define mosaicv mccana1_mosaicv
#define r0 mccana1_r0
#define Q mccana1_Q
#define color mccana1_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'ana2'. */
#define mccompcurname "ana2"
#define zmin mccana2_zmin
#define zmax mccana2_zmax
#define ymin mccana2_ymin
#define ymax mccana2_ymax
#define mosaich mccana2_mosaich
#define mosaicv mccana2_mosaicv
#define r0 mccana2_r0
#define Q mccana2_Q
#define color mccana2_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'ana3'. */
#define mccompcurname "ana3"
#define zmin mccana3_zmin
#define zmax mccana3_zmax
#define ymin mccana3_ymin
#define ymax mccana3_ymax
#define mosaich mccana3_mosaich
#define mosaicv mccana3_mosaicv
#define r0 mccana3_r0
#define Q mccana3_Q
#define color mccana3_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'ana4'. */
#define mccompcurname "ana4"
#define zmin mccana4_zmin
#define zmax mccana4_zmax
#define ymin mccana4_ymin
#define ymax mccana4_ymax
#define mosaich mccana4_mosaich
#define mosaicv mccana4_mosaicv
#define r0 mccana4_r0
#define Q mccana4_Q
#define color mccana4_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'ana5'. */
#define mccompcurname "ana5"
#define zmin mccana5_zmin
#define zmax mccana5_zmax
#define ymin mccana5_ymin
#define ymax mccana5_ymax
#define mosaich mccana5_mosaich
#define mosaicv mccana5_mosaicv
#define r0 mccana5_r0
#define Q mccana5_Q
#define color mccana5_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'ana6'. */
#define mccompcurname "ana6"
#define zmin mccana6_zmin
#define zmax mccana6_zmax
#define ymin mccana6_ymin
#define ymax mccana6_ymax
#define mosaich mccana6_mosaich
#define mosaicv mccana6_mosaicv
#define r0 mccana6_r0
#define Q mccana6_Q
#define color mccana6_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'ana7'. */
#define mccompcurname "ana7"
#define zmin mccana7_zmin
#define zmax mccana7_zmax
#define ymin mccana7_ymin
#define ymax mccana7_ymax
#define mosaich mccana7_mosaich
#define mosaicv mccana7_mosaicv
#define r0 mccana7_r0
#define Q mccana7_Q
#define color mccana7_color
#line 29 "examples/prisma2.instr"
#define DIV_CUTOFF 2 /* ~ 10^-5 cutoff. */
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
/* User declarations for component 'mon5'. */
#define mccompcurname "mon5"
#define xmin mccmon5_xmin
#define xmax mccmon5_xmax
#define ymin mccmon5_ymin
#define ymax mccmon5_ymax
#define Nsum mccmon5_Nsum
#define psum mccmon5_psum
#define p2sum mccmon5_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'mon6'. */
#define mccompcurname "mon6"
#define xmin mccmon6_xmin
#define xmax mccmon6_xmax
#define ymin mccmon6_ymin
#define ymax mccmon6_ymax
#define Nsum mccmon6_Nsum
#define psum mccmon6_psum
#define p2sum mccmon6_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'psd'. */
#define mccompcurname "psd"
#define xmin mccpsd_xmin
#define xmax mccpsd_xmax
#define ymin mccpsd_ymin
#define ymax mccpsd_ymax
#define nx mccpsd_nx
#define ny mccpsd_ny
#define filename mccpsd_filename
#define PSD_N mccpsd_PSD_N
#define PSD_p mccpsd_PSD_p
#define PSD_p2 mccpsd_PSD_p2
#line 40 "/usr/users/batman/mc01/lib/mcstas/PSD_monitor.comp"
int PSD_N[nx][ny];
double PSD_p[nx][ny];
double PSD_p2[nx][ny];
#undef PSD_p2
#undef PSD_p
#undef PSD_N
#undef filename
#undef ny
#undef nx
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'detector'. */
#define mccompcurname "detector"
#define xmin mccdetector_xmin
#define xmax mccdetector_xmax
#define ymin mccdetector_ymin
#define ymax mccdetector_ymax
#define nchan mccdetector_nchan
#define dt mccdetector_dt
#define filename mccdetector_filename
#define maxcolor mccdetector_maxcolor
#define TOF_N mccdetector_TOF_N
#define TOF_p mccdetector_TOF_p
#define TOF_p2 mccdetector_TOF_p2
#line 99 "examples/prisma2.instr"
int TOF_N[maxcolor+1][nchan];
double TOF_p[maxcolor+1][nchan];
double TOF_p2[maxcolor+1][nchan];
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef maxcolor
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User declarations for component 'mon9'. */
#define mccompcurname "mon9"
#define xmin mccmon9_xmin
#define xmax mccmon9_xmax
#define ymin mccmon9_ymin
#define ymax mccmon9_ymax
#define Nsum mccmon9_Nsum
#define psum mccmon9_psum
#define p2sum mccmon9_p2sum
#line 36 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
int Nsum;
double psum, p2sum;
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
Coords mcposamod, mcposrmod;
Rotation mcrotamod, mcrotrmod;
Coords mcposamodslit, mcposrmodslit;
Rotation mcrotamodslit, mcrotrmodslit;
Coords mcposatof_test, mcposrtof_test;
Rotation mcrotatof_test, mcrotrtof_test;
Coords mcposamon1, mcposrmon1;
Rotation mcrotamon1, mcrotrmon1;
Coords mcposaslit1, mcposrslit1;
Rotation mcrotaslit1, mcrotrslit1;
Coords mcposaslit2, mcposrslit2;
Rotation mcrotaslit2, mcrotrslit2;
Coords mcposamon2, mcposrmon2;
Rotation mcrotamon2, mcrotrmon2;
Coords mcposasample, mcposrsample;
Rotation mcrotasample, mcrotrsample;
Coords mcposaa2, mcposra2;
Rotation mcrotaa2, mcrotra2;
Coords mcposamon3, mcposrmon3;
Rotation mcrotamon3, mcrotrmon3;
Coords mcposacoll2, mcposrcoll2;
Rotation mcrotacoll2, mcrotrcoll2;
Coords mcposamon4, mcposrmon4;
Rotation mcrotamon4, mcrotrmon4;
Coords mcposarita_ana, mcposrrita_ana;
Rotation mcrotarita_ana, mcrotrrita_ana;
Coords mcposaana1, mcposrana1;
Rotation mcrotaana1, mcrotrana1;
Coords mcposaana2, mcposrana2;
Rotation mcrotaana2, mcrotrana2;
Coords mcposaana3, mcposrana3;
Rotation mcrotaana3, mcrotrana3;
Coords mcposaana4, mcposrana4;
Rotation mcrotaana4, mcrotrana4;
Coords mcposaana5, mcposrana5;
Rotation mcrotaana5, mcrotrana5;
Coords mcposaana6, mcposrana6;
Rotation mcrotaana6, mcrotrana6;
Coords mcposaana7, mcposrana7;
Rotation mcrotaana7, mcrotrana7;
Coords mcposaa3, mcposra3;
Rotation mcrotaa3, mcrotra3;
Coords mcposamon5, mcposrmon5;
Rotation mcrotamon5, mcrotrmon5;
Coords mcposamon6, mcposrmon6;
Rotation mcrotamon6, mcrotrmon6;
Coords mcposapsd, mcposrpsd;
Rotation mcrotapsd, mcrotrpsd;
Coords mcposadetector, mcposrdetector;
Rotation mcrotadetector, mcrotrdetector;
Coords mcposamon9, mcposrmon9;
Rotation mcrotamon9, mcrotrmon9;
MCNUM mcnx, mcny, mcnz, mcnvx, mcnvy, mcnvz, mcnt, mcns1, mcns2, mcnp;
void mcinit(void) {
#define TT mcipTT
#define PHA mcipPHA
#define PHA1 mcipPHA1
#define PHA2 mcipPHA2
#define PHA3 mcipPHA3
#define PHA4 mcipPHA4
#define PHA5 mcipPHA5
#define PHA6 mcipPHA6
#define PHA7 mcipPHA7
#define TTA mcipTTA
#line 184 "examples/prisma2.instr"
{
focus_x = 0.52 * sin(TT*DEG2RAD);
focus_z = 0.52 * cos(TT*DEG2RAD);
/* Rita-style analyser. */
{
double l = 0.0125;
apos1 = -3*l;
apos2 = -2*l;
apos3 = -1*l;
apos4 = 0*l;
apos5 = 1*l;
apos6 = 2*l;
apos7 = 3*l;
}
}
#undef TTA
#undef PHA7
#undef PHA6
#undef PHA5
#undef PHA4
#undef PHA3
#undef PHA2
#undef PHA1
#undef PHA
#undef TT
/* Component initializations. */
/* Initializations for component mod. */
#define mccompcurname "mod"
#define radius mccmod_radius
#define E0 mccmod_E0
#define E1 mccmod_E1
#define dist mccmod_dist
#define xw mccmod_xw
#define yh mccmod_yh
#define t0 mccmod_t0
#define Ec mccmod_Ec
#define gam mccmod_gam
#line 38 "/usr/users/batman/mc01/lib/mcstas/Moderator.comp"
{
hdiv = atan(xw/(2.0*dist));
vdiv = atan(yh/(2.0*dist));
moder_v0 = SE2V*sqrt(E0);
moder_dv = SE2V*sqrt(E1) - moder_v0;
p_in = (4*hdiv*vdiv)/(4*PI);
}
#undef gam
#undef Ec
#undef t0
#undef yh
#undef xw
#undef dist
#undef E1
#undef E0
#undef radius
#undef mccompcurname
/* Initializations for component modslit. */
/* Initializations for component tof_test. */
#define mccompcurname "tof_test"
#define xmin mcctof_test_xmin
#define xmax mcctof_test_xmax
#define ymin mcctof_test_ymin
#define ymax mcctof_test_ymax
#define nchan mcctof_test_nchan
#define dt mcctof_test_dt
#define filename mcctof_test_filename
#define TOF_N mcctof_test_TOF_N
#define TOF_p mcctof_test_TOF_p
#define TOF_p2 mcctof_test_TOF_p2
#line 44 "/usr/users/batman/mc01/lib/mcstas/TOF_monitor.comp"
{
int i;
for (i=0; i<nchan; i++)
{
TOF_N[i] = 0;
TOF_p[i] = 0;
TOF_p2[i] = 0;
}
}
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component mon1. */
#define mccompcurname "mon1"
#define xmin mccmon1_xmin
#define xmax mccmon1_xmax
#define ymin mccmon1_ymin
#define ymax mccmon1_ymax
#define Nsum mccmon1_Nsum
#define psum mccmon1_psum
#define p2sum mccmon1_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component slit1. */
/* Initializations for component slit2. */
/* Initializations for component mon2. */
#define mccompcurname "mon2"
#define xmin mccmon2_xmin
#define xmax mccmon2_xmax
#define ymin mccmon2_ymin
#define ymax mccmon2_ymax
#define Nsum mccmon2_Nsum
#define psum mccmon2_psum
#define p2sum mccmon2_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component sample. */
mccsample_target_x = focus_x;
mccsample_target_y = 0;
mccsample_target_z = focus_z;
/* Initializations for component a2. */
/* Initializations for component mon3. */
#define mccompcurname "mon3"
#define xmin mccmon3_xmin
#define xmax mccmon3_xmax
#define ymin mccmon3_ymin
#define ymax mccmon3_ymax
#define Nsum mccmon3_Nsum
#define psum mccmon3_psum
#define p2sum mccmon3_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component coll2. */
#define mccompcurname "coll2"
#define xmin mcccoll2_xmin
#define xmax mcccoll2_xmax
#define ymin mcccoll2_ymin
#define ymax mcccoll2_ymax
#define len mcccoll2_len
#define divergence mcccoll2_divergence
#define slope mcccoll2_slope
#line 39 "/usr/users/batman/mc01/lib/mcstas/Soller.comp"
{
slope = tan(MIN2RAD*divergence);
}
#undef slope
#undef divergence
#undef len
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component mon4. */
#define mccompcurname "mon4"
#define xmin mccmon4_xmin
#define xmax mccmon4_xmax
#define ymin mccmon4_ymin
#define ymax mccmon4_ymax
#define Nsum mccmon4_Nsum
#define psum mccmon4_psum
#define p2sum mccmon4_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component rita_ana. */
/* Initializations for component ana1. */
/* Initializations for component ana2. */
/* Initializations for component ana3. */
/* Initializations for component ana4. */
/* Initializations for component ana5. */
/* Initializations for component ana6. */
/* Initializations for component ana7. */
/* Initializations for component a3. */
/* Initializations for component mon5. */
#define mccompcurname "mon5"
#define xmin mccmon5_xmin
#define xmax mccmon5_xmax
#define ymin mccmon5_ymin
#define ymax mccmon5_ymax
#define Nsum mccmon5_Nsum
#define psum mccmon5_psum
#define p2sum mccmon5_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component mon6. */
#define mccompcurname "mon6"
#define xmin mccmon6_xmin
#define xmax mccmon6_xmax
#define ymin mccmon6_ymin
#define ymax mccmon6_ymax
#define Nsum mccmon6_Nsum
#define psum mccmon6_psum
#define p2sum mccmon6_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component psd. */
#define mccompcurname "psd"
#define xmin mccpsd_xmin
#define xmax mccpsd_xmax
#define ymin mccpsd_ymin
#define ymax mccpsd_ymax
#define nx mccpsd_nx
#define ny mccpsd_ny
#define filename mccpsd_filename
#define PSD_N mccpsd_PSD_N
#define PSD_p mccpsd_PSD_p
#define PSD_p2 mccpsd_PSD_p2
#line 45 "/usr/users/batman/mc01/lib/mcstas/PSD_monitor.comp"
{
int i,j;
for (i=0; i<nx; i++)
for (j=0; j<ny; j++)
{
PSD_N[i][j] = 0;
PSD_p[i][j] = 0;
PSD_p2[i][j] = 0;
}
}
#undef PSD_p2
#undef PSD_p
#undef PSD_N
#undef filename
#undef ny
#undef nx
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component detector. */
#define mccompcurname "detector"
#define xmin mccdetector_xmin
#define xmax mccdetector_xmax
#define ymin mccdetector_ymin
#define ymax mccdetector_ymax
#define nchan mccdetector_nchan
#define dt mccdetector_dt
#define filename mccdetector_filename
#define maxcolor mccdetector_maxcolor
#define TOF_N mccdetector_TOF_N
#define TOF_p mccdetector_TOF_p
#define TOF_p2 mccdetector_TOF_p2
#line 104 "examples/prisma2.instr"
{
int i,c;
for (i=0; i<nchan; i++)
for (c=0; c<=maxcolor; c++)
{
TOF_N[c][i] = 0;
TOF_p[c][i] = 0;
TOF_p2[c][i] = 0;
}
}
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef maxcolor
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Initializations for component mon9. */
#define mccompcurname "mon9"
#define xmin mccmon9_xmin
#define xmax mccmon9_xmax
#define ymin mccmon9_ymin
#define ymax mccmon9_ymax
#define Nsum mccmon9_Nsum
#define psum mccmon9_psum
#define p2sum mccmon9_p2sum
#line 40 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
psum = 0;
p2sum = 0;
Nsum = 0;
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* Computation of coordinate transformations. */
{
Coords mctc1, mctc2;
Rotation mctr1, mctr2;
mcDEBUG_INSTR()
/* Component mod. */
rot_set_rotation(mcrotamod, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_copy(mcrotrmod, mcrotamod);
mcposamod = coords_set(0, 0, 0);
mctc1 = coords_neg(mcposamod);
mcposrmod = rot_apply(mcrotamod, mctc1);
mcDEBUG_COMPONENT("mod", mcposamod, mcrotamod)
/* Component modslit. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotamodslit);
rot_transpose(mcrotamod, mctr1);
rot_mul(mcrotamodslit, mctr1, mcrotrmodslit);
mctc1 = coords_set(0, 0, 0);
rot_transpose(mcrotamod, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamodslit = coords_add(mcposamod, mctc2);
mctc1 = coords_sub(mcposamod, mcposamodslit);
mcposrmodslit = rot_apply(mcrotamodslit, mctc1);
mcDEBUG_COMPONENT("modslit", mcposamodslit, mcrotamodslit)
/* Component tof_test. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotatof_test);
rot_transpose(mcrotamodslit, mctr1);
rot_mul(mcrotatof_test, mctr1, mcrotrtof_test);
mctc1 = coords_set(0, 0, 0.005);
rot_transpose(mcrotamod, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposatof_test = coords_add(mcposamod, mctc2);
mctc1 = coords_sub(mcposamodslit, mcposatof_test);
mcposrtof_test = rot_apply(mcrotatof_test, mctc1);
mcDEBUG_COMPONENT("tof_test", mcposatof_test, mcrotatof_test)
/* Component mon1. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotamon1);
rot_transpose(mcrotatof_test, mctr1);
rot_mul(mcrotamon1, mctr1, mcrotrmon1);
mctc1 = coords_set(0, 0, 0.01);
rot_transpose(mcrotamod, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon1 = coords_add(mcposamod, mctc2);
mctc1 = coords_sub(mcposatof_test, mcposamon1);
mcposrmon1 = rot_apply(mcrotamon1, mctc1);
mcDEBUG_COMPONENT("mon1", mcposamon1, mcrotamon1)
/* Component slit1. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotaslit1);
rot_transpose(mcrotamon1, mctr1);
rot_mul(mcrotaslit1, mctr1, mcrotrslit1);
mctc1 = coords_set(0, 0, 1.7);
rot_transpose(mcrotamod, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaslit1 = coords_add(mcposamod, mctc2);
mctc1 = coords_sub(mcposamon1, mcposaslit1);
mcposrslit1 = rot_apply(mcrotaslit1, mctc1);
mcDEBUG_COMPONENT("slit1", mcposaslit1, mcrotaslit1)
/* Component slit2. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotaslit2);
rot_transpose(mcrotaslit1, mctr1);
rot_mul(mcrotaslit2, mctr1, mcrotrslit2);
mctc1 = coords_set(0, 0, 7);
rot_transpose(mcrotaslit1, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaslit2 = coords_add(mcposaslit1, mctc2);
mctc1 = coords_sub(mcposaslit1, mcposaslit2);
mcposrslit2 = rot_apply(mcrotaslit2, mctc1);
mcDEBUG_COMPONENT("slit2", mcposaslit2, mcrotaslit2)
/* Component mon2. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotamon2);
rot_transpose(mcrotaslit2, mctr1);
rot_mul(mcrotamon2, mctr1, mcrotrmon2);
mctc1 = coords_set(0, 0, 9);
rot_transpose(mcrotamod, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon2 = coords_add(mcposamod, mctc2);
mctc1 = coords_sub(mcposaslit2, mcposamon2);
mcposrmon2 = rot_apply(mcrotamon2, mctc1);
mcDEBUG_COMPONENT("mon2", mcposamon2, mcrotamon2)
/* Component sample. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotamod, mcrotasample);
rot_transpose(mcrotamon2, mctr1);
rot_mul(mcrotasample, mctr1, mcrotrsample);
mctc1 = coords_set(0, 0, 9.035);
rot_transpose(mcrotamod, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposasample = coords_add(mcposamod, mctc2);
mctc1 = coords_sub(mcposamon2, mcposasample);
mcposrsample = rot_apply(mcrotasample, mctc1);
mcDEBUG_COMPONENT("sample", mcposasample, mcrotasample)
/* Component a2. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipTT)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotasample, mcrotaa2);
rot_transpose(mcrotasample, mctr1);
rot_mul(mcrotaa2, mctr1, mcrotra2);
mctc1 = coords_set(0, 0, 0);
rot_transpose(mcrotasample, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaa2 = coords_add(mcposasample, mctc2);
mctc1 = coords_sub(mcposasample, mcposaa2);
mcposra2 = rot_apply(mcrotaa2, mctc1);
mcDEBUG_COMPONENT("a2", mcposaa2, mcrotaa2)
/* Component mon3. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa2, mcrotamon3);
rot_transpose(mcrotaa2, mctr1);
rot_mul(mcrotamon3, mctr1, mcrotrmon3);
mctc1 = coords_set(0, 0, 0.39);
rot_transpose(mcrotaa2, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon3 = coords_add(mcposaa2, mctc2);
mctc1 = coords_sub(mcposaa2, mcposamon3);
mcposrmon3 = rot_apply(mcrotamon3, mctc1);
mcDEBUG_COMPONENT("mon3", mcposamon3, mcrotamon3)
/* Component coll2. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa2, mcrotacoll2);
rot_transpose(mcrotamon3, mctr1);
rot_mul(mcrotacoll2, mctr1, mcrotrcoll2);
mctc1 = coords_set(0, 0, 0.4);
rot_transpose(mcrotaa2, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposacoll2 = coords_add(mcposaa2, mctc2);
mctc1 = coords_sub(mcposamon3, mcposacoll2);
mcposrcoll2 = rot_apply(mcrotacoll2, mctc1);
mcDEBUG_COMPONENT("coll2", mcposacoll2, mcrotacoll2)
/* Component mon4. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa2, mcrotamon4);
rot_transpose(mcrotacoll2, mctr1);
rot_mul(mcrotamon4, mctr1, mcrotrmon4);
mctc1 = coords_set(0, 0, 0.521);
rot_transpose(mcrotaa2, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon4 = coords_add(mcposaa2, mctc2);
mctc1 = coords_sub(mcposacoll2, mcposamon4);
mcposrmon4 = rot_apply(mcrotamon4, mctc1);
mcDEBUG_COMPONENT("mon4", mcposamon4, mcrotamon4)
/* Component rita_ana. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa2, mcrotarita_ana);
rot_transpose(mcrotamon4, mctr1);
rot_mul(mcrotarita_ana, mctr1, mcrotrrita_ana);
mctc1 = coords_set(0, 0, 0.58);
rot_transpose(mcrotaa2, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposarita_ana = coords_add(mcposaa2, mctc2);
mctc1 = coords_sub(mcposamon4, mcposarita_ana);
mcposrrita_ana = rot_apply(mcrotarita_ana, mctc1);
mcDEBUG_COMPONENT("rita_ana", mcposarita_ana, mcrotarita_ana)
/* Component ana1. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA1)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana1);
rot_transpose(mcrotarita_ana, mctr1);
rot_mul(mcrotaana1, mctr1, mcrotrana1);
mctc1 = coords_set(0, 0, apos1);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana1 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposarita_ana, mcposaana1);
mcposrana1 = rot_apply(mcrotaana1, mctc1);
mcDEBUG_COMPONENT("ana1", mcposaana1, mcrotaana1)
/* Component ana2. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA2)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana2);
rot_transpose(mcrotaana1, mctr1);
rot_mul(mcrotaana2, mctr1, mcrotrana2);
mctc1 = coords_set(0, 0, apos2);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana2 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana1, mcposaana2);
mcposrana2 = rot_apply(mcrotaana2, mctc1);
mcDEBUG_COMPONENT("ana2", mcposaana2, mcrotaana2)
/* Component ana3. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA3)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana3);
rot_transpose(mcrotaana2, mctr1);
rot_mul(mcrotaana3, mctr1, mcrotrana3);
mctc1 = coords_set(0, 0, apos3);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana3 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana2, mcposaana3);
mcposrana3 = rot_apply(mcrotaana3, mctc1);
mcDEBUG_COMPONENT("ana3", mcposaana3, mcrotaana3)
/* Component ana4. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA4)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana4);
rot_transpose(mcrotaana3, mctr1);
rot_mul(mcrotaana4, mctr1, mcrotrana4);
mctc1 = coords_set(0, 0, apos4);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana4 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana3, mcposaana4);
mcposrana4 = rot_apply(mcrotaana4, mctc1);
mcDEBUG_COMPONENT("ana4", mcposaana4, mcrotaana4)
/* Component ana5. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA5)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana5);
rot_transpose(mcrotaana4, mctr1);
rot_mul(mcrotaana5, mctr1, mcrotrana5);
mctc1 = coords_set(0, 0, apos5);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana5 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana4, mcposaana5);
mcposrana5 = rot_apply(mcrotaana5, mctc1);
mcDEBUG_COMPONENT("ana5", mcposaana5, mcrotaana5)
/* Component ana6. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA6)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana6);
rot_transpose(mcrotaana5, mctr1);
rot_mul(mcrotaana6, mctr1, mcrotrana6);
mctc1 = coords_set(0, 0, apos6);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana6 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana5, mcposaana6);
mcposrana6 = rot_apply(mcrotaana6, mctc1);
mcDEBUG_COMPONENT("ana6", mcposaana6, mcrotaana6)
/* Component ana7. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHA7)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotarita_ana, mcrotaana7);
rot_transpose(mcrotaana6, mctr1);
rot_mul(mcrotaana7, mctr1, mcrotrana7);
mctc1 = coords_set(0, 0, apos7);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaana7 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana6, mcposaana7);
mcposrana7 = rot_apply(mcrotaana7, mctc1);
mcDEBUG_COMPONENT("ana7", mcposaana7, mcrotaana7)
/* Component a3. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipTTA)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa2, mcrotaa3);
rot_transpose(mcrotaana7, mctr1);
rot_mul(mcrotaa3, mctr1, mcrotra3);
mctc1 = coords_set(0, 0, 0);
rot_transpose(mcrotarita_ana, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposaa3 = coords_add(mcposarita_ana, mctc2);
mctc1 = coords_sub(mcposaana7, mcposaa3);
mcposra3 = rot_apply(mcrotaa3, mctc1);
mcDEBUG_COMPONENT("a3", mcposaa3, mcrotaa3)
/* Component mon5. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa3, mcrotamon5);
rot_transpose(mcrotaa3, mctr1);
rot_mul(mcrotamon5, mctr1, mcrotrmon5);
mctc1 = coords_set(0, 0, 0.06);
rot_transpose(mcrotaa3, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon5 = coords_add(mcposaa3, mctc2);
mctc1 = coords_sub(mcposaa3, mcposamon5);
mcposrmon5 = rot_apply(mcrotamon5, mctc1);
mcDEBUG_COMPONENT("mon5", mcposamon5, mcrotamon5)
/* Component mon6. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa3, mcrotamon6);
rot_transpose(mcrotamon5, mctr1);
rot_mul(mcrotamon6, mctr1, mcrotrmon6);
mctc1 = coords_set(0, 0, 0.161);
rot_transpose(mcrotaa3, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon6 = coords_add(mcposaa3, mctc2);
mctc1 = coords_sub(mcposamon5, mcposamon6);
mcposrmon6 = rot_apply(mcrotamon6, mctc1);
mcDEBUG_COMPONENT("mon6", mcposamon6, mcrotamon6)
/* Component psd. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa3, mcrotapsd);
rot_transpose(mcrotamon6, mctr1);
rot_mul(mcrotapsd, mctr1, mcrotrpsd);
mctc1 = coords_set(0, 0, 0.2);
rot_transpose(mcrotaa3, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposapsd = coords_add(mcposaa3, mctc2);
mctc1 = coords_sub(mcposamon6, mcposapsd);
mcposrpsd = rot_apply(mcrotapsd, mctc1);
mcDEBUG_COMPONENT("psd", mcposapsd, mcrotapsd)
/* Component detector. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa3, mcrotadetector);
rot_transpose(mcrotapsd, mctr1);
rot_mul(mcrotadetector, mctr1, mcrotrdetector);
mctc1 = coords_set(0, 0, 0.2);
rot_transpose(mcrotaa3, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposadetector = coords_add(mcposaa3, mctc2);
mctc1 = coords_sub(mcposapsd, mcposadetector);
mcposrdetector = rot_apply(mcrotadetector, mctc1);
mcDEBUG_COMPONENT("detector", mcposadetector, mcrotadetector)
/* Component mon9. */
rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
rot_mul(mctr1, mcrotaa3, mcrotamon9);
rot_transpose(mcrotadetector, mctr1);
rot_mul(mcrotamon9, mctr1, mcrotrmon9);
mctc1 = coords_set(0, 0, 0.01);
rot_transpose(mcrotadetector, mctr1);
mctc2 = rot_apply(mctr1, mctc1);
mcposamon9 = coords_add(mcposadetector, mctc2);
mctc1 = coords_sub(mcposadetector, mcposamon9);
mcposrmon9 = rot_apply(mcrotamon9, mctc1);
mcDEBUG_COMPONENT("mon9", mcposamon9, mcrotamon9)
mcDEBUG_INSTR_END()
}
}
void mcraytrace(void) {
/* Copy neutron state to local variables. */
MCNUM mcnlx = mcnx;
MCNUM mcnly = mcny;
MCNUM mcnlz = mcnz;
MCNUM mcnlvx = mcnvx;
MCNUM mcnlvy = mcnvy;
MCNUM mcnlvz = mcnvz;
MCNUM mcnlt = mcnt;
MCNUM mcnls1 = mcns1;
MCNUM mcnls2 = mcns2;
MCNUM mcnlp = mcnp;
mcDEBUG_ENTER()
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mod. */
mcDEBUG_COMP("mod")
mccoordschange(mcposrmod, mcrotrmod,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mod"
#define radius mccmod_radius
#define E0 mccmod_E0
#define E1 mccmod_E1
#define dist mccmod_dist
#define xw mccmod_xw
#define yh mccmod_yh
#define t0 mccmod_t0
#define Ec mccmod_Ec
#define gam mccmod_gam
#line 46 "/usr/users/batman/mc01/lib/mcstas/Moderator.comp"
{
double theta0,phi0,chi,theta,phi,v,r,tauv,E;
p=p_in;
z=0;
chi = 2*PI*rand01(); /* Choose point on source */
r = sqrt(rand01())*radius; /* with uniform distribution. */
x = r*cos(chi);
y = r*sin(chi);
theta0 = -atan(x/dist); /* Angles to aim at target centre */
phi0 = -atan(y/dist);
theta = theta0 + hdiv*randpm1(); /* Small angle approx. */
phi = phi0 + vdiv*randpm1();
v = moder_v0 + moder_dv*rand01(); /* Assume linear distribution */
vz = v*cos(phi)*cos(theta); /* Small angle approx. */
vy = v*sin(phi);
vx = v*cos(phi)*sin(theta);
E = VS2E*v*v;
if(E < Ec)
{
tauv = t0;
}
else
{
double tmp = ((E - Ec) / gam);
tauv = t0 / (1 + (tmp*tmp));
}
t = -tauv*log(rand01())*1E-6;
}
#undef gam
#undef Ec
#undef t0
#undef yh
#undef xw
#undef dist
#undef E1
#undef E0
#undef radius
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component modslit. */
mcDEBUG_COMP("modslit")
mccoordschange(mcposrmodslit, mcrotrmodslit,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "modslit"
#define xmin mccmodslit_xmin
#define xmax mccmodslit_xmax
#define ymin mccmodslit_ymin
#define ymax mccmodslit_ymax
#line 28 "/usr/users/batman/mc01/lib/mcstas/Slit.comp"
{
PROP_Z0;
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB;
}
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component tof_test. */
mcDEBUG_COMP("tof_test")
mccoordschange(mcposrtof_test, mcrotrtof_test,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "tof_test"
#define xmin mcctof_test_xmin
#define xmax mcctof_test_xmax
#define ymin mcctof_test_ymin
#define ymax mcctof_test_ymax
#define nchan mcctof_test_nchan
#define dt mcctof_test_dt
#define filename mcctof_test_filename
#define TOF_N mcctof_test_TOF_N
#define TOF_p mcctof_test_TOF_p
#define TOF_p2 mcctof_test_TOF_p2
#line 55 "/usr/users/batman/mc01/lib/mcstas/TOF_monitor.comp"
{
int i;
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
i = floor(1E6*t/dt); /* Bin number */
if(i >= nchan) i = nchan;
if(i < 0)
{
printf("FATAL ERROR: negative time-of-flight.\n");
exit(1);
}
TOF_N[i]++;
TOF_p[i] += p;
TOF_p2[i] += p*p;
}
}
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon1. */
mcDEBUG_COMP("mon1")
mccoordschange(mcposrmon1, mcrotrmon1,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon1"
#define xmin mccmon1_xmin
#define xmax mccmon1_xmax
#define ymin mccmon1_ymin
#define ymax mccmon1_ymax
#define Nsum mccmon1_Nsum
#define psum mccmon1_psum
#define p2sum mccmon1_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component slit1. */
mcDEBUG_COMP("slit1")
mccoordschange(mcposrslit1, mcrotrslit1,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "slit1"
#define xmin mccslit1_xmin
#define xmax mccslit1_xmax
#define ymin mccslit1_ymin
#define ymax mccslit1_ymax
#line 28 "/usr/users/batman/mc01/lib/mcstas/Slit.comp"
{
PROP_Z0;
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB;
}
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component slit2. */
mcDEBUG_COMP("slit2")
mccoordschange(mcposrslit2, mcrotrslit2,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "slit2"
#define xmin mccslit2_xmin
#define xmax mccslit2_xmax
#define ymin mccslit2_ymin
#define ymax mccslit2_ymax
#line 28 "/usr/users/batman/mc01/lib/mcstas/Slit.comp"
{
PROP_Z0;
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB;
}
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon2. */
mcDEBUG_COMP("mon2")
mccoordschange(mcposrmon2, mcrotrmon2,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon2"
#define xmin mccmon2_xmin
#define xmax mccmon2_xmax
#define ymin mccmon2_ymin
#define ymax mccmon2_ymax
#define Nsum mccmon2_Nsum
#define psum mccmon2_psum
#define p2sum mccmon2_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component sample. */
mcDEBUG_COMP("sample")
mccoordschange(mcposrsample, mcrotrsample,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "sample"
#define radius_i mccsample_radius_i
#define radius_o mccsample_radius_o
#define h mccsample_h
#define pack mccsample_pack
#define focus_r mccsample_focus_r
#define target_x mccsample_target_x
#define target_y mccsample_target_y
#define target_z mccsample_target_z
#line 49 "/usr/users/batman/mc01/lib/mcstas/V_sample.comp"
{
double t0, t3; /* Entry/exit time for outer cylinder */
double t1, t2; /* Entry/exit time for inner cylinder */
double v; /* Neutron velocity */
double dt0, dt1, dt2, dt; /* Flight times through sample */
double l_full; /* Flight path length for non-scattered neutron */
double l_i, l_o; /* Flight path lenght in/out for scattered neutron */
double my_a; /* Velocity-dependent attenuation factor */
double solid_angle; /* Solid angle of target as seen from scattering point */
double aim_x, aim_y, aim_z; /* Position of target relative to scattering point */
if(cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius_o, h))
{
if(t0 < 0)
ABSORB;
/* Neutron enters at t=t0. */
if(!cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius_i, h))
t1 = t2 = t3;
dt0 = t1-t0; /* Time in sample, ingoing */
dt1 = t2-t1; /* Time in hole */
dt2 = t3-t2; /* Time in sample, outgoing */
v = sqrt(vx*vx + vy*vy + vz*vz);
l_full = v * (dt0 + dt2); /* Length of full path through sample */
dt = rand01()*(dt0+dt2); /* Time of scattering (relative to t0) */
l_i = v*dt; /* Penetration in sample */
if (dt > dt0)
dt += dt1;
PROP_DT(dt+t0); /* Point of scattering */
aim_x = target_x-x; /* Vector pointing at target (anal./det.) */
aim_y = target_y-y;
aim_z = target_z-z;
randvec_target_sphere(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r);
NORM(vx, vy, vz);
vx *= v;
vy *= v;
vz *= v;
if(!cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius_o, h))
{
/* ??? did not hit cylinder */
printf("FATAL ERROR: Did not hit cylinder from inside.\n");
exit(1);
}
dt = t3;
if(cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius_i, h) &&
t2 > 0)
dt -= (t2-t1); /* Subtract hollow part */
l_o = v*dt;
my_a = V_my_a_v/v;
p *= l_full*V_my_s*exp(-(my_a+V_my_s)*(l_i+l_o));
p /= 4*PI/solid_angle;
}
else
ABSORB;
}
#undef target_z
#undef target_y
#undef target_x
#undef focus_r
#undef pack
#undef h
#undef radius_o
#undef radius_i
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component a2. */
mcDEBUG_COMP("a2")
mccoordschange(mcposra2, mcrotra2,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define mccompcurname "a2"
#undef mccompcurname
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon3. */
mcDEBUG_COMP("mon3")
mccoordschange(mcposrmon3, mcrotrmon3,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon3"
#define xmin mccmon3_xmin
#define xmax mccmon3_xmax
#define ymin mccmon3_ymin
#define ymax mccmon3_ymax
#define Nsum mccmon3_Nsum
#define psum mccmon3_psum
#define p2sum mccmon3_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component coll2. */
mcDEBUG_COMP("coll2")
mccoordschange(mcposrcoll2, mcrotrcoll2,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "coll2"
#define xmin mcccoll2_xmin
#define xmax mcccoll2_xmax
#define ymin mcccoll2_ymin
#define ymax mcccoll2_ymax
#define len mcccoll2_len
#define divergence mcccoll2_divergence
#define slope mcccoll2_slope
#line 43 "/usr/users/batman/mc01/lib/mcstas/Soller.comp"
{
double phi, dt;
PROP_Z0;
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB;
dt = len/vz;
PROP_DT(dt);
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB;
if(slope > 0.0)
{
phi = fabs(vx/vz);
if (phi > slope)
ABSORB;
else
p *= 1.0 - phi/slope;
}
}
#undef slope
#undef divergence
#undef len
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon4. */
mcDEBUG_COMP("mon4")
mccoordschange(mcposrmon4, mcrotrmon4,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon4"
#define xmin mccmon4_xmin
#define xmax mccmon4_xmax
#define ymin mccmon4_ymin
#define ymax mccmon4_ymax
#define Nsum mccmon4_Nsum
#define psum mccmon4_psum
#define p2sum mccmon4_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component rita_ana. */
mcDEBUG_COMP("rita_ana")
mccoordschange(mcposrrita_ana, mcrotrrita_ana,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define mccompcurname "rita_ana"
#undef mccompcurname
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana1. */
mcDEBUG_COMP("ana1")
mccoordschange(mcposrana1, mcrotrana1,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana1"
#define zmin mccana1_zmin
#define zmax mccana1_zmax
#define ymin mccana1_ymin
#define ymax mccana1_ymax
#define mosaich mccana1_mosaich
#define mosaicv mccana1_mosaicv
#define r0 mccana1_r0
#define Q mccana1_Q
#define color mccana1_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana2. */
mcDEBUG_COMP("ana2")
mccoordschange(mcposrana2, mcrotrana2,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana2"
#define zmin mccana2_zmin
#define zmax mccana2_zmax
#define ymin mccana2_ymin
#define ymax mccana2_ymax
#define mosaich mccana2_mosaich
#define mosaicv mccana2_mosaicv
#define r0 mccana2_r0
#define Q mccana2_Q
#define color mccana2_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana3. */
mcDEBUG_COMP("ana3")
mccoordschange(mcposrana3, mcrotrana3,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana3"
#define zmin mccana3_zmin
#define zmax mccana3_zmax
#define ymin mccana3_ymin
#define ymax mccana3_ymax
#define mosaich mccana3_mosaich
#define mosaicv mccana3_mosaicv
#define r0 mccana3_r0
#define Q mccana3_Q
#define color mccana3_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana4. */
mcDEBUG_COMP("ana4")
mccoordschange(mcposrana4, mcrotrana4,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana4"
#define zmin mccana4_zmin
#define zmax mccana4_zmax
#define ymin mccana4_ymin
#define ymax mccana4_ymax
#define mosaich mccana4_mosaich
#define mosaicv mccana4_mosaicv
#define r0 mccana4_r0
#define Q mccana4_Q
#define color mccana4_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana5. */
mcDEBUG_COMP("ana5")
mccoordschange(mcposrana5, mcrotrana5,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana5"
#define zmin mccana5_zmin
#define zmax mccana5_zmax
#define ymin mccana5_ymin
#define ymax mccana5_ymax
#define mosaich mccana5_mosaich
#define mosaicv mccana5_mosaicv
#define r0 mccana5_r0
#define Q mccana5_Q
#define color mccana5_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana6. */
mcDEBUG_COMP("ana6")
mccoordschange(mcposrana6, mcrotrana6,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana6"
#define zmin mccana6_zmin
#define zmax mccana6_zmax
#define ymin mccana6_ymin
#define ymax mccana6_ymax
#define mosaich mccana6_mosaich
#define mosaicv mccana6_mosaicv
#define r0 mccana6_r0
#define Q mccana6_Q
#define color mccana6_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component ana7. */
mcDEBUG_COMP("ana7")
mccoordschange(mcposrana7, mcrotrana7,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "ana7"
#define zmin mccana7_zmin
#define zmax mccana7_zmax
#define ymin mccana7_ymin
#define ymax mccana7_ymax
#define mosaich mccana7_mosaich
#define mosaicv mccana7_mosaicv
#define r0 mccana7_r0
#define Q mccana7_Q
#define color mccana7_color
#line 33 "examples/prisma2.instr"
{
double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
double old_x = x, old_y = y, old_z = z, old_t = t;
double dt;
if(vx != 0.0 && (dt = -x/vx) >= 0.0)
{
y += vy*dt; z += vz*dt; t += dt; x = 0.0;
if (z>zmin && z<zmax && y>ymin && y<ymax)
{
/* First: scattering in plane */
theta0 = atan2(vx,vz); /* neutron angle to slab */
v = sqrt(vx*vx+vy*vy+vz*vz);
theta = asin(Q2V*Q/(2.0*v)); /* Bragg's law */
if(theta0 < 0)
theta = -theta;
tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
if(tmp3 > DIV_CUTOFF)
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
else
{
p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
tmp1 = 2*theta;
cs = cos(tmp1);
sn = sin(tmp1);
tmp2 = cs*vx - sn*vz;
vy = vy;
vz = cs*vz + sn*vx;
vx = tmp2;
/* Second: scatering out of plane.
Approximation is that Debye-Scherrer cone is a plane */
phi = atan2(vy,vz); /* out-of plane angle */
dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm(); /* MC choice: */
/* Vertical angle of the crystallite */
vy = vz*tan(phi+2*dphi*sin(theta));
vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
vz = vz*vratio;
vy = vy*vratio; /* Renormalize v */
vx = vx*vratio;
neu_color = color;
}
}
else
{
x = old_x; y = old_y; z = old_z; t = old_t;
}
}
}
#undef color
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component a3. */
mcDEBUG_COMP("a3")
mccoordschange(mcposra3, mcrotra3,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define mccompcurname "a3"
#undef mccompcurname
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon5. */
mcDEBUG_COMP("mon5")
mccoordschange(mcposrmon5, mcrotrmon5,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon5"
#define xmin mccmon5_xmin
#define xmax mccmon5_xmax
#define ymin mccmon5_ymin
#define ymax mccmon5_ymax
#define Nsum mccmon5_Nsum
#define psum mccmon5_psum
#define p2sum mccmon5_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon6. */
mcDEBUG_COMP("mon6")
mccoordschange(mcposrmon6, mcrotrmon6,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon6"
#define xmin mccmon6_xmin
#define xmax mccmon6_xmax
#define ymin mccmon6_ymin
#define ymax mccmon6_ymax
#define Nsum mccmon6_Nsum
#define psum mccmon6_psum
#define p2sum mccmon6_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component psd. */
mcDEBUG_COMP("psd")
mccoordschange(mcposrpsd, mcrotrpsd,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "psd"
#define xmin mccpsd_xmin
#define xmax mccpsd_xmax
#define ymin mccpsd_ymin
#define ymax mccpsd_ymax
#define nx mccpsd_nx
#define ny mccpsd_ny
#define filename mccpsd_filename
#define PSD_N mccpsd_PSD_N
#define PSD_p mccpsd_PSD_p
#define PSD_p2 mccpsd_PSD_p2
#line 57 "/usr/users/batman/mc01/lib/mcstas/PSD_monitor.comp"
{
int i,j;
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
i = floor((x - xmin)*nx/(xmax - xmin));
j = floor((y - ymin)*ny/(ymax - ymin));
PSD_N[i][j]++;
PSD_p[i][j] += p;
PSD_p2[i][j] += p*p;
}
}
#undef PSD_p2
#undef PSD_p
#undef PSD_N
#undef filename
#undef ny
#undef nx
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component detector. */
mcDEBUG_COMP("detector")
mccoordschange(mcposrdetector, mcrotrdetector,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "detector"
#define xmin mccdetector_xmin
#define xmax mccdetector_xmax
#define ymin mccdetector_ymin
#define ymax mccdetector_ymax
#define nchan mccdetector_nchan
#define dt mccdetector_dt
#define filename mccdetector_filename
#define maxcolor mccdetector_maxcolor
#define TOF_N mccdetector_TOF_N
#define TOF_p mccdetector_TOF_p
#define TOF_p2 mccdetector_TOF_p2
#line 116 "examples/prisma2.instr"
{
int i;
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
i = floor(1E6*t/dt); /* Bin number */
if(i >= nchan) i = nchan;
if(i < 0)
{
printf("FATAL ERROR: negative time-of-flight.\n");
exit(1);
}
if(neu_color < 0 || neu_color > maxcolor)
{
printf("FATAL ERROR: wrong color neutron.\n");
exit(1);
}
TOF_N[neu_color][i]++;
TOF_p[neu_color][i] += p;
TOF_p2[neu_color][i] += p*p;
}
}
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef maxcolor
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Component mon9. */
mcDEBUG_COMP("mon9")
mccoordschange(mcposrmon9, mcrotrmon9,
&mcnlx, &mcnly, &mcnlz,
&mcnlvx, &mcnlvy, &mcnlvz,
&mcnlt, &mcnls1, &mcnls2);
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnls1
#define s2 mcnls2
#define p mcnlp
#define mccompcurname "mon9"
#define xmin mccmon9_xmin
#define xmax mccmon9_xmax
#define ymin mccmon9_ymin
#define ymax mccmon9_ymax
#define Nsum mccmon9_Nsum
#define psum mccmon9_psum
#define p2sum mccmon9_p2sum
#line 46 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
PROP_Z0;
if (x>xmin && x<xmax && y>ymin && y<ymax)
{
Nsum++;
psum += p;
p2sum += p*p;
}
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
#undef p
#undef s2
#undef s1
#undef t
#undef vz
#undef vy
#undef vx
#undef z
#undef y
#undef x
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
mcabsorb:
mcDEBUG_LEAVE()
mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnls1,mcnls2, mcnlp)
/* Copy neutron state to global variables. */
mcnx = mcnlx;
mcny = mcnly;
mcnz = mcnlz;
mcnvx = mcnlvx;
mcnvy = mcnlvy;
mcnvz = mcnlvz;
mcnt = mcnlt;
mcns1 = mcnls1;
mcns2 = mcnls2;
mcnp = mcnlp;
}
void mcfinally(void) {
/* User component FINALLY code. */
/* User FINALLY code for component 'tof_test'. */
#define mccompcurname "tof_test"
#define xmin mcctof_test_xmin
#define xmax mcctof_test_xmax
#define ymin mcctof_test_ymin
#define ymax mcctof_test_ymax
#define nchan mcctof_test_nchan
#define dt mcctof_test_dt
#define filename mcctof_test_filename
#define TOF_N mcctof_test_TOF_N
#define TOF_p mcctof_test_TOF_p
#define TOF_p2 mcctof_test_TOF_p2
#line 74 "/usr/users/batman/mc01/lib/mcstas/TOF_monitor.comp"
{
int i, Nsum;
double psum, p2sum;
FILE *outfile;
Nsum = psum = p2sum = 0;
outfile=fopen(filename,"w");
if(!outfile)
{
fprintf(stderr,
"FATAL ERROR: could not open output file '%s'\n", filename);
exit(1);
}
for (i=0; i<nchan; i++)
{
Nsum += TOF_N[i];
psum += TOF_p[i];
p2sum += TOF_p2[i];
if(TOF_p[i] != 0.0 || TOF_N[i] != 0)
fprintf(outfile,"%g %d %g %g\n",
(double)i*dt, TOF_N[i], TOF_p[i], TOF_p2[i]);
}
fclose(outfile);
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon1'. */
#define mccompcurname "mon1"
#define xmin mccmon1_xmin
#define xmax mccmon1_xmax
#define ymin mccmon1_ymin
#define ymax mccmon1_ymax
#define Nsum mccmon1_Nsum
#define psum mccmon1_psum
#define p2sum mccmon1_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon2'. */
#define mccompcurname "mon2"
#define xmin mccmon2_xmin
#define xmax mccmon2_xmax
#define ymin mccmon2_ymin
#define ymax mccmon2_ymax
#define Nsum mccmon2_Nsum
#define psum mccmon2_psum
#define p2sum mccmon2_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon3'. */
#define mccompcurname "mon3"
#define xmin mccmon3_xmin
#define xmax mccmon3_xmax
#define ymin mccmon3_ymin
#define ymax mccmon3_ymax
#define Nsum mccmon3_Nsum
#define psum mccmon3_psum
#define p2sum mccmon3_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon4'. */
#define mccompcurname "mon4"
#define xmin mccmon4_xmin
#define xmax mccmon4_xmax
#define ymin mccmon4_ymin
#define ymax mccmon4_ymax
#define Nsum mccmon4_Nsum
#define psum mccmon4_psum
#define p2sum mccmon4_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon5'. */
#define mccompcurname "mon5"
#define xmin mccmon5_xmin
#define xmax mccmon5_xmax
#define ymin mccmon5_ymin
#define ymax mccmon5_ymax
#define Nsum mccmon5_Nsum
#define psum mccmon5_psum
#define p2sum mccmon5_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon6'. */
#define mccompcurname "mon6"
#define xmin mccmon6_xmin
#define xmax mccmon6_xmax
#define ymin mccmon6_ymin
#define ymax mccmon6_ymax
#define Nsum mccmon6_Nsum
#define psum mccmon6_psum
#define p2sum mccmon6_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'psd'. */
#define mccompcurname "psd"
#define xmin mccpsd_xmin
#define xmax mccpsd_xmax
#define ymin mccpsd_ymin
#define ymax mccpsd_ymax
#define nx mccpsd_nx
#define ny mccpsd_ny
#define filename mccpsd_filename
#define PSD_N mccpsd_PSD_N
#define PSD_p mccpsd_PSD_p
#define PSD_p2 mccpsd_PSD_p2
#line 71 "/usr/users/batman/mc01/lib/mcstas/PSD_monitor.comp"
{
int i,j;
FILE *outfile;
int Nsum;
double Psum, P2sum;
Nsum = Psum = P2sum = 0;
outfile=fopen(filename,"w");
for (j=0; j<ny; j++)
{
for (i=0; i<nx; i++)
{
fprintf(outfile,"%g ", PSD_p[i][j]);
Nsum += PSD_N[i][j];
Psum += PSD_p[i][j];
P2sum += PSD_p2[i][j];
}
fprintf(outfile,"\n");
}
DETECTOR_OUT(Nsum, Psum, P2sum);
fclose(outfile);
}
#undef PSD_p2
#undef PSD_p
#undef PSD_N
#undef filename
#undef ny
#undef nx
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'detector'. */
#define mccompcurname "detector"
#define xmin mccdetector_xmin
#define xmax mccdetector_xmax
#define ymin mccdetector_ymin
#define ymax mccdetector_ymax
#define nchan mccdetector_nchan
#define dt mccdetector_dt
#define filename mccdetector_filename
#define maxcolor mccdetector_maxcolor
#define TOF_N mccdetector_TOF_N
#define TOF_p mccdetector_TOF_p
#define TOF_p2 mccdetector_TOF_p2
#line 140 "examples/prisma2.instr"
{
int i,c;
FILE *outfile;
outfile=fopen(filename,"w");
for (i=0; i<nchan; i++)
{
int empty = 1;
for(c=0; c<=maxcolor; c++)
if(TOF_p[c][i] != 0.0)
{
empty = 0;
break;
}
if(empty)
continue;
fprintf(outfile,"%g", (double)i*dt);
for(c=0; c<=maxcolor; c++)
fprintf(outfile," %g", TOF_p[c][i]);
fprintf(outfile,"\n");
}
fclose(outfile);
}
#undef TOF_p2
#undef TOF_p
#undef TOF_N
#undef maxcolor
#undef filename
#undef dt
#undef nchan
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
/* User FINALLY code for component 'mon9'. */
#define mccompcurname "mon9"
#define xmin mccmon9_xmin
#define xmax mccmon9_xmax
#define ymin mccmon9_ymin
#define ymax mccmon9_ymax
#define Nsum mccmon9_Nsum
#define psum mccmon9_psum
#define p2sum mccmon9_p2sum
#line 56 "/usr/users/batman/mc01/lib/mcstas/Monitor.comp"
{
DETECTOR_OUT(Nsum, psum, p2sum);
}
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname
}
More information about the mcstas-users
mailing list