/***************************************************************************
*
* Component 3He_simp
*
* Written by Trefor Roberts, ILL, March 1999
*
* Simple cylindrical polarised 3He neutron spin filter cell.
* The glass container for the cell is not included in the model.
*
* Input parameters:
*
* radius: radius of the cylinder (m)
* length: length of the cylinder (m)
* pressure: pressure of the gas in the cell (bar)
* p3he: polarisation of the 3He gas (-1 to +1)
* bx: x component of field at the cell (tesla)
* by: y component of field at the cell (tesla)
* bz: z component of field at the cell (tesla)
*
***************************************************************************/
DEFINE COMPONENT 3He_simp
DEFINITION PARAMETERS (radius,length,pressure,p3he,bx,by,bz)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)
DECLARE
%{
#define opac_cnv 7.33 /* opacity conversion factor so that
the opacity can be expressed in
bar.m.angstroms */
#define gyro 1.832e8 /* absolute value of the gyromagnetic
ratio of the neutron (s-1T-1) */
%}
INITIALIZE
%{
%}
TRACE
%{
double t0,t1; /* time that neutron enters and leaves gas (s) */
double v,lamda; /* neutron velocity and wavelength (ms-1, Anstroms) */
double l_full; /* path length of neutron through gas (m) */
double dt0; /* time neutron spends in the gas (s) */
double opacity; /* opacity of the gas for this neutron (dimless) */
double omega; /* angle through which polarisation precesses
over the path through the cell (radians) */
double px,py,pz; /* components of propagated polarisation (-1 to +1) */
/* calculate the intersection times with the volume of gas, if the neutron
goes through the cell, continue with calculation otherwise all done.
Note that y and z are swapped - this is because the cylindrical axis of
a 3He cell lies along the beam. */
if(cylinder_intersect(&t0,&t1,x,z,y,vx,vz,vy,radius,length))
{
/* Calculate the neutron velocity and wavelength */
v=sqrt(vx*vx+vy*vy+vz*vz);
lamda=2*PI/(V2K*v);
/* Calculate the path length of the neutron through the gas */
dt0=t1-t0;
l_full=v*dt0;
/* Calculate the opacity of the cell for the path length travelled */
opacity=pressure*l_full*lamda*opac_cnv;
/* propagate the polarisation accross the cell (assuming a constant
magnetic field). The actual interaction point is not taken into account
as only the parallel and perpendicular components are important. */
omega=dt0*gyro*sqrt(bx*bx+by*by+bz*bz);
rotate(px,py,pz,sx,sy,sz,omega,bx,by,bz);
sx=px;
sy=py;
sz=pz;
/* adjust the neutron weight according to spin state relative to the
3He nuclei - antiparallel spins are as good as absorbed, whereas parallel
spins are transmitted (depending upon degree of polarisation of gas */
if(scalar_prod(sx,sy,sz,bx,by,bz)>0)
{
p*=0.5*exp(-opacity*(1.0-p3he));
} else {
p*=0.5*exp(-opacity*(1.0+p3he));
}
}
%}
MCDISPLAY
%{
magnify("xyz");
circle("xy",0.0,0.0,-length/2.0,radius);
circle("xy",0.0,0.0,length/2.0,radius);
line(0.0,radius,-length/2.0,0.0,radius,length/2.0);
line(radius,0.0,-length/2.0,radius,0.0,length/2.0);
line(0.0,-radius,-length/2.0,0.0,-radius,length/2.0);
line(-radius,0.0,-length/2.0,-radius,0.0,length/2.0);
%}
END
/* Set polarisation to random value. */
DEFINE COMPONENT pol_set
DEFINITION PARAMETERS ()
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)
TRACE
%{
double norm;
do
{
sx = randpm1();
sy = randpm1();
sz = randpm1();
norm = sx*sx + sy*sy + sz*sz;
} while (norm > 1 || norm == 0);
norm = sqrt(norm);
sx /= norm;
sy /= norm;
sz /= norm;
%}
END
DEFINE COMPONENT pol_monitor
DEFINITION PARAMETERS (nx, ny, filename)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (PSD_N, PSD_p, PSD_p2)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)
DECLARE
%{
int PSD_N[nx][ny];
double PSD_p[nx][ny];
double PSD_p2[nx][ny];
%}
INITIALIZE
%{
int i,j;
for (i=0; i<nx; i++)
for (j=0; j<ny; j++)
{
PSD_N[i][j] = 0;
PSD_p[i][j] = 0;
PSD_p2[i][j] = 0;
}
%}
TRACE
%{
double t0, t1, phi,theta;
int i,j;
phi = atan2(sz,sx);
theta = (sy / (2*sqrt(sx*sx + sy*sy + sz*sz)));
i = floor(nx*(phi/(2*PI) + 0.5));
if(i == nx)
i--; /* Special case for phi = PI. */
else if(i < 0)
y = 0;
j = floor(ny*(theta + 0.5));
if(j == ny)
j--;
else if(j < 0)
j = 0;
PSD_N[i][j]++;
PSD_p[i][j] += p;
PSD_p2[i][j] += p*p;
%}
FINALLY
%{
DETECTOR_OUT_2D(
"Polarisation monitor",
"Longitude [deg]",
"Lattitude [deg]",
-180, 180, -90, 90,
nx, ny,
&PSD_N[0][0],&PSD_p[0][0],&PSD_p2[0][0],
filename);
%}
MCDISPLAY
%{
magnify("");
%}
END
DEFINE INSTRUMENT 3He_test(en,den,pol,pa,p,l)
TRACE
COMPONENT source = Source_flat(
radius = 0.060,
dist = 3.288,
xw = 0.042, yh = 0.082,
E0 = en,
dE = den)
AT (0,0,0) ABSOLUTE
COMPONENT init_pol = pol_set()
AT (0,0,0) RELATIVE source
COMPONENT soller1 = Soller(
xmin = -0.01,
xmax = 0.01,
ymin = -0.01,
ymax = 0.01,
len = 0.09,
divergence = 40)
AT (0,0,0.3) RELATIVE source
COMPONENT soller2 = Soller(
xmin = -0.01,
xmax = 0.01,
ymin = -0.01,
ymax = 0.01,
len = 0.09,
divergence = 40)
AT (0,0,0.4) RELATIVE source ROTATED (0,0,90.0) RELATIVE source
COMPONENT PSD1 = PSD_monitor(
xmin=-0.05,xmax=0.05,ymin=-0.075,ymax=0.075,
nx=51,ny=51,filename="psd1.dat")
AT (0,0,0.501) RELATIVE source
COMPONENT polariser = 3He_simp(
radius=0.03,
length=l,
pressure=p,
p3he=pol,
bx=0,
by=1e-3,
bz=0)
AT (0,0,0.610) RELATIVE source
COMPONENT PSD2 = PSD_monitor(
xmin=-0.05,xmax=0.05,ymin=-0.075,ymax=0.075,
nx=51,ny=51,filename="psd2.dat")
AT (0,0,0.72) RELATIVE source
COMPONENT POL1=pol_monitor(
nx=50,
ny=100,
filename="pol1.dat")
AT (0,0,0.721) RELATIVE source
COMPONENT DIV1=Divergence_monitor(
xmin=-0.05,
xmax=0.05,
ymin=-0.075,
ymax=0.075,
nv=51,
nh=51,
v_maxdiv=2.0,
h_maxdiv=2.0,
filename="div1.dat")
AT (0,0,0.722) RELATIVE source
COMPONENT slit2 = Slit(xmin = -0.01, xmax = 0.01,
ymin = -0.01, ymax = 0.01)
AT (0,0,0.8) RELATIVE source
COMPONENT PSD3 = PSD_monitor(
xmin=-0.05,xmax=0.05,ymin=-0.075,ymax=0.075,
nx=51,ny=51,filename="psd3.dat")
AT (0,0,0.801) RELATIVE source
COMPONENT analyser = 3He_simp(
radius=0.03,
length=0.1,
pressure=100.0,
p3he=pa,
bx=0,
by=1e-3,
bz=0)
AT (0,0,0.855) RELATIVE source
COMPONENT PSD4 = PSD_monitor(
xmin=-0.05,xmax=0.05,ymin=-0.075,ymax=0.075,
nx=51,ny=51,filename="psd4.dat")
AT (0,0,0.910) RELATIVE source
COMPONENT POL2=pol_monitor(
nx=50,
ny=100,
filename="pol2.dat")
AT (0,0,0.911) RELATIVE source
END
Legal notice