Mcstas 1.1

Bertrand.Roessli at psi.ch Bertrand.Roessli at psi.ch
Thu Jun 10 13:53:56 CEST 1999


Dear Kristian,

please find below the linup-1.c file which
has got undefined variables.

regards,

B. Roessli

/* 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 335 "linup-1.c"

#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". */

#line 1084 "linup-1.c"
#ifdef MC_TRACE_ENABLED
int mctraceenabled = 1;
#else
int mctraceenabled = 0;
#endif
int mcdefaultmain = 1;
char mcinstrument_name[] = "TAS1";
char mcinstrument_source[] = "linup-1.instr";
int main(int argc, char *argv[]){return mcstas_main(argc, argv);}
void mcinit(void);
void mcraytrace(void);
void mcfinally(void);
void mcdisplay(void);

/* Instrument parameters. */
MCNUM mcipPHM;
MCNUM mcipTTM;
MCNUM mcipC1;

#define mcNUMIPAR 3
int mcnumipar = 3;
struct mcinputtable_struct mcinputtable[mcNUMIPAR+1] = {
  "PHM", &mcipPHM,
  "TTM", &mcipTTM,
  "C1", &mcipC1,
  NULL, NULL
};

/* User declarations from instrument definition. */
#define PHM mcipPHM
#define TTM mcipTTM
#define C1 mcipC1
#line 5 "linup-1.instr"
/* Mosaicity used on monochromator and analysator */
double tas1_mono_mosaic = 45;	/* Measurements indicate its really 45' */
/* Q vector for bragg scattering with monochromator and analysator */
double tas1_mono_q = 3.354;	/* Fake 2nd order scattering for 20meV */
double tas1_mono_r0 = 0.6;

double mpos0, mpos1, mpos2, mpos3, mpos4, mpos5, mpos6, mpos7;
double mrot0, mrot1, mrot2, mrot3, mrot4, mrot5, mrot6, mrot7;
#line 1126 "linup-1.c"
#undef C1
#undef TTM
#undef PHM

/* Declarations of component definition and setting parameters. */

/* Definition parameters for component 'source'. */
#define mccsource_radius 0.06
#define mccsource_dist 3.288
#define mccsource_xw 0.042
#define mccsource_yh 0.082
#define mccsource_E0 20
#define mccsource_dE 0.82

/* Definition parameters for component 'slit1'. */
#define mccslit1_xmin -0.02
#define mccslit1_xmax 0.065
#define mccslit1_ymin -0.075
#define mccslit1_ymax 0.075

/* Definition parameters for component 'slit2'. */
#define mccslit2_xmin -0.02
#define mccslit2_xmax 0.02
#define mccslit2_ymin -0.04
#define mccslit2_ymax 0.04

/* Definition parameters for component 'slit3'. */
#define mccslit3_xmin -0.021
#define mccslit3_xmax 0.021
#define mccslit3_ymin -0.041
#define mccslit3_ymax 0.041

/* Definition parameters for component 'm0'. */
#define mccm0_zmin -0.0375
#define mccm0_zmax 0.0375
#define mccm0_ymin -0.006
#define mccm0_ymax 0.006
#define mccm0_mosaich tas1_mono_mosaic
#define mccm0_mosaicv tas1_mono_mosaic
#define mccm0_r0 tas1_mono_r0
#define mccm0_Q tas1_mono_q

/* Definition parameters for component 'm1'. */
#define mccm1_zmin -0.0375
#define mccm1_zmax 0.0375
#define mccm1_ymin -0.006
#define mccm1_ymax 0.006
#define mccm1_mosaich tas1_mono_mosaic
#define mccm1_mosaicv tas1_mono_mosaic
#define mccm1_r0 tas1_mono_r0
#define mccm1_Q tas1_mono_q

/* Definition parameters for component 'm2'. */
#define mccm2_zmin -0.0375
#define mccm2_zmax 0.0375
#define mccm2_ymin -0.006
#define mccm2_ymax 0.006
#define mccm2_mosaich tas1_mono_mosaic
#define mccm2_mosaicv tas1_mono_mosaic
#define mccm2_r0 tas1_mono_r0
#define mccm2_Q tas1_mono_q

/* Definition parameters for component 'm3'. */
#define mccm3_zmin -0.0375
#define mccm3_zmax 0.0375
#define mccm3_ymin -0.006
#define mccm3_ymax 0.006
#define mccm3_mosaich tas1_mono_mosaic
#define mccm3_mosaicv tas1_mono_mosaic
#define mccm3_r0 tas1_mono_r0
#define mccm3_Q tas1_mono_q

/* Definition parameters for component 'm4'. */
#define mccm4_zmin -0.0375
#define mccm4_zmax 0.0375
#define mccm4_ymin -0.006
#define mccm4_ymax 0.006
#define mccm4_mosaich tas1_mono_mosaic
#define mccm4_mosaicv tas1_mono_mosaic
#define mccm4_r0 tas1_mono_r0
#define mccm4_Q tas1_mono_q

/* Definition parameters for component 'm5'. */
#define mccm5_zmin -0.0375
#define mccm5_zmax 0.0375
#define mccm5_ymin -0.006
#define mccm5_ymax 0.006
#define mccm5_mosaich tas1_mono_mosaic
#define mccm5_mosaicv tas1_mono_mosaic
#define mccm5_r0 tas1_mono_r0
#define mccm5_Q tas1_mono_q

/* Definition parameters for component 'm6'. */
#define mccm6_zmin -0.0375
#define mccm6_zmax 0.0375
#define mccm6_ymin -0.006
#define mccm6_ymax 0.006
#define mccm6_mosaich tas1_mono_mosaic
#define mccm6_mosaicv tas1_mono_mosaic
#define mccm6_r0 tas1_mono_r0
#define mccm6_Q tas1_mono_q

/* Definition parameters for component 'm7'. */
#define mccm7_zmin -0.0375
#define mccm7_zmax 0.0375
#define mccm7_ymin -0.006
#define mccm7_ymax 0.006
#define mccm7_mosaich tas1_mono_mosaic
#define mccm7_mosaicv tas1_mono_mosaic
#define mccm7_r0 tas1_mono_r0
#define mccm7_Q tas1_mono_q

/* Definition parameters for component 'slitMS1'. */
#define mccslitMS1_xmin -0.0105
#define mccslitMS1_xmax 0.0105
#define mccslitMS1_ymin -0.035
#define mccslitMS1_ymax 0.035

/* Definition parameters for component 'slitMS2'. */
#define mccslitMS2_xmin -0.0105
#define mccslitMS2_xmax 0.0105
#define mccslitMS2_ymin -0.035
#define mccslitMS2_ymax 0.035

/* Definition parameters for component 'c1'. */
#define mccc1_xmin -0.02
#define mccc1_xmax 0.02
#define mccc1_ymin -0.0375
#define mccc1_ymax 0.0375
#define mccc1_len 0.25
#define mccc1_divergence mcipC1

/* Definition parameters for component 'slitMS3'. */
#define mccslitMS3_radius 0.025

/* Definition parameters for component 'slitMS4'. */
#define mccslitMS4_radius 0.025

/* Definition parameters for component 'slitMS5'. */
#define mccslitMS5_radius 0.0275

/* Definition parameters for component 'emon1'. */
#define mccemon1_xmin -0.01
#define mccemon1_xmax 0.01
#define mccemon1_ymin -0.1
#define mccemon1_ymax 0.1
#define mccemon1_Emin 19.25
#define mccemon1_Emax 20.75
#define mccemon1_nchan 35
#define mccemon1_filename "sim/linup_1_1.emon"

/* Definition parameters for component 'slitS'. */
#define mccslitS_xmin -0.00525
#define mccslitS_xmax 0.00525
#define mccslitS_ymin -0.02025
#define mccslitS_ymax 0.02025

/* Definition parameters for component 'sng'. */
#define mccsng_xmin -0.025
#define mccsng_xmax 0.025
#define mccsng_ymin -0.0375
#define mccsng_ymax 0.0375

/* User component declarations. */

/* User declarations for component 'source'. */
#define mccompcurname "source"
#define radius mccsource_radius
#define dist mccsource_dist
#define xw mccsource_xw
#define yh mccsource_yh
#define E0 mccsource_E0
#define dE mccsource_dE
#define hdiv mccsource_hdiv
#define vdiv mccsource_vdiv
#define p_in mccsource_p_in
#line 36 "/data/lnslib/lib/mcstas/Source_flat.comp"
 double hdiv,vdiv;
 double p_in;
#line 1306 "linup-1.c"
#undef p_in
#undef vdiv
#undef hdiv
#undef dE
#undef E0
#undef yh
#undef xw
#undef dist
#undef radius
#undef mccompcurname

/* User declarations for component 'm0'. */
#define mccompcurname "m0"
#define zmin mccm0_zmin
#define zmax mccm0_zmax
#define ymin mccm0_ymin
#define ymax mccm0_ymax
#define mosaich mccm0_mosaich
#define mosaicv mccm0_mosaicv
#define r0 mccm0_r0
#define Q mccm0_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1330 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm1'. */
#define mccompcurname "m1"
#define zmin mccm1_zmin
#define zmax mccm1_zmax
#define ymin mccm1_ymin
#define ymax mccm1_ymax
#define mosaich mccm1_mosaich
#define mosaicv mccm1_mosaicv
#define r0 mccm1_r0
#define Q mccm1_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1353 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm2'. */
#define mccompcurname "m2"
#define zmin mccm2_zmin
#define zmax mccm2_zmax
#define ymin mccm2_ymin
#define ymax mccm2_ymax
#define mosaich mccm2_mosaich
#define mosaicv mccm2_mosaicv
#define r0 mccm2_r0
#define Q mccm2_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1376 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm3'. */
#define mccompcurname "m3"
#define zmin mccm3_zmin
#define zmax mccm3_zmax
#define ymin mccm3_ymin
#define ymax mccm3_ymax
#define mosaich mccm3_mosaich
#define mosaicv mccm3_mosaicv
#define r0 mccm3_r0
#define Q mccm3_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1399 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm4'. */
#define mccompcurname "m4"
#define zmin mccm4_zmin
#define zmax mccm4_zmax
#define ymin mccm4_ymin
#define ymax mccm4_ymax
#define mosaich mccm4_mosaich
#define mosaicv mccm4_mosaicv
#define r0 mccm4_r0
#define Q mccm4_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1422 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm5'. */
#define mccompcurname "m5"
#define zmin mccm5_zmin
#define zmax mccm5_zmax
#define ymin mccm5_ymin
#define ymax mccm5_ymax
#define mosaich mccm5_mosaich
#define mosaicv mccm5_mosaicv
#define r0 mccm5_r0
#define Q mccm5_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1445 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm6'. */
#define mccompcurname "m6"
#define zmin mccm6_zmin
#define zmax mccm6_zmax
#define ymin mccm6_ymin
#define ymax mccm6_ymax
#define mosaich mccm6_mosaich
#define mosaicv mccm6_mosaicv
#define r0 mccm6_r0
#define Q mccm6_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1468 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'm7'. */
#define mccompcurname "m7"
#define zmin mccm7_zmin
#define zmax mccm7_zmax
#define ymin mccm7_ymin
#define ymax mccm7_ymax
#define mosaich mccm7_mosaich
#define mosaicv mccm7_mosaicv
#define r0 mccm7_r0
#define Q mccm7_Q
#line 37 "/data/lnslib/lib/mcstas/Monochromator.comp"
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
#line 1491 "linup-1.c"
#undef Q
#undef r0
#undef mosaicv
#undef mosaich
#undef ymax
#undef ymin
#undef zmax
#undef zmin
#undef mccompcurname

/* User declarations for component 'c1'. */
#define mccompcurname "c1"
#define xmin mccc1_xmin
#define xmax mccc1_xmax
#define ymin mccc1_ymin
#define ymax mccc1_ymax
#define len mccc1_len
#define divergence mccc1_divergence
#define slope mccc1_slope
#line 36 "/data/lnslib/lib/mcstas/Soller.comp"
  double slope;
#line 1513 "linup-1.c"
#undef slope
#undef divergence
#undef len
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

/* User declarations for component 'emon1'. */
#define mccompcurname "emon1"
#define xmin mccemon1_xmin
#define xmax mccemon1_xmax
#define ymin mccemon1_ymin
#define ymax mccemon1_ymax
#define Emin mccemon1_Emin
#define Emax mccemon1_Emax
#define nchan mccemon1_nchan
#define filename mccemon1_filename
#define E_N mccemon1_E_N
#define E_p mccemon1_E_p
#define E_p2 mccemon1_E_p2
#line 40 "/data/lnslib/lib/mcstas/E_monitor.comp"
    int E_N[nchan];
    double E_p[nchan], E_p2[nchan];
#line 1539 "linup-1.c"
#undef E_p2
#undef E_p
#undef E_N
#undef filename
#undef nchan
#undef Emax
#undef Emin
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

/* User declarations for component 'sng'. */
#define mccompcurname "sng"
#define xmin mccsng_xmin
#define xmax mccsng_xmax
#define ymin mccsng_ymin
#define ymax mccsng_ymax
#define Nsum mccsng_Nsum
#define psum mccsng_psum
#define p2sum mccsng_p2sum
#line 36 "/data/lnslib/lib/mcstas/Monitor.comp"
    int Nsum;
    double psum, p2sum;
#line 1565 "linup-1.c"
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

Coords mcposaa1, mcposra1;
Rotation mcrotaa1, mcrotra1;
Coords mcposasource, mcposrsource;
Rotation mcrotasource, mcrotrsource;
Coords mcposaslit1, mcposrslit1;
Rotation mcrotaslit1, mcrotrslit1;
Coords mcposaslit2, mcposrslit2;
Rotation mcrotaslit2, mcrotrslit2;
Coords mcposaslit3, mcposrslit3;
Rotation mcrotaslit3, mcrotrslit3;
Coords mcposafocus_mono, mcposrfocus_mono;
Rotation mcrotafocus_mono, mcrotrfocus_mono;
Coords mcposam0, mcposrm0;
Rotation mcrotam0, mcrotrm0;
Coords mcposam1, mcposrm1;
Rotation mcrotam1, mcrotrm1;
Coords mcposam2, mcposrm2;
Rotation mcrotam2, mcrotrm2;
Coords mcposam3, mcposrm3;
Rotation mcrotam3, mcrotrm3;
Coords mcposam4, mcposrm4;
Rotation mcrotam4, mcrotrm4;
Coords mcposam5, mcposrm5;
Rotation mcrotam5, mcrotrm5;
Coords mcposam6, mcposrm6;
Rotation mcrotam6, mcrotrm6;
Coords mcposam7, mcposrm7;
Rotation mcrotam7, mcrotrm7;
Coords mcposaa2, mcposra2;
Rotation mcrotaa2, mcrotra2;
Coords mcposaslitMS1, mcposrslitMS1;
Rotation mcrotaslitMS1, mcrotrslitMS1;
Coords mcposaslitMS2, mcposrslitMS2;
Rotation mcrotaslitMS2, mcrotrslitMS2;
Coords mcposac1, mcposrc1;
Rotation mcrotac1, mcrotrc1;
Coords mcposaslitMS3, mcposrslitMS3;
Rotation mcrotaslitMS3, mcrotrslitMS3;
Coords mcposaslitMS4, mcposrslitMS4;
Rotation mcrotaslitMS4, mcrotrslitMS4;
Coords mcposaslitMS5, mcposrslitMS5;
Rotation mcrotaslitMS5, mcrotrslitMS5;
Coords mcposaemon1, mcposremon1;
Rotation mcrotaemon1, mcrotremon1;
Coords mcposaa3, mcposra3;
Rotation mcrotaa3, mcrotra3;
Coords mcposaslitS, mcposrslitS;
Rotation mcrotaslitS, mcrotrslitS;
Coords mcposasng, mcposrsng;
Rotation mcrotasng, mcrotrsng;

MCNUM mcnx, mcny, mcnz, mcnvx, mcnvy, mcnvz, mcnt, mcnsx, mcnsy, mcnsz, mcnp;

void mcinit(void) {
#define PHM mcipPHM
#define TTM mcipTTM
#define C1 mcipC1
#line 16 "linup-1.instr"
{
  double d = 0.0125;		/* 12.5 mm between slab centers. */
  double phi = 0.5443;		/* Rotation between adjacent slabs. */
  mpos0 = -3.5*d; mrot0 = -3.5*phi;
  mpos1 = -2.5*d; mrot1 = -2.5*phi;
  mpos2 = -1.5*d; mrot2 = -1.5*phi;
  mpos3 = -0.5*d; mrot3 = -0.5*phi;
  mpos4 =  0.5*d; mrot4 =  0.5*phi;
  mpos5 =  1.5*d; mrot5 =  1.5*phi;
  mpos6 =  2.5*d; mrot6 =  2.5*phi;
  mpos7 =  3.5*d; mrot7 =  3.5*phi;
}
#line 1645 "linup-1.c"
#undef C1
#undef TTM
#undef PHM
  /* Component initializations. */
  /* Initializations for component a1. */


  /* Initializations for component source. */

#define mccompcurname "source"
#define radius mccsource_radius
#define dist mccsource_dist
#define xw mccsource_xw
#define yh mccsource_yh
#define E0 mccsource_E0
#define dE mccsource_dE
#define hdiv mccsource_hdiv
#define vdiv mccsource_vdiv
#define p_in mccsource_p_in
#line 40 "/data/lnslib/lib/mcstas/Source_flat.comp"
{
 hdiv = atan(xw/(2.0*dist));
 vdiv = atan(yh/(2.0*dist));
 p_in = (4*hdiv*vdiv)/(4*PI); /* Small angle approx. */
}
#line 1671 "linup-1.c"
#undef p_in
#undef vdiv
#undef hdiv
#undef dE
#undef E0
#undef yh
#undef xw
#undef dist
#undef radius
#undef mccompcurname

  /* Initializations for component slit1. */


  /* Initializations for component slit2. */


  /* Initializations for component slit3. */


  /* Initializations for component focus_mono. */


  /* Initializations for component m0. */


  /* Initializations for component m1. */


  /* Initializations for component m2. */


  /* Initializations for component m3. */


  /* Initializations for component m4. */


  /* Initializations for component m5. */


  /* Initializations for component m6. */


  /* Initializations for component m7. */


  /* Initializations for component a2. */


  /* Initializations for component slitMS1. */


  /* Initializations for component slitMS2. */


  /* Initializations for component c1. */

#define mccompcurname "c1"
#define xmin mccc1_xmin
#define xmax mccc1_xmax
#define ymin mccc1_ymin
#define ymax mccc1_ymax
#define len mccc1_len
#define divergence mccc1_divergence
#define slope mccc1_slope
#line 39 "/data/lnslib/lib/mcstas/Soller.comp"
{
  slope = tan(MIN2RAD*divergence);
}
#line 1742 "linup-1.c"
#undef slope
#undef divergence
#undef len
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

  /* Initializations for component slitMS3. */


  /* Initializations for component slitMS4. */


  /* Initializations for component slitMS5. */


  /* Initializations for component emon1. */

#define mccompcurname "emon1"
#define xmin mccemon1_xmin
#define xmax mccemon1_xmax
#define ymin mccemon1_ymin
#define ymax mccemon1_ymax
#define Emin mccemon1_Emin
#define Emax mccemon1_Emax
#define nchan mccemon1_nchan
#define filename mccemon1_filename
#define E_N mccemon1_E_N
#define E_p mccemon1_E_p
#define E_p2 mccemon1_E_p2
#line 44 "/data/lnslib/lib/mcstas/E_monitor.comp"
{
    int i;

    for (i=0; i<nchan; i++)
    {
      E_N[i] = 0;
      E_p[i] = 0;
      E_p2[i] = 0;
    }
}
#line 1786 "linup-1.c"
#undef E_p2
#undef E_p
#undef E_N
#undef filename
#undef nchan
#undef Emax
#undef Emin
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

  /* Initializations for component a3. */


  /* Initializations for component slitS. */


  /* Initializations for component sng. */

#define mccompcurname "sng"
#define xmin mccsng_xmin
#define xmax mccsng_xmax
#define ymin mccsng_ymin
#define ymax mccsng_ymax
#define Nsum mccsng_Nsum
#define psum mccsng_psum
#define p2sum mccsng_p2sum
#line 40 "/data/lnslib/lib/mcstas/Monitor.comp"
{
    psum = 0;
    p2sum = 0;
    Nsum = 0;
}
#line 1822 "linup-1.c"
#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 a1. */
    rot_set_rotation(mcrotaa1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_copy(mcrotra1, mcrotaa1);
    mcposaa1 = coords_set(0, 0, 0);
    mctc1 = coords_neg(mcposaa1);
    mcposra1 = rot_apply(mcrotaa1, mctc1);
    mcDEBUG_COMPONENT("a1", mcposaa1, mcrotaa1)
    /* Component source. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa1, mcrotasource);
    rot_transpose(mcrotaa1, mctr1);
    rot_mul(mcrotasource, mctr1, mcrotrsource);
    mctc1 = coords_set(0, 0, 0);
    rot_transpose(mcrotaa1, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposasource = coords_add(mcposaa1, mctc2);
    mctc1 = coords_sub(mcposaa1, mcposasource);
    mcposrsource = rot_apply(mcrotasource, mctc1);
    mcDEBUG_COMPONENT("source", mcposasource, mcrotasource)
    /* Component slit1. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa1, mcrotaslit1);
    rot_transpose(mcrotasource, mctr1);
    rot_mul(mcrotaslit1, mctr1, mcrotrslit1);
    mctc1 = coords_set(0, 0, 1.1215);
    rot_transpose(mcrotaa1, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslit1 = coords_add(mcposaa1, mctc2);
    mctc1 = coords_sub(mcposasource, 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, mcrotaa1, mcrotaslit2);
    rot_transpose(mcrotaslit1, mctr1);
    rot_mul(mcrotaslit2, mctr1, mcrotrslit2);
    mctc1 = coords_set(0, 0, 1.9);
    rot_transpose(mcrotaa1, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslit2 = coords_add(mcposaa1, mctc2);
    mctc1 = coords_sub(mcposaslit1, mcposaslit2);
    mcposrslit2 = rot_apply(mcrotaslit2, mctc1);
    mcDEBUG_COMPONENT("slit2", mcposaslit2, mcrotaslit2)
    /* Component slit3. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa1, mcrotaslit3);
    rot_transpose(mcrotaslit2, mctr1);
    rot_mul(mcrotaslit3, mctr1, mcrotrslit3);
    mctc1 = coords_set(0, 0, 3.288);
    rot_transpose(mcrotaa1, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslit3 = coords_add(mcposaa1, mctc2);
    mctc1 = coords_sub(mcposaslit2, mcposaslit3);
    mcposrslit3 = rot_apply(mcrotaslit3, mctc1);
    mcDEBUG_COMPONENT("slit3", mcposaslit3, mcrotaslit3)
    /* Component focus_mono. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipPHM)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa1, mcrotafocus_mono);
    rot_transpose(mcrotaslit3, mctr1);
    rot_mul(mcrotafocus_mono, mctr1, mcrotrfocus_mono);
    mctc1 = coords_set(0, 0, 3.56);
    rot_transpose(mcrotaa1, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposafocus_mono = coords_add(mcposaa1, mctc2);
    mctc1 = coords_sub(mcposaslit3, mcposafocus_mono);
    mcposrfocus_mono = rot_apply(mcrotafocus_mono, mctc1);
    mcDEBUG_COMPONENT("focus_mono", mcposafocus_mono, mcrotafocus_mono)
    /* Component m0. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot0)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam0);
    rot_transpose(mcrotafocus_mono, mctr1);
    rot_mul(mcrotam0, mctr1, mcrotrm0);
    mctc1 = coords_set(0, mpos0, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam0 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposafocus_mono, mcposam0);
    mcposrm0 = rot_apply(mcrotam0, mctc1);
    mcDEBUG_COMPONENT("m0", mcposam0, mcrotam0)
    /* Component m1. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot1)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam1);
    rot_transpose(mcrotam0, mctr1);
    rot_mul(mcrotam1, mctr1, mcrotrm1);
    mctc1 = coords_set(0, mpos1, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam1 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam0, mcposam1);
    mcposrm1 = rot_apply(mcrotam1, mctc1);
    mcDEBUG_COMPONENT("m1", mcposam1, mcrotam1)
    /* Component m2. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot2)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam2);
    rot_transpose(mcrotam1, mctr1);
    rot_mul(mcrotam2, mctr1, mcrotrm2);
    mctc1 = coords_set(0, mpos2, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam2 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam1, mcposam2);
    mcposrm2 = rot_apply(mcrotam2, mctc1);
    mcDEBUG_COMPONENT("m2", mcposam2, mcrotam2)
    /* Component m3. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot3)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam3);
    rot_transpose(mcrotam2, mctr1);
    rot_mul(mcrotam3, mctr1, mcrotrm3);
    mctc1 = coords_set(0, mpos3, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam3 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam2, mcposam3);
    mcposrm3 = rot_apply(mcrotam3, mctc1);
    mcDEBUG_COMPONENT("m3", mcposam3, mcrotam3)
    /* Component m4. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot4)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam4);
    rot_transpose(mcrotam3, mctr1);
    rot_mul(mcrotam4, mctr1, mcrotrm4);
    mctc1 = coords_set(0, mpos4, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam4 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam3, mcposam4);
    mcposrm4 = rot_apply(mcrotam4, mctc1);
    mcDEBUG_COMPONENT("m4", mcposam4, mcrotam4)
    /* Component m5. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot5)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam5);
    rot_transpose(mcrotam4, mctr1);
    rot_mul(mcrotam5, mctr1, mcrotrm5);
    mctc1 = coords_set(0, mpos5, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam5 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam4, mcposam5);
    mcposrm5 = rot_apply(mcrotam5, mctc1);
    mcDEBUG_COMPONENT("m5", mcposam5, mcrotam5)
    /* Component m6. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot6)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam6);
    rot_transpose(mcrotam5, mctr1);
    rot_mul(mcrotam6, mctr1, mcrotrm6);
    mctc1 = coords_set(0, mpos6, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam6 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam5, mcposam6);
    mcposrm6 = rot_apply(mcrotam6, mctc1);
    mcDEBUG_COMPONENT("m6", mcposam6, mcrotam6)
    /* Component m7. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (mrot7)*DEG2RAD);
    rot_mul(mctr1, mcrotafocus_mono, mcrotam7);
    rot_transpose(mcrotam6, mctr1);
    rot_mul(mcrotam7, mctr1, mcrotrm7);
    mctc1 = coords_set(0, mpos7, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposam7 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam6, mcposam7);
    mcposrm7 = rot_apply(mcrotam7, mctc1);
    mcDEBUG_COMPONENT("m7", mcposam7, mcrotam7)
    /* Component a2. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (mcipTTM)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa1, mcrotaa2);
    rot_transpose(mcrotam7, mctr1);
    rot_mul(mcrotaa2, mctr1, mcrotra2);
    mctc1 = coords_set(0, 0, 0);
    rot_transpose(mcrotafocus_mono, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaa2 = coords_add(mcposafocus_mono, mctc2);
    mctc1 = coords_sub(mcposam7, mcposaa2);
    mcposra2 = rot_apply(mcrotaa2, mctc1);
    mcDEBUG_COMPONENT("a2", mcposaa2, mcrotaa2)
    /* Component slitMS1. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaslitMS1);
    rot_transpose(mcrotaa2, mctr1);
    rot_mul(mcrotaslitMS1, mctr1, mcrotrslitMS1);
    mctc1 = coords_set(0, 0, 0.565);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslitMS1 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaa2, mcposaslitMS1);
    mcposrslitMS1 = rot_apply(mcrotaslitMS1, mctc1);
    mcDEBUG_COMPONENT("slitMS1", mcposaslitMS1, mcrotaslitMS1)
    /* Component slitMS2. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaslitMS2);
    rot_transpose(mcrotaslitMS1, mctr1);
    rot_mul(mcrotaslitMS2, mctr1, mcrotrslitMS2);
    mctc1 = coords_set(0, 0, 0.855);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslitMS2 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaslitMS1, mcposaslitMS2);
    mcposrslitMS2 = rot_apply(mcrotaslitMS2, mctc1);
    mcDEBUG_COMPONENT("slitMS2", mcposaslitMS2, mcrotaslitMS2)
    /* Component c1. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotac1);
    rot_transpose(mcrotaslitMS2, mctr1);
    rot_mul(mcrotac1, mctr1, mcrotrc1);
    mctc1 = coords_set(0, 0, 0.87);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposac1 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaslitMS2, mcposac1);
    mcposrc1 = rot_apply(mcrotac1, mctc1);
    mcDEBUG_COMPONENT("c1", mcposac1, mcrotac1)
    /* Component slitMS3. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaslitMS3);
    rot_transpose(mcrotac1, mctr1);
    rot_mul(mcrotaslitMS3, mctr1, mcrotrslitMS3);
    mctc1 = coords_set(0, 0, 1.13);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslitMS3 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposac1, mcposaslitMS3);
    mcposrslitMS3 = rot_apply(mcrotaslitMS3, mctc1);
    mcDEBUG_COMPONENT("slitMS3", mcposaslitMS3, mcrotaslitMS3)
    /* Component slitMS4. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaslitMS4);
    rot_transpose(mcrotaslitMS3, mctr1);
    rot_mul(mcrotaslitMS4, mctr1, mcrotrslitMS4);
    mctc1 = coords_set(0, 0, 1.18);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslitMS4 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaslitMS3, mcposaslitMS4);
    mcposrslitMS4 = rot_apply(mcrotaslitMS4, mctc1);
    mcDEBUG_COMPONENT("slitMS4", mcposaslitMS4, mcrotaslitMS4)
    /* Component slitMS5. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaslitMS5);
    rot_transpose(mcrotaslitMS4, mctr1);
    rot_mul(mcrotaslitMS5, mctr1, mcrotrslitMS5);
    mctc1 = coords_set(0, 0, 1.23);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslitMS5 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaslitMS4, mcposaslitMS5);
    mcposrslitMS5 = rot_apply(mcrotaslitMS5, mctc1);
    mcDEBUG_COMPONENT("slitMS5", mcposaslitMS5, mcrotaslitMS5)
    /* Component emon1. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaemon1);
    rot_transpose(mcrotaslitMS5, mctr1);
    rot_mul(mcrotaemon1, mctr1, mcrotremon1);
    mctc1 = coords_set(0, 0, 1.5);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaemon1 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaslitMS5, mcposaemon1);
    mcposremon1 = rot_apply(mcrotaemon1, mctc1);
    mcDEBUG_COMPONENT("emon1", mcposaemon1, mcrotaemon1)
    /* Component a3. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa2, mcrotaa3);
    rot_transpose(mcrotaemon1, mctr1);
    rot_mul(mcrotaa3, mctr1, mcrotra3);
    mctc1 = coords_set(0, 0, 1.565);
    rot_transpose(mcrotaa2, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaa3 = coords_add(mcposaa2, mctc2);
    mctc1 = coords_sub(mcposaemon1, mcposaa3);
    mcposra3 = rot_apply(mcrotaa3, mctc1);
    mcDEBUG_COMPONENT("a3", mcposaa3, mcrotaa3)
    /* Component slitS. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa3, mcrotaslitS);
    rot_transpose(mcrotaa3, mctr1);
    rot_mul(mcrotaslitS, mctr1, mcrotrslitS);
    mctc1 = coords_set(0, 0, 0);
    rot_transpose(mcrotaa3, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposaslitS = coords_add(mcposaa3, mctc2);
    mctc1 = coords_sub(mcposaa3, mcposaslitS);
    mcposrslitS = rot_apply(mcrotaslitS, mctc1);
    mcDEBUG_COMPONENT("slitS", mcposaslitS, mcrotaslitS)
    /* Component sng. */
    rot_set_rotation(mctr1, (0)*DEG2RAD, (0)*DEG2RAD, (0)*DEG2RAD);
    rot_mul(mctr1, mcrotaa3, mcrotasng);
    rot_transpose(mcrotaslitS, mctr1);
    rot_mul(mcrotasng, mctr1, mcrotrsng);
    mctc1 = coords_set(0, 0, 0.02);
    rot_transpose(mcrotaa3, mctr1);
    mctc2 = rot_apply(mctr1, mctc1);
    mcposasng = coords_add(mcposaa3, mctc2);
    mctc1 = coords_sub(mcposaslitS, mcposasng);
    mcposrsng = rot_apply(mcrotasng, mctc1);
    mcDEBUG_COMPONENT("sng", mcposasng, mcrotasng)
    if(mcdotrace) mcdisplay();
    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 mcnlsx = mcnsx;
  MCNUM mcnlsy = mcnsy;
  MCNUM mcnlsz = mcnsz;
  MCNUM mcnlp = mcnp;

  mcDEBUG_ENTER()
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
  /* Component a1. */
  mcDEBUG_COMP("a1")
  mccoordschange(mcposra1, mcrotra1,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define mccompcurname "a1"
#undef mccompcurname
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)

  /* Component source. */
  mcDEBUG_COMP("source")
  mccoordschange(mcposrsource, mcrotrsource,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "source"
#define radius mccsource_radius
#define dist mccsource_dist
#define xw mccsource_xw
#define yh mccsource_yh
#define E0 mccsource_E0
#define dE mccsource_dE
#define hdiv mccsource_hdiv
#define vdiv mccsource_vdiv
#define p_in mccsource_p_in
#line 46 "/data/lnslib/lib/mcstas/Source_flat.comp"
{
 double theta0,phi0,chi,theta,phi,E,v,r;

 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();

 E=E0+dE*randpm1();                  /* Assume linear distribution */
 v=sqrt(E)*SE2V;

 vz=v*cos(phi)*cos(theta);
 vy=v*sin(phi);
 vx=v*cos(phi)*sin(theta); 
}
#line 2218 "linup-1.c"
#undef p_in
#undef vdiv
#undef hdiv
#undef dE
#undef E0
#undef yh
#undef xw
#undef dist
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slit1. */
  mcDEBUG_COMP("slit1")
  mccoordschange(mcposrslit1, mcrotrslit1,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slit1"
#define xmin mccslit1_xmin
#define xmax mccslit1_xmax
#define ymin mccslit1_ymin
#define ymax mccslit1_ymax
#line 28 "/data/lnslib/lib/mcstas/Slit.comp"
{
    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax)
      ABSORB;
}
#line 2269 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slit2. */
  mcDEBUG_COMP("slit2")
  mccoordschange(mcposrslit2, mcrotrslit2,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slit2"
#define xmin mccslit2_xmin
#define xmax mccslit2_xmax
#define ymin mccslit2_ymin
#define ymax mccslit2_ymax
#line 28 "/data/lnslib/lib/mcstas/Slit.comp"
{
    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax)
      ABSORB;
}
#line 2315 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slit3. */
  mcDEBUG_COMP("slit3")
  mccoordschange(mcposrslit3, mcrotrslit3,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slit3"
#define xmin mccslit3_xmin
#define xmax mccslit3_xmax
#define ymin mccslit3_ymin
#define ymax mccslit3_ymax
#line 28 "/data/lnslib/lib/mcstas/Slit.comp"
{
    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax)
      ABSORB;
}
#line 2361 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component focus_mono. */
  mcDEBUG_COMP("focus_mono")
  mccoordschange(mcposrfocus_mono, mcrotrfocus_mono,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define mccompcurname "focus_mono"
#undef mccompcurname
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)

  /* Component m0. */
  mcDEBUG_COMP("m0")
  mccoordschange(mcposrm0, mcrotrm0,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m0"
#define zmin mccm0_zmin
#define zmax mccm0_zmax
#define ymin mccm0_ymin
#define ymax mccm0_ymax
#define mosaich mccm0_mosaich
#define mosaicv mccm0_mosaicv
#define r0 mccm0_r0
#define Q mccm0_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 2470 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m1. */
  mcDEBUG_COMP("m1")
  mccoordschange(mcposrm1, mcrotrm1,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m1"
#define zmin mccm1_zmin
#define zmax mccm1_zmax
#define ymin mccm1_ymin
#define ymax mccm1_ymax
#define mosaich mccm1_mosaich
#define mosaicv mccm1_mosaicv
#define r0 mccm1_r0
#define Q mccm1_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 2572 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m2. */
  mcDEBUG_COMP("m2")
  mccoordschange(mcposrm2, mcrotrm2,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m2"
#define zmin mccm2_zmin
#define zmax mccm2_zmax
#define ymin mccm2_ymin
#define ymax mccm2_ymax
#define mosaich mccm2_mosaich
#define mosaicv mccm2_mosaicv
#define r0 mccm2_r0
#define Q mccm2_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 2674 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m3. */
  mcDEBUG_COMP("m3")
  mccoordschange(mcposrm3, mcrotrm3,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m3"
#define zmin mccm3_zmin
#define zmax mccm3_zmax
#define ymin mccm3_ymin
#define ymax mccm3_ymax
#define mosaich mccm3_mosaich
#define mosaicv mccm3_mosaicv
#define r0 mccm3_r0
#define Q mccm3_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 2776 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m4. */
  mcDEBUG_COMP("m4")
  mccoordschange(mcposrm4, mcrotrm4,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m4"
#define zmin mccm4_zmin
#define zmax mccm4_zmax
#define ymin mccm4_ymin
#define ymax mccm4_ymax
#define mosaich mccm4_mosaich
#define mosaicv mccm4_mosaicv
#define r0 mccm4_r0
#define Q mccm4_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 2878 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m5. */
  mcDEBUG_COMP("m5")
  mccoordschange(mcposrm5, mcrotrm5,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m5"
#define zmin mccm5_zmin
#define zmax mccm5_zmax
#define ymin mccm5_ymin
#define ymax mccm5_ymax
#define mosaich mccm5_mosaich
#define mosaicv mccm5_mosaicv
#define r0 mccm5_r0
#define Q mccm5_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 2980 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m6. */
  mcDEBUG_COMP("m6")
  mccoordschange(mcposrm6, mcrotrm6,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m6"
#define zmin mccm6_zmin
#define zmax mccm6_zmax
#define ymin mccm6_ymin
#define ymax mccm6_ymax
#define mosaich mccm6_mosaich
#define mosaicv mccm6_mosaicv
#define r0 mccm6_r0
#define Q mccm6_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 3082 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component m7. */
  mcDEBUG_COMP("m7")
  mccoordschange(mcposrm7, mcrotrm7,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "m7"
#define zmin mccm7_zmin
#define zmax mccm7_zmax
#define ymin mccm7_ymin
#define ymax mccm7_ymax
#define mosaich mccm7_mosaich
#define mosaicv mccm7_mosaicv
#define r0 mccm7_r0
#define Q mccm7_Q
#line 41 "/data/lnslib/lib/mcstas/Monochromator.comp"
{
    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;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
}
#line 3184 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component a2. */
  mcDEBUG_COMP("a2")
  mccoordschange(mcposra2, mcrotra2,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define mccompcurname "a2"
#undef mccompcurname
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)

  /* Component slitMS1. */
  mcDEBUG_COMP("slitMS1")
  mccoordschange(mcposrslitMS1, mcrotrslitMS1,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slitMS1"
#define xmin mccslitMS1_xmin
#define xmax mccslitMS1_xmax
#define ymin mccslitMS1_ymin
#define ymax mccslitMS1_ymax
#line 28 "/data/lnslib/lib/mcstas/Slit.comp"
{
    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax)
      ABSORB;
}
#line 3245 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slitMS2. */
  mcDEBUG_COMP("slitMS2")
  mccoordschange(mcposrslitMS2, mcrotrslitMS2,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slitMS2"
#define xmin mccslitMS2_xmin
#define xmax mccslitMS2_xmax
#define ymin mccslitMS2_ymin
#define ymax mccslitMS2_ymax
#line 28 "/data/lnslib/lib/mcstas/Slit.comp"
{
    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax)
      ABSORB;
}
#line 3291 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component c1. */
  mcDEBUG_COMP("c1")
  mccoordschange(mcposrc1, mcrotrc1,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "c1"
#define xmin mccc1_xmin
#define xmax mccc1_xmax
#define ymin mccc1_ymin
#define ymax mccc1_ymax
#define len mccc1_len
#define divergence mccc1_divergence
#define slope mccc1_slope
#line 43 "/data/lnslib/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;
    }
}
#line 3355 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slitMS3. */
  mcDEBUG_COMP("slitMS3")
  mccoordschange(mcposrslitMS3, mcrotrslitMS3,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slitMS3"
#define radius mccslitMS3_radius
#line 22 "/data/lnslib/lib/mcstas/Circular_slit.comp"
{
  PROP_Z0;
  if(x*x + y*y > radius*radius)
    ABSORB;
}
#line 3401 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slitMS4. */
  mcDEBUG_COMP("slitMS4")
  mccoordschange(mcposrslitMS4, mcrotrslitMS4,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slitMS4"
#define radius mccslitMS4_radius
#line 22 "/data/lnslib/lib/mcstas/Circular_slit.comp"
{
  PROP_Z0;
  if(x*x + y*y > radius*radius)
    ABSORB;
}
#line 3441 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component slitMS5. */
  mcDEBUG_COMP("slitMS5")
  mccoordschange(mcposrslitMS5, mcrotrslitMS5,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slitMS5"
#define radius mccslitMS5_radius
#line 22 "/data/lnslib/lib/mcstas/Circular_slit.comp"
{
  PROP_Z0;
  if(x*x + y*y > radius*radius)
    ABSORB;
}
#line 3481 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component emon1. */
  mcDEBUG_COMP("emon1")
  mccoordschange(mcposremon1, mcrotremon1,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "emon1"
#define xmin mccemon1_xmin
#define xmax mccemon1_xmax
#define ymin mccemon1_ymin
#define ymax mccemon1_ymax
#define Emin mccemon1_Emin
#define Emax mccemon1_Emax
#define nchan mccemon1_nchan
#define filename mccemon1_filename
#define E_N mccemon1_E_N
#define E_p mccemon1_E_p
#define E_p2 mccemon1_E_p2
#line 55 "/data/lnslib/lib/mcstas/E_monitor.comp"
{
    int i;
    double E;

    PROP_Z0;
    if (x>xmin && x<xmax && y>ymin && y<ymax)
    {
      E = VS2E*(vx*vx + vy*vy + vz*vz);

      i = floor((E-Emin)*nchan/(Emax-Emin));
      if(i >= 0 && i < nchan)
      {
        E_N[i]++;
        E_p[i] += p;
        E_p2[i] += p*p;
      }
    }
}
#line 3544 "linup-1.c"
#undef E_p2
#undef E_p
#undef E_N
#undef filename
#undef nchan
#undef Emax
#undef Emin
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component a3. */
  mcDEBUG_COMP("a3")
  mccoordschange(mcposra3, mcrotra3,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define mccompcurname "a3"
#undef mccompcurname
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)

  /* Component slitS. */
  mcDEBUG_COMP("slitS")
  mccoordschange(mcposrslitS, mcrotrslitS,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "slitS"
#define xmin mccslitS_xmin
#define xmax mccslitS_xmax
#define ymin mccslitS_ymin
#define ymax mccslitS_ymax
#line 28 "/data/lnslib/lib/mcstas/Slit.comp"
{
    PROP_Z0;
    if (x<xmin || x>xmax || y<ymin || y>ymax)
      ABSORB;
}
#line 3608 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

  /* Component sng. */
  mcDEBUG_COMP("sng")
  mccoordschange(mcposrsng, mcrotrsng,
    &mcnlx, &mcnly, &mcnlz,
    &mcnlvx, &mcnlvy, &mcnlvz,
    &mcnlt, &mcnlsx, &mcnlsy);
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
#define x mcnlx
#define y mcnly
#define z mcnlz
#define vx mcnlvx
#define vy mcnlvy
#define vz mcnlvz
#define t mcnlt
#define s1 mcnlsx
#define s2 mcnlsy
#define p mcnlp
#define mccompcurname "sng"
#define xmin mccsng_xmin
#define xmax mccsng_xmax
#define ymin mccsng_ymin
#define ymax mccsng_ymax
#define Nsum mccsng_Nsum
#define psum mccsng_psum
#define p2sum mccsng_p2sum
#line 46 "/data/lnslib/lib/mcstas/Monitor.comp"
{
    PROP_Z0;
    if (x>xmin && x<xmax && y>ymin && y<ymax)
    {
      Nsum++;
      psum += p;
      p2sum += p*p;
    }
}
#line 3661 "linup-1.c"
#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,mcnlsx,mcnlsy, mcnlp)

 mcabsorb:
  mcDEBUG_LEAVE()
  mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz,mcnlt,mcnlsx,mcnlsy, mcnlp)
  /* Copy neutron state to global variables. */
  mcnx = mcnlx;
  mcny = mcnly;
  mcnz = mcnlz;
  mcnvx = mcnlvx;
  mcnvy = mcnlvy;
  mcnvz = mcnlvz;
  mcnt = mcnlt;
  mcnsx = mcnlsx;
  mcnsy = mcnlsy;
  mcnsz = mcnlsz;
  mcnp = mcnlp;
}

void mcfinally(void) {
  /* User component FINALLY code. */

  /* User FINALLY code for component 'emon1'. */
#define mccompcurname "emon1"
#define xmin mccemon1_xmin
#define xmax mccemon1_xmax
#define ymin mccemon1_ymin
#define ymax mccemon1_ymax
#define Emin mccemon1_Emin
#define Emax mccemon1_Emax
#define nchan mccemon1_nchan
#define filename mccemon1_filename
#define E_N mccemon1_E_N
#define E_p mccemon1_E_p
#define E_p2 mccemon1_E_p2
#line 74 "/data/lnslib/lib/mcstas/E_monitor.comp"
{
    int i;
    FILE *outfile;

    outfile=fopen(filename,"w");
    if(outfile == NULL)
      fprintf(stderr, "Error: could not open output file '%s'\n", filename);
    else
    {
      for (i=0; i<nchan; i++)
        fprintf(outfile,"%g %d %g %g\n", 
                Emin + (double)i/nchan*(Emax-Emin), E_N[i], E_p[i], E_p2[i]);
      fclose(outfile);
    }
}
#line 3731 "linup-1.c"
#undef E_p2
#undef E_p
#undef E_N
#undef filename
#undef nchan
#undef Emax
#undef Emin
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

  /* User FINALLY code for component 'sng'. */
#define mccompcurname "sng"
#define xmin mccsng_xmin
#define xmax mccsng_xmax
#define ymin mccsng_ymin
#define ymax mccsng_ymax
#define Nsum mccsng_Nsum
#define psum mccsng_psum
#define p2sum mccsng_p2sum
#line 56 "/data/lnslib/lib/mcstas/Monitor.comp"
{
    DETECTOR_OUT(Nsum, psum, p2sum);
}
#line 3758 "linup-1.c"
#undef p2sum
#undef psum
#undef Nsum
#undef ymax
#undef ymin
#undef xmax
#undef xmin
#undef mccompcurname

}
#define magnify mcdis_magnify
#define line mcdis_line
#define multiline mcdis_multiline
#define circle mcdis_circle
void mcdisplay(void) {
  printf("MCDISPLAY: start\n");
  /* Component MCDISPLAY code. */

  printf("MCDISPLAY: end\n");
}
#undef magnify
#undef line
#undef multiline
#undef circle



More information about the mcstas-users mailing list