1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
|
/******************************************************************
*
* McStas, version 3.0
*
* Component: vertical_T0
*
* %Identification
* Written by: Garrett Granroth
* Date: 2 NOV 2004
* Origin: SNS Oak Ridge,TN
* Version: 0.4
* %Parameters
* Input Parameters:
*
* len: length of slot (m)
* w1: center width (m)
* w2: edgewidth
* nu: frequency (Hz)
* delta: time from edge of chopper to center Phase angle (sec)
* tc: time when desired neutron is at the center of the chopper (sec)
* ymin: Lower y bound (m)
* ymax: Upper y bound (m)
*
*
* %End
*******************************************************************/
DEFINE COMPONENT Vertical_T0a
SETTING PARAMETERS (len, w1,w2, nu, delta, tc, ymin, ymax)
SHARE
%{
#ifndef FERMI_CHOP_DEFS
#define FERMI_CHOP_DEFS
/* routine to calculate acos in proper quadrant range = 0 to 2PI*/
#pragma acc routine
double acos0_2pi(double x,double y)
{
if (y>0.0){
return acos(x);
}
return 2.0*PI-acos(x);
}
/*routine to calculate x and y positions of a neutron in a fermi chopper */
#pragma acc routine
void neutxypos(double *x, double *y, double phi, double inrad, double* c)
{
*x=c[0]+inrad*cos(phi);
*y=c[1]+inrad*sin(phi);
}
/* routine to calculate the origin of a circle that describes the neutron path through the chopper */
#pragma acc routine
void calccenter(double* c, double* zr, double* xr){
double denom, A,B,C,D,a,b;
denom=2*(-zr[0]*xr[2] +zr[0]*xr[1]+ zr[1]*xr[2]+xr[0]*zr[2]-xr[0]*zr[1] - xr[1]*zr[2]);
A=xr[1]-xr[2];B=xr[0]-xr[1];C=zr[2]-zr[1];D=zr[1]-zr[0];
a=zr[0]*zr[0]-zr[1]*zr[1]+xr[0]*xr[0]-xr[1]*xr[1];
b=zr[2]*zr[2]-zr[1]*zr[1]+xr[2]*xr[2]-xr[1]*xr[1];
c[0]=1.0/denom*(A*a+B*b);
c[1]=1.0/denom*(C*a+D*b);
}
#endif
/* function to calculate if the neutron is in the channel or not
* return 0 if neutron does not transmit return 1 if neutron will pass*/
int t0checkabsorb(double phi, double inrad,double inw1, double inw2, double* c){
double xtmp,neuzr,neuxr;
neutxypos(&neuzr,&neuxr,phi,inrad,c);
// printf("xr:%g zr:%g phi: %g r: %g c[0]: %g c[1]: %g\n",neuxr,neuzr,phi,inrad,c[0],c[1]);
if (fabs(neuxr)>inw1/2.0+(inw2-inw1)/(inrad/2.0)*fabs(neuzr)) // check if neutron x position is outside of channel
return 0;
return 1;
}
%}
DECLARE
%{
double omega;
double off;
double splen;
double rad;
double sw;
%}
INITIALIZE
%{
splen=len/2.0;
omega=2.0*PI*nu;
rad=sqrt(w2*w2/4.0+splen*splen); //radius of cylinder containing slit package.
%}
TRACE
%{
double t0,t1,dphi,dt2,tneuzr,tneuxr,nrad;
double phivec[200],tpt[3],xpt[3],ypt[3],zpt[3],zr[3],xr[3],yr[3],theta[3],c[2];
int chan_num,chan_num0,idx1,idx3;
if (cylinder_intersect (&t0, &t1, x, y, z, vx, vy, vz, rad, ymax-ymin)){
if (t0 < 0) /*Neutron started inside cylinder */
ABSORB;
dt2=t1-t0;
PROP_DT(t0); /*propagate neutron to edge of chopper*/
/*calculate neutron position and velocity in chopper frame
calculate 3 points in the instrument frame and put them into the
chopper frame inorder to determine the radius and center of a circle
that describes the path of the neutron in the chopper frame. */
tpt[1]=t;
tpt[2]=t+dt2;
tpt[0]=t+dt2/2.0;
//set local 0 in time as tc and calculate angle of rotation for each point
for(idx3=0;idx3<3;idx3++){
theta[idx3]=(tpt[idx3]-tc)*omega;
}
zpt[1]=-sqrt(rad*rad-x*x); xpt[1]=x; ypt[1]=y; /* point where neutron intersects chopper */
zpt[2]=zpt[1]+vz*(dt2); xpt[2]=xpt[1]+vx*(dt2); ypt[2]=ypt[1]+vy*(dt2); /* point where neutron leaves the chopper */
xpt[0]=xpt[1]+vx*(dt2/2.0); ypt[0]=ypt[1]+vy*(dt2/2.0); zpt[0]=zpt[1]+vz*(dt2/2.0); /*point half way between in time */
/* do the rotation */
for(idx3=0;idx3<3;idx3++){
rotate(xr[idx3],yr[idx3],zr[idx3],xpt[idx3],ypt[idx3],zpt[idx3],theta[idx3],0,1,0);
}
calccenter(c,zr,xr); /* calculate the center */
nrad=sqrt((zr[0]-c[0])*(zr[0]-c[0])+(xr[0]-c[1])*(xr[0]-c[1])); /*calculate the radius of curvature for the neutron path */
/* calculate points along path of neutron through cylinder quit on absorption
* or transmit neutron if 200 points are calculated
* calculate phi for first and last points */
phivec[0]=acos0_2pi((zr[1]-c[0])/nrad,xr[1]-c[1]);phivec[1]=acos0_2pi((zr[2]-c[0])/nrad,xr[2]-c[1]);
neutxypos(&tneuzr,&tneuxr,phivec[0],nrad,c);
/* reset phi[0] and phi[1] to match the length of the slit package rather than cylinder radius*/
if(tneuzr<-splen){
phivec[0]=acos0_2pi((-c[0]-splen)/nrad,-c[1]);
}
neutxypos(&tneuzr,&tneuxr,phivec[1],nrad,c);
if(tneuzr>splen){
phivec[1]=acos0_2pi((-c[0]+splen/2.0)/nrad,-c[1]);
}
dphi=phivec[1]-phivec[0]; /* initial dphi */
idx1=2;
phivec[idx1]=phivec[0]+dphi/2.0; /* calculate center point */
if (!t0checkabsorb(phivec[idx1],nrad,w1,w2,c))
ABSORB;
while (idx1<129){
dphi=phivec[1]-phivec[idx1];
idx1++;
phivec[idx1]=phivec[0]+dphi/2.0;
if (!t0checkabsorb(phivec[idx1],nrad,w1,w2,c))
ABSORB;
if (dphi>0){
while ((phivec[idx1]<phivec[1])&&(idx1<129)){
/* printf("phivec[%i]: %g dphi: %g phivec[1]: %g\n", idx1,phivec[idx1],dphi,phivec[1]);*/
idx1++;
phivec[idx1]=phivec[idx1-1]+dphi;
if (!t0checkabsorb(phivec[idx1],nrad,w1,w2,c))
ABSORB;
}
if (phivec[idx1]>=phivec[1]) idx1--; //remove the point that is beyond phivec[1]
}
else if (dphi<0){
while ((phivec[idx1]>phivec[1])&&(idx1<129)){
/* printf("phivec[%i]: %g\n", idx1,phivec[idx1]);*/
idx1++;
phivec[idx1]=phivec[idx1-1]+dphi;
if (!t0checkabsorb(phivec[idx1],nrad,w1,w2,c))
ABSORB;
}
if (phivec[idx1]<=phivec[1]) idx1--; //remove the point that is beyond phivec[1]
}
else
ABSORB; /* dphi =0? */
}
}
else /* The neutron failed to even hit the chopper */
ABSORB;
%}
MCDISPLAY
%{
double zstep,x1,x2,x3,x4,z1,z2;
int idx, idx2;
line(w2/2.0,ymin,splen,w2/2.0,ymax,splen);
line(w2/2.0,ymin,-splen,w2/2.0,ymax,-splen);
line(-w2/2.0,ymin,splen,-w2/2.0,ymax,splen);
line(-w2/2.0,ymin,-splen,-w2/2.0,ymax,-splen);
line(w2/2.0,ymax,splen,w1/2.0,ymax,0);
line(w1/2.0,ymax,0,w2/2.0,ymax,-splen);
line(-w2/2.0,ymax,splen,-w1/2.0,ymax,0);
line(-w1/2.0,ymax,0,-w2/2.0,ymax,-splen);
line(w2/2.0,ymin,splen,w1/2.0,ymin,0);
line(w1/2.0,ymin,0,w2/2.0,ymin,-splen);
line(-w2/2.0,ymin,splen,-w1/2.0,ymin,0);
line(-w1/2.0,ymin,0,-w2/2.0,ymin,-splen);
circle("zx",0,ymin,0,rad);
circle("zx",0,ymax,0,rad);
zstep=2.0*splen/10.0;
%}
END
|
