/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2008, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: FermiChopper.
*
* %Identification
*
* Written by: M. Poehlmann, C. Carbogno, H. Schober, E. Farhi
* Date: May 2002
* Origin: ILL Grenoble / TU Muenchen
* Modified by: K.Lieutenant, June 2005: added phase parameter. Comp validation.
* Modified by: EF, Nov 2005: completely rewrote comp.
* Modified by: EF, Oct 2008: fix chopper orientation
* Modified by: EF, Mar 2009: improved intersection algorithm
*
* Fermi Chopper with rotating frame.
*
* %D
* Models a fermi chopper with optional supermirror coated blades
* supermirror facilities may be disabled by setting m = 0, R0=0
* Slit packages are straight. Chopper slits are separated by an infinitely
* thin absorbing material. The effective transmission (resulting from fraction
* of the transparent material and its transmission) may be specified.
* The chopper slit package width may be specified through the total width 'xwidth'
* of the full package or the width 'w' of each single slit. The other parameter
* is calculated by: xwidth = nslit*w. The slit package may be made curved and use
* super-mirror coating.
*
* Example:
* FermiChopper(phase=-50.0, radius=0.04, nu=100, yheight=0.08, w=0.00022475, nslit=200.0, R0=0.0, Qc=0.02176, alpha=2.33, m=0.0, length=0.012, eff=0.95)
*
* %VALIDATION
* Apr 2005: extensive external test, most problems solved (cf. 'Bugs')
* Validated by: K. Lieutenant, E. Farhi
*
* limitations:
* no absorbing blade width used
*
* %Parameters
* INPUT PARAMETERS:
*
* Geometrical chopper constants:
* radius: [m] chopper cylinder radius
* nslit: [1] number of chopper slits
* length: [m] channel length of the Fermi chopper
* w: [m] width of one chopper slit
* xwidth: [m] optional total width of slit package
* yheight: [m] height of slit package
* nu: [Hz] chopper frequency. Omega=2*PI*nu in rad/s, nu*60 in rpm. Positive value corresponds to a counter-clockwise rotation around y.
* eff: [1] efficiency = transmission x fraction of transparent material
* verbose: [1] set to 1,2 or 3 gives debugging information
* curvature: [m-1] Curvature of slits (1/radius of curvature).
*
* Supermirror constants:
* m: [1] m-value of material. Zero means completely absorbing.
* alpha: [AA] slope of reflectivity
* Qc: [AA-1] critical scattering vector
* W: [AA-1] width of supermirror cut-off
* R0: [1] low-angle reflectivity
*
* Constants to reset time of flight:
* zero_time: [1] set time to zero: 0=no, 1=once per half cycle, 2=auto adjust phase
* phase: [deg] chopper phase at t=0
* delay: [s] sets phase so that transmision is centered on 'delay'
*
* OUTPUT PARAMETERS:
* FCVars : [] structure
*
* %L
* Vitess_ChopperFermi component by
* G. Zsigmond, imported from Vitess by K. Lieutenant.
*
* %End
*****************************************************************************/
/* NOTE:
* The initial component version (McStas version <= 1.12) was written
* with an inverted coordinate frame orientation. This corresponds
* to inverting the frequency and phase sign.
*/
DEFINE COMPONENT FermiChopper
DEFINITION PARAMETERS ()
SETTING PARAMETERS (phase=0, radius=0.04, nu=100,
w=0.00022475, nslit=200, R0=0.0,
Qc=0.02176, alpha=2.33, m=0.0, W=2e-3, length=0.012, eff=0.95,
zero_time=0, xwidth=0, verbose=0, yheight=0.08,
curvature=0,delay=0)
OUTPUT PARAMETERS()
SHARE
%{
%include "ref-lib"
#ifndef FermiChopper_TimeAccuracy
#define FermiChopper_TimeAccuracy 1e-9
#define FermiChopper_MAXITER 100
/* Definition of internal variable structure: all counters */
struct FermiChopper_struct
{
double omega; /* chopper rotation */
double ph0; /* chopper rotation */
double t0; /* chopper rotation */
double C_slit; /* slit curvature radius in [m] */
double L_slit; /* slit package length [m] */
double sum_t;
double sum_v;
double sum_N;
double sum_N_pass;
/* events */
long absorb_alreadyinside;
long absorb_topbottom;
long absorb_cylentrance;
long absorb_sideentrance;
long absorb_notreachentrance;
long absorb_packentrance;
long absorb_slitcoating;
long warn_notreachslitwall;
long absorb_exitslitpack;
long absorb_maxiterations;
long absorb_wrongdirection;
long absorb_nocontrol;
long absorb_cylexit;
long warn_notreachslitoutput;
char compcurname[256];
};
/*****************************************************************************
* FC_zrot: returns Z' in rotating frame, from X,Z and t,omega,ph0
****************************************************************************/
#pragma acc routine seq
double FC_zrot(double X, double Z, double T, struct FermiChopper_struct FCs)
{
double omega =FCs.omega;
double ph0 =FCs.ph0;
return( Z*cos(omega*T+ph0)-X*sin(omega*T+ph0) );
}
/*****************************************************************************
* FC_xrot: returns X' in rotating frame, from X,Z and omega,t,ph0
* additional coordinate shift in case of curved slits
****************************************************************************/
#pragma acc routine seq
double FC_xrot(double X, double Z, double T, struct FermiChopper_struct FCs)
{
double omega =FCs.omega;
double ph0 =FCs.ph0;
double C_slit=FCs.C_slit;
double ret, tmp;
ret = X*cos(omega*T+ph0)+Z*sin(omega*T+ph0);
if (C_slit) {
tmp = fabs(FC_zrot(X, Z, T, FCs));
if (tmp < FCs.L_slit/2) {
tmp = (FCs.L_slit/2 - tmp)*C_slit;
ret += (1-sqrt(1-tmp*tmp))/C_slit;
}
}
return( ret );
}
/*****************************************************************************
* FC_xzrot_dt(x,z,vx,vz, t,dt, type='x' or 'z', FCs)
* returns X' or Z' in rotating frame, from X,Z and t,omega,ph0
* taking into account propagation with velocity during time dt
****************************************************************************/
#pragma acc routine seq
double FC_xzrot_dt(double x, double z, double vx, double vz,
double t, double dt, char type, struct FermiChopper_struct FCs)
{
if (dt) /* with propagation */
return( (type == 'x' ? FC_xrot(x+vx*dt, z+vz*dt, t+dt, FCs)
: FC_zrot(x+vx*dt, z+vz*dt, t+dt, FCs)) );
else /* without propagation */
return( (type == 'x' ? FC_xrot(x,z,t,FCs)
: FC_zrot(x,z,t,FCs)) );
}
/*****************************************************************************
* FC_xzbrent(x,z,vx,vz, t,dt, type='x' or 'z', d, FCs)
* solves X'=d and Z'=d with Brent algorithm in time interval [0, dt].
* Returns time within [0,dt], from NumRecip in C, chap 9, p360 (zbrent)
* ERRORS: return -1 not used
* -2 if exceed MAX iteration
* -3 no sign change in range
****************************************************************************/
#pragma acc routine seq
double FC_xzbrent(double x, double z, double vx, double vz,
double t, double dt,
char type, double d, struct FermiChopper_struct FCs)
{
int iter;
double a=0,b=dt;
double c,dd,e,min1,min2;
double tol=FermiChopper_TimeAccuracy;
double EPS=FermiChopper_TimeAccuracy;
double fa=FC_xzrot_dt(x,z,vx,vz, t,a, type, FCs) - d;
double fb=FC_xzrot_dt(x,z,vx,vz, t,b, type, FCs) - d;
double fc,p,q,r,s,tol1,xm;
if (fb*fa > 0.0) return -3;
fc=fb;
for (iter=1;iter<=FermiChopper_MAXITER;iter++) {
if (fb*fc > 0.0) {
c=a;
fc=fa;
e=dd=b-a;
}
if (fabs(fc) < fabs(fb)) {
a=b;
b=c;
c=a;
fa=fb;
fb=fc;
fc=fa;
}
tol1=2.0*EPS*fabs(b)+0.5*tol;
xm=0.5*(c-b);
if (fabs(xm) <= tol1 || fb == 0.0) return b;
if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
s=fb/fa;
if (a == c) {
p=2.0*xm*s;
q=1.0-s;
} else {
q=fa/fc;
r=fb/fc;
p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
q=(q-1.0)*(r-1.0)*(s-1.0);
}
if (p > 0.0) q = -q;
p=fabs(p);
min1=3.0*xm*q-fabs(tol1*q);
min2=fabs(e*q);
if (2.0*p < (min1 < min2 ? min1 : min2)) {
e=dd;
dd=p/q;
} else {
dd=xm;
e=dd;
}
} else {
dd=xm;
e=dd;
}
a=b;
fa=fb;
if (fabs(dd) > tol1)
b += dd;
else
b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
fb=FC_xzrot_dt(x,z,vx,vz, t,b, type, FCs) - d;
}
return -2;
} /* FC_xzbrent */
/*****************************************************************************
* Wrappers to intersection algorithms
****************************************************************************/
#pragma acc routine seq
double FC_xintersect(double x, double z, double vx, double vz,
double t, double dt,
double d, struct FermiChopper_struct FCs)
{
return(FC_xzbrent(x, z, vx, vz, t, dt, 'x', d, FCs));
}
#pragma acc routine seq
double FC_zintersect(double x, double z, double vx, double vz,
double t, double dt,
double d, struct FermiChopper_struct FCs)
{
return(FC_xzbrent(x, z, vx, vz, t, dt, 'z', d, FCs));
}
#endif
%}
DECLARE
%{
struct FermiChopper_struct FCVars;
%}
INITIALIZE
%{
/************************ CALCULATION CONSTANTS *****************************/
strcpy(FCVars.compcurname, NAME_CURRENT_COMP);
FCVars.omega = 2*PI*nu;
if (!phase && delay) {
FCVars.ph0= fmod(-delay*nu*360,360)*DEG2RAD;
} else FCVars.ph0 = phase*DEG2RAD;
FCVars.sum_t=FCVars.sum_v=FCVars.sum_N=FCVars.sum_N_pass=0;
/* check of input parameters */
if (nslit < 1) nslit=1;
if (yheight <= 0) exit(printf("FermiChopper: %s: FATAL: unrealistic cylinder yheight =%g [m]\n", NAME_CURRENT_COMP, yheight));
if (m <= 0) { m=0; R0=0; }
if (radius <= 0) {
printf("FermiChopper: %s: FATAL: Unrealistic cylinder radius radius=%g [m]\n", NAME_CURRENT_COMP, radius);
exit(-1);
}
if (xwidth > 0 && xwidth < radius*2 && nslit > 0) {
w = xwidth/nslit;
}
if (w <= 0) {
printf("FermiChopper: %s: FATAL: Slits in the package have unrealistic width w=%g [m]\n", NAME_CURRENT_COMP, w);
exit(-1);
}
if (nslit*w > radius*2) {
nslit = floor(radius/w);
printf("FermiChopper: %s: Too many slits to fit in the cylinder\n"
" Adjusting nslit=%f\n", NAME_CURRENT_COMP, nslit);
}
if (length > radius*2) {
length = 2*sqrt(radius*radius - nslit*w*nslit*w/4);
printf("FermiChopper: %s: Slit package is longer than the whole\n"
" chopper cylinder. Adjusting length=%g [m]\n", NAME_CURRENT_COMP, length);
}
if (eff <= 0 || eff > 1) {
eff = 0.95;
printf("FermiChopper: %s: Efficiency is unrealistic\n"
" Adjusting eff=%f\n", NAME_CURRENT_COMP, eff);
}
if (Qc <= 0) { Qc = 0.02176; m = 0; R0 = 0; }
if (W <= 0) W=1e-6;
if (curvature) {
FCVars.C_slit = curvature;
if (1 < fabs(radius*curvature))
exit(printf("FermiChopper: %s: Slit curvature is unrealistic\n",
NAME_CURRENT_COMP));
}
FCVars.L_slit = length;
if (verbose && nu)
printf("FermiChopper: %s: Frequency nu=%g [Hz] %g [rpm], time frame=%g [s] phase=%g [deg]\n"
, NAME_CURRENT_COMP, nu, nu*60, 2/nu, FCVars.ph0*RAD2DEG);
FCVars.absorb_alreadyinside = 0;
FCVars.absorb_topbottom = 0;
FCVars.absorb_cylentrance = 0;
FCVars.absorb_sideentrance = 0;
FCVars.absorb_notreachentrance = 0;
FCVars.absorb_packentrance = 0;
FCVars.absorb_slitcoating = 0;
FCVars.warn_notreachslitwall = 0;
FCVars.absorb_exitslitpack = 0;
FCVars.absorb_maxiterations = 0;
FCVars.absorb_wrongdirection = 0;
FCVars.absorb_nocontrol = 0;
FCVars.absorb_cylexit = 0;
FCVars.warn_notreachslitoutput = 0;
/* fix for the wrong coordinate frame orientation to come back to McStas XYZ system */
FCVars.omega *= -1;
FCVars.ph0 *= -1;
FCVars.t0 = -FCVars.ph0/FCVars.omega;
%}
TRACE
%{
/** local CALCULATION VARIABLES**************************************/
/** Interaction with slit package ***************************/
double slit_input; /* length of the slits */
/** Variables for calculating interaction with blades ***************/
double xp1, zp1, vxp1;
/** Reflections ***********************************************/
double n1;
/** Multiple Reflections ******************************/
int loopcounter=0; /* How many reflections happen? */
/** Time variables *********************************/
double t3=0; /* interaction time at n3 position (slit wall) */
double t1=0,t2=0; /* cylinder intersection time (at entry and exit of slit pack - n1 n2 - or cylinder) */
double dt; /* interaction intervals (e.g. time for exiting slit pack) */
double X[5],Z[5]; /* position of interaction locations in the Fermi chopper rotating frame
[0]=cylinder input
[1]=slit input
[2]=slit wall reflection
[3]=slit exit
[4]=cylinder exit
*/
double ref = 0;
double par[5] = {R0, Qc, alpha, m, W};
/************************ TIME OF FLIGHT RESET ************************/
if (zero_time == 1 && nu)
t -= floor( (t+1/(4*nu))*(2*nu) )/(2*nu);
/* zero arrays used to store positions */
for (loopcounter=0; loopcounter < 5; loopcounter++)
X[loopcounter] = Z[loopcounter]=0;
/************** test, if the neutron interacts with the cylinder ***/
if (cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight))
{
if (t1 <= 0) { /* Neutron started inside the cylinder */
if (verbose > 0 && FCVars.absorb_alreadyinside 2)
printf("FermiChopper: %s: t1=%8.3g t2=%8.3g xyz=[%8.3g %8.3g %8.3g] v=[%8.3g %8.3g %8.3g] t=%8.3g (init).\n",
NAME_CURRENT_COMP, t1, t2, x,y,z,vx,vy,vz,t);
dt=t2-t1; /* total time of flight inside the cylinder */
PROP_DT(t1); /* Propagates neutron to entrance of the cylinder */
SCATTER;
if (verbose > 2)
printf("FermiChopper: %s: PROP_DT t1=%8.3g t2=%8.3g xyz=[%8.3g %8.3g %8.3g] v=[%8.3g %8.3g %8.3g] t=%8.3g (IN cyl).\n",
NAME_CURRENT_COMP, t1, t2, x,y,z,vx,vy,vz,t);
/* neutron must not enter or leave from top or bottom of cylinder. */
if (fabs(y) >= yheight/2.0 || fabs(y+vy*dt) >= yheight/2.0) {
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (enter).\n",
NAME_CURRENT_COMP, y);
#pragma acc atomic
FCVars.absorb_topbottom = FCVars.absorb_topbottom + 1;
ABSORB;
}
vxp1 = sqrt(vx*vx+vy*vy+vz*vz);
double tmp=p*vxp1;
#pragma acc atomic
FCVars.sum_v = FCVars.sum_v + tmp;
tmp = p*t;
#pragma acc atomic
FCVars.sum_t = FCVars.sum_t + tmp;
#pragma acc atomic
FCVars.sum_N = FCVars.sum_N + p;
if (zero_time > 1 && FCVars.sum_N) { /* automatic phase adjustment */
double mean_t, mean_phase;
mean_t = FCVars.sum_t/FCVars.sum_N;
mean_phase = fmod(mean_t*nu*2*PI, 2*PI);
/* now we shift the phase so that the neutron pulse is centered on the slit pack */
mean_phase+= radius/vxp1*2*PI*nu;
#pragma acc atomic write
FCVars.ph0 = mean_phase;
}
/* neutron must enter the cylinder opening: |X'| < full package width*/
xp1 = FC_xrot(x,z, t,FCVars); /* X'(t) */
if (fabs(xp1) >= nslit*w/2) {
if (verbose > 2)
printf("FermiChopper: %s: ABSORB Neutron X is outside cylinder aperture, x'=%8.3g > %g (enter).\n",
NAME_CURRENT_COMP, xp1, nslit*w/2);
#pragma acc atomic
FCVars.absorb_cylentrance = FCVars.absorb_cylentrance + 1;
ABSORB;
}
/*********************** PROPAGATE TO SLIT PACKAGE **************************/
/* zp1 = Z' at entrance of cylinder Z'(t) */
zp1 = FC_zrot(x,z, t, FCVars);
X[0] = xp1; Z[0] = zp1;
/* here we should have sqrt(x*x+z*z) == sqrt(xp1*xp1+zp1*zp1) */
t3 = sqrt(x*x+z*z) / sqrt(xp1*xp1+zp1*zp1);
if (t3 < 0.99 || 1.01 < t3) {
if (verbose > 1 && FCVars.absorb_cylentrance < FermiChopper_MAXITER)
printf("FermiChopper: %s: ABSORB Neutron radius on cylinder in static frame r=%g does not match that of rotating frame r'=%g.\n",
NAME_CURRENT_COMP, sqrt(x*x+z*z), sqrt(xp1*xp1+zp1*zp1));
#pragma acc atomic
FCVars.absorb_cylentrance = FCVars.absorb_cylentrance +1 ;
ABSORB;
}
/* Checking on which side of the Chopper the Neutron enters: sign(Z') */
slit_input = (zp1 > 0 ? length/2 : -length/2);
/* time shift to reach slit package in [0,time to exit cylinder]: Z'=slit_input */
/* t3 is used here as a tmp variable, will be redefined in for loop */
t3 = FC_zintersect(x,z,vx,vz, t,dt, slit_input, FCVars);
if( (t3 < 0)||(t3 > dt) ) {
if (verbose > 2 && FCVars.absorb_notreachentrance < FermiChopper_MAXITER) {
printf("FermiChopper: %s: Can not reach entrance of slits. dt=%8.3g t3=%8.3g (intersection:1).\n",
NAME_CURRENT_COMP, dt, t3);
if (t3 == -3)
printf(" No sign change to determine intersection\n");
else if (t3 == -2)
printf(" Max iterations reached\n");
else if (t3 < 0)
printf(" Error when solving intersection\n");
}
#pragma acc atomic
FCVars.absorb_notreachentrance = FCVars.absorb_notreachentrance + 1;
ABSORB; /* neutron can not reach slit entrance */
}
/* Propagating to the slit package entrance */
PROP_DT(t3); /* dt = t2-t1: time in cylinder */
dt -= t3; /* remaining time from slit pack entry to exit of cylinder */
xp1 = FC_xrot(x,z, t, FCVars); /* should be slit_input */
zp1 = FC_zrot(x,z, t, FCVars);
X[1] = xp1; Z[1] = zp1;
#ifndef OPENACC
if (mcdotrace) {
/* indicate position of neutron in mcdisplay */
double xp2 = x; double zp2 = z; x = xp1; z=zp1; SCATTER; x=xp2; z=zp2;
} else
#endif
SCATTER;
if (verbose > 2)
printf("FermiChopper: %s: PROP_DT t=%8.3g dt=%8.3g x'=%8.3g z'=%8.3g length=%g (slit enter).\n",
NAME_CURRENT_COMP, t, dt, xp1, zp1, slit_input);
/* must have X'< slit package width at package Z'=slit_input */
if (fabs(xp1) >= nslit*w/2) {
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package, x'=%8.3g > %g (enter).\n",
NAME_CURRENT_COMP, xp1, nslit*w/2);
#pragma acc atomic
FCVars.absorb_packentrance = FCVars.absorb_packentrance + 1;
ABSORB;
}
/* solve Z'=-slit_input for time of exit of slit package */
/* t3 is used here as a tmp variable, will be redefined in for loop */
t3 = FC_zintersect(x,z,vx,vz, t,dt*1.1, -slit_input, FCVars);
if((t3 < FermiChopper_TimeAccuracy)||(t3 > dt)) {
if (verbose > 1 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) {
printf("FermiChopper: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g (intersection:2).\n",
NAME_CURRENT_COMP, dt, t3);
if (t3 == -3)
printf(" No sign change to determine intersection\n");
else if (t3 == -2)
printf(" Max iterations reached\n");
else if (t3 < 0)
printf(" Error when solving intersection\n");
}
#pragma acc atomic
FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1;
/* we estimate analytically the time to the slit exit */
t3 = length/vxp1;
}
dt = t3; /* reduce time interval to [0, time of slit exit] */
/* here we should have dt*v = length (slit entrance -> exit) */
t3 = fabs(FC_xzrot_dt(x,z,vx,vz, t, 0 , 'z', FCVars) - FC_xzrot_dt(x,z,vx,vz, t, dt, 'z', FCVars))/length;
if (fabs(t3-1) > 0.02) {
if (verbose > 0 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER)
printf("FermiChopper: %s: ABSORB Neutron propagation time v*dt/length=%g in slit does not match its length=%g (slit exit expected).\n",
NAME_CURRENT_COMP, t3, length);
#pragma acc atomic
FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1;
ABSORB;
}
/* dt= time shift to go to exit of slit package (or exit of cylinder in case of error) */
/*
|---------|
| / |
|o * (dt)
| |
|---------|
*/
/*********************PROPAGATION INSIDE THE SLIT PACKAGE *******************/
/* index of slit hit at entrance n1=-N/2 (xp1=-) ... N/2 (xp1=+) */
n1 = floor(xp1/w);
/******************* BEGIN LOOP FOR MULTIPLE REFLECTIONS ********************/
for (loopcounter=0; loopcounter<=FermiChopper_MAXITER; loopcounter++) {
double dt_to_tangent=0; /* time shift to go to tangent intersection */
double n2,n3; /* slit indices */
double xp2, zp2,xp3=0,vxp2=0,vzp1=0; /* position, velocity */
double q; /* used for calculating new velocities after reflection */
int i;
/* compute trajectory tangents: m1=Vz'+w.X'(t), m2=Vz'+w.X'(t+dt) */
xp1 = FC_xrot (x,z, t, FCVars); /* X'(t) current position */
xp2 = FC_xzrot_dt(x,z,vx,vz, t, dt, 'x', FCVars); /* X'(t+dt) slit exit */
zp2 = FC_xzrot_dt(x,z,vx,vz, t, dt, 'z', FCVars); /* Z'(t+dt) slit exit */
/* slit index at the end of the slit: */
n2 = floor(xp2/w);
/* quick exit for absorbing walls when neutron changes slit index */
if (n2 != n1 && (m <= 0 || R0 <= 0 || Qc <= 0)) {
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (change slit).\n",
NAME_CURRENT_COMP);
#pragma acc atomic
FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
ABSORB;
}
/* compute transversal velocity to determine their intersection */
vxp1= FC_xrot (vx+z*FCVars.omega,vz-x*FCVars.omega,
t, FCVars); /* dX'(t)/dt slope at current position*/
vxp2= FC_xrot (vx+(z+vz*dt)*FCVars.omega,vz-(x+vx*dt)*FCVars.omega,
t+dt,FCVars); /* dX'(t+dt)/dt slope at slit exit */
/* absolute time at tangent intersection, changed to time shift below */
dt_to_tangent = (vxp1 - vxp2 ? (xp2 - xp1 - dt*vxp2)/(vxp1 - vxp2) : -1);
/* If algorithm fails, take the middle of the interval*/
if (dt_to_tangent < 0 || dt_to_tangent > dt) {
if (verbose > 1)
printf("FermiChopper: %s: WARNING could not determine tangent intersection dt=%g. Using middle.\n",
NAME_CURRENT_COMP, dt_to_tangent);
dt_to_tangent=dt*0.5;
}
/*
*(dt_to_tangent, xp3)
|---------|
| / \ |
xp1 |o \|(dt) xp2
| |
|---------|
*/
/* point coordinates at tangent intersection/middle point (max deviation from optical axis) */
xp3 = FC_xzrot_dt(x,z,vx,vz, t, dt_to_tangent, 'x', FCVars); /* X'(t+dt_to_tangent) */
/* slit index at the tangent intersection/middle point */
n3 = floor(xp3/w);
if (verbose > 2)
printf("FermiChopper: %s: t3=%8.3g slit_indices=[%g %g %g] (time at tangent intersection).\n",
NAME_CURRENT_COMP, dt_to_tangent, n1, n2, n3);
/* change slit means there is a reflection/intersection inside */
if ( n2!=n1 || n3!=n1 ) {
double t3a, t3b, distance_Wa, distance_Wb;
if (m <= 0 || R0 <= 0 || Qc <= 0) {
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (change slit).\n",
NAME_CURRENT_COMP);
#pragma acc atomic
FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
ABSORB;
}
/* Choosing the first time it isn't in the slit anymore */
if(n2 == n1 && n3 != n1){
n2 = n3;
}
/* get position of slit wall towards which neutron is propagating */
if (n2 > n1) { /* X' positive side of slit is in principle the first intersection to test*/
distance_Wa = n1*w+w;
distance_Wb = n1*w;
} else { /* X' negative side of slit */
distance_Wb = n1*w+w;
distance_Wa = n1*w;
}
/* time shift to reach slit wall point in [0,dt_to_tangent]: X'=distance_W slit_wall */
for (i=0; i< 2; i++) {
/* first attempt: [0,dt_to_tangent] (to max deviation location)
second attempt: [0, dt] (to slit exit) */
double dt_search = (i == 0 ? dt_to_tangent : dt);
t3a = FC_xintersect(x,z,vx,vz, t,dt_to_tangent, distance_Wa, FCVars);
t3b = FC_xintersect(x,z,vx,vz, t,dt_to_tangent, distance_Wb, FCVars);
if (t3b < 0) t3 = t3a;
else if (t3a < 0 && t3b >= 0) t3 = t3b;
else t3 = (t3a < t3b ? t3a : t3b);
if (FermiChopper_TimeAccuracy < t3 && t3 < dt_search)
break; /* we found the intersection location */
}
/* handle case where intersection search fails */
if (t3 < FermiChopper_TimeAccuracy || t3 >= dt_to_tangent) {
if (verbose > 1 && FCVars.warn_notreachslitwall < FermiChopper_MAXITER) {
printf("FermiChopper: %s: Can not reach slit wall (iteration %i). dt=%8.3g t3=%8.3g (intersection:3).\n",
NAME_CURRENT_COMP, loopcounter, dt_to_tangent, t3);
if (t3 == -3)
printf(" No sign change to determine intersection\n");
else if (t3 == -2)
printf(" Max iterations reached\n");
else if (t3 < 0 || t3 >= dt_to_tangent)
printf(" Error when solving intersection\n");
}
/* neutron can not reach slit wall. */
#pragma acc atomic
FCVars.warn_notreachslitwall = FCVars.warn_notreachslitwall + 1;
ABSORB;
}
/* Propagate to slit wall point (t+t3) on slit n3 wall */
PROP_DT(t3); dt -= t3; /* dt: time remaining to slit exit after propagation */
xp1 = FC_xrot(x,z, t, FCVars); /* X'(t+t3) : on slit wall */
zp1 = FC_zrot(x,z, t, FCVars); /* Z'(t+t3) : on slit wall */
X[2] = xp1; Z[2] = zp1;
if (verbose > 2)
printf("FermiChopper: %s: PROP_DT t3=%8.3g dt=%8.3g xyz=[%8.3g %8.3g %8.3g] (on wall). z'=%g\n",
NAME_CURRENT_COMP, t3, dt, x,y,z, zp1);
/* check if intersection point is still in the slit package, else exit loop */
if (fabs(zp1) > length/2 || dt <= 0) {
if (verbose > 2)
printf("FermiChopper: %s: Neutron is outside slit pack (on slit wall).\n",
NAME_CURRENT_COMP);
break;
}
/* here
|---o-----|
| / \ |
|/ \ |
| \ |(dt)
|---------|
*/
/* get velocity in rotating frame, on slit wall */
vxp1 = FC_xrot(vx,vz, t, FCVars);
vzp1 = FC_zrot(vx,vz, t, FCVars);
q = 2*V2Q*(fabs(vxp1));
{
/* double ref = 0;
double par[] = {R0, Qc, alpha, m, W};*/
StdReflecFunc(q, par, &ref);
if (ref>0) p *= ref;
else {
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits absorbing coating (on slit wall).\n",
NAME_CURRENT_COMP);
#pragma acc atomic
FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
ABSORB;
} /* Cutoff ~ 1E-10 */
}
#ifndef OPENACC
if (mcdotrace) {
double xp2 = x; double zp2 = z;
/* indicate position of neutron in mcdisplay */
x = FC_xrot(x,z, t,FCVars); z= FC_zrot(x,z, t,FCVars); SCATTER; x=xp2; z=zp2;
} else
#endif
SCATTER;
/* reflect perpendicular velocity and compute new velocity in static frame */
vxp1 *= -1;
/* apply transposed Transformation matrix */
vx = FC_xrot( vxp1,-vzp1, t,FCVars);
vz = FC_zrot(-vxp1, vzp1, t,FCVars);
/* recompute time to slit exit */
/* solve Z'=-slit_input for time of exit of slit package */
t3 = FC_zintersect(x,z,vx,vz, t,dt, -slit_input, FCVars);
if(t3 < 0 || t3 > dt) {
if (verbose > 1 && FCVars.warn_notreachslitoutput < FermiChopper_MAXITER) {
printf("FermiChopper: %s: Can not reach exit of slits. dt=%8.3g t3=%8.3g (intersection:4).\n",
NAME_CURRENT_COMP, dt, t3);
if (t3 == -3)
printf(" No sign change to determine intersection\n");
else if (t3 == -2)
printf(" Max iterations reached\n");
else if (t3 < 0)
printf(" Error when solving intersection\n");
}
#pragma acc atomic
FCVars.warn_notreachslitoutput = FCVars.warn_notreachslitoutput + 1;
ABSORB;
/* neutron can not reach slit output. */
} else dt = t3; /* reduce time interval to [0, time of slit exit] */
} /* end if changed slit index */
else { /* neutron remains in the same slit: direct flight */
if (dt > 0) PROP_DT(dt); /* go to slit exit */
break;
}
} /* end for */
xp1 = FC_xrot(x,z, t,FCVars);
zp1 = FC_zrot(x,z, t,FCVars);
X[3] = xp1; Z[3] = zp1; /* slit exit */
if (fabs(xp1) >= nslit*w/2)
{
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package, x'=%8.3g (slit exit).\n",
NAME_CURRENT_COMP, xp1);
#pragma acc atomic
FCVars.absorb_slitcoating = FCVars.absorb_slitcoating + 1;
ABSORB;
}
if (loopcounter >= FermiChopper_MAXITER)
{
if (verbose > 1 && FCVars.absorb_maxiterations < FermiChopper_MAXITER)
printf("FermiChopper: %s: Max iterations %i reached inside slit. Absorb.\n",
NAME_CURRENT_COMP, FermiChopper_MAXITER);
#pragma acc atomic
FCVars.absorb_maxiterations = FCVars.absorb_maxiterations + 1;
ABSORB;
}
/********************* EXIT SLIT PACKAGE ********************************/
/* propagation times to cylinder. Should use t2 to exit */
if (!cylinder_intersect (&t1, &t2, x, y, z, vx, vy, vz, radius, yheight)) {
/* must be inside the cylinder */
if (verbose > 1) printf("FermiChopper: %s: ABSORB Neutron has unexpectidely exited cylinder ! (exiting)\n",
NAME_CURRENT_COMP);
#pragma acc atomic
FCVars.absorb_exitslitpack = FCVars.absorb_exitslitpack + 1;
ABSORB;
}
if (t1 > 0)
{
if (verbose > 1 && FCVars.absorb_wrongdirection < FermiChopper_MAXITER)
printf("FermiChopper: %s: Neutrons are leaving chopper\n"
" in the wrong direction. Absorb.\n", NAME_CURRENT_COMP);
#pragma acc atomic
FCVars.absorb_wrongdirection = FCVars.absorb_wrongdirection + 1;
ABSORB;
}
if (t2 <= 0 && FCVars.absorb_nocontrol < FermiChopper_MAXITER)
{
if (verbose > 1)
printf("FermiChopper: %s: Neutrons are leaving chopper\n"
" without any control. Absorb.\n", NAME_CURRENT_COMP);
#pragma acc atomic
FCVars.absorb_nocontrol = FCVars.absorb_nocontrol + 1;
ABSORB;
}
/* propagate to cylinder surface (exit) */
PROP_DT(t2);
SCATTER;
xp1 = FC_xrot(x,z, t,FCVars);
zp1 = FC_zrot(x,z, t,FCVars);
X[4] = xp1; Z[4] = zp1;
if (verbose > 2)
printf("FermiChopper: %s: t1=%8.3g PROP_DT t2=%8.3g xyz=[%8.3g %8.3g %8.3g] (OUT cyl). z'=%g\n",
NAME_CURRENT_COMP, t1, t2, x,y,z, zp1);
/* Check if the neutron left the cylinder by its top or bottom */
if (fabs(y) >= yheight/2)
{
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron hits top/bottom of cylinder, y=%8.3g (exiting)\n",
NAME_CURRENT_COMP, y);
#pragma acc atomic
FCVars.absorb_topbottom = FCVars.absorb_topbottom + 1;
ABSORB;
}
/* must have X'< slit package width at package Z'=cylinder output */
if (fabs(xp1) >= nslit*w/2)
{
if (verbose > 2) printf("FermiChopper: %s: ABSORB Neutron X is outside slit package cylinder, xp1=%8.3g (exiting).\n",
NAME_CURRENT_COMP, xp1);
#pragma acc atomic
FCVars.absorb_cylexit = FCVars.absorb_cylexit + 1;
ABSORB;
}
/* Transmission coefficent */
p = p*eff; //finite cross section + transmission
#pragma acc atomic
FCVars.sum_N_pass = FCVars.sum_N_pass + p;
} /* end if cylinder_intersect */
%}
SAVE
%{
double mean_k, mean_v, mean_t, mean_w=0, mean_L=0;
if (FCVars.sum_N) {
mean_v = FCVars.sum_v/FCVars.sum_N;
mean_t = FCVars.sum_t/FCVars.sum_N;
mean_k = V2K*mean_v;
if (mean_k) mean_L = 2*PI/mean_k;
mean_w = VS2E*mean_v*mean_v;
/* base opening time */
double div, mean_phase;
div = atan2(w,length)/PI*180;
mean_phase = fmod(mean_t*nu*360, 360);
mean_phase+=radius/mean_v*360*nu;
if (mean_phase > 180) mean_phase -= 360;
if (!FCVars.sum_N_pass)
printf("FermiChopper: %s: No neutron can pass the chopper.\n", NAME_CURRENT_COMP);
if (!FCVars.sum_N_pass || verbose)
printf("FermiChopper: %s\n"
" Mean velocity v = %g [m/s]\n"
" Mean wavelength lambda= %g [Angs]\n"
" Mean energy omega = %g [meV]\n"
" Mean arrival time t = %g [s]\n"
" Mean phase = %g [deg] (%s)\n"
" Slit pack divergence = %g [deg] (full width)\n"
" Opening time dt = %g [s]\n"
" Intensity reaching FC = %g [n/s]\n"
" Intensity passing FC = %g [n/s]\n"
, NAME_CURRENT_COMP,
mean_v, mean_L, mean_w, mean_t, mean_phase,
(zero_time > 1 ? "set automatically" : "use phase=-(this value) to optimize"),
2*div,
(nu ? fabs(div/PI/nu) : 1),
FCVars.sum_N,
FCVars.sum_N_pass);
if (!FCVars.sum_N_pass || verbose) {
printf("FermiChopper: %s: Lost events anaylsis\n"
" Already inside: %li\n"
" By Top/Bottom of cylinder: %li\n"
" At cylinder entrance: %li\n"
" Hit cyl. entrance sides: %li\n"
" Can't prop. to slit pack: %li\n"
" At slit pack entrance: %li\n"
" On absorbing slit coating: %li\n"
" Exiting slit pack: %li\n"
" Too many iterations: %li\n"
" Prop. in wrong direction : %li\n"
" Mad neutron (no control): %li\n"
" At cylinder exit: %li\n"
, NAME_CURRENT_COMP,
FCVars.absorb_alreadyinside,
FCVars.absorb_topbottom,
FCVars.absorb_cylentrance,
FCVars.absorb_sideentrance,
FCVars.absorb_notreachentrance,
FCVars.absorb_packentrance,
FCVars.absorb_slitcoating,
FCVars.absorb_exitslitpack,
FCVars.absorb_maxiterations,
FCVars.absorb_wrongdirection,
FCVars.absorb_nocontrol,
FCVars.absorb_cylexit);
if (FCVars.warn_notreachslitwall || FCVars.warn_notreachslitoutput)
printf("Warning: Can not reach slit wall: %li\n"
"Warning: Can not reach slit output: %li\n",
FCVars.warn_notreachslitwall,
FCVars.warn_notreachslitoutput);
}
} else {
printf("FermiChopper: %s: No neutron can reach the chopper.\n", NAME_CURRENT_COMP);
}
%}
MCDISPLAY
%{
double index_x=0;
double index_z=0;
double xpos, zpos;
double Nz,Nx;
double ymin=-yheight/2.0;
double ymax= yheight/2.0;
double omega=FCVars.omega;
FCVars.omega=0;
// FCVars.ph0 =0;
Nz = (FCVars.C_slit ? 4 : 1);
Nx = (nslit > 11 ? 11 : nslit);
FCVars.C_slit *= -1;
/* cylinder top/center/bottom */
circle("xz", 0,ymax,0,radius);
circle("xz", 0,0 ,0,radius);
circle("xz", 0,ymin,0,radius);
/* vertical lines to make a kind of volume */
line( 0 ,ymin,-radius, 0 ,ymax,-radius);
line( 0 ,ymin, radius, 0 ,ymax, radius);
line(-radius,ymin, 0 ,-radius,ymax, 0 );
line( radius,ymin, 0 , radius,ymax, 0 );
/* slit package */
for (index_x = -Nx/2; index_x < Nx/2; index_x++) {
for (index_z = -Nz/2; index_z < Nz/2; index_z++) {
double xs1, zs1, zs2;
double xp1, xp2, zp1, zp2;
zs1 = length*index_z/Nz;
zs2 = length*(index_z+1)/Nz;
xs1 = w*nslit*index_x/Nx;
xp1 = FC_xrot(xs1, zs1, 0, FCVars);
xp2 = FC_xrot(xs1, zs2, 0, FCVars);
zp1 = FC_zrot(xs1, zs1, 0, FCVars);
zp2 = FC_zrot(xs1, zs2, 0, FCVars);
multiline(5, xp1, ymin, zp1,
xp1, ymax, zp1,
xp2, ymax, zp2,
xp2, ymin, zp2,
xp1, ymin, zp1);
}
}
/* cylinder inner sides containing slit package */
double xp1, xp2, zp1, zp2;
xpos = nslit*w/2;
zpos = sqrt(radius*radius - xpos*xpos);
xp1 = FC_xrot(xpos, -zpos, 0, FCVars);
xp2 = FC_xrot(xpos, +zpos, 0, FCVars);
zp1 = FC_zrot(xpos, -zpos, 0, FCVars);
zp2 = FC_zrot(xpos, +zpos, 0, FCVars);
multiline(5, xp1, ymin, zp1,
xp1, ymax, zp1,
xp2, ymax, zp2,
xp2, ymin, zp2,
xp1, ymin, zp1);
xpos *= -1;
xp1 = FC_xrot(xpos, -zpos, 0, FCVars);
xp2 = FC_xrot(xpos, +zpos, 0, FCVars);
zp1 = FC_zrot(xpos, -zpos, 0, FCVars);
zp2 = FC_zrot(xpos, +zpos, 0, FCVars);
multiline(5, xp1, ymin, zp1,
xp1, ymax, zp1,
xp2, ymax, zp2,
xp2, ymin, zp2,
xp1, ymin, zp1);
FCVars.omega=omega;
%}
END