/*******************************************************************************
*
* McXtrace, Xray-tracing package
* Copyright 1997-2003, All rights reserved
* Technical University of Denmark, Kgs. Lyngby, Denmark
* Institut Laue Langevin, Grenoble, France
* University of Copenhagen, Copenhagen, Denmark
*
* Component: Air
*
* %I
* Written by: M. B. Nielsen
* Date: 05.02.2015
* Version: $Revision: beta$
* Origin: DTU Fysik, NBI
*
* Component simulating atmospheric air.
*
* %D
* The component simulates air and can be inserted as if it was just some extra sample placed somewhere in the beam line.
* The air component is intended to be used in all kinds of setups where air may introduce background.
* Code structure in this component is based on the component Saxs_spheres.
* The shape of the sample may be a filled box with dimensions
* xwidth, yheight, zdepth, a filled cylinder with dimensions radius and yheight,
* a filled sphere with radius R.
* (NB: As we assume air to be an ideal gas, the volume fractions of the elements in the gas are merely the mole fractional part of the given element.
* From this the number density of atoms/molecules is calculated)
* The air is dry and assumed to be made of nitrogen, oxygen and argon - all other constituents are neglected.
*
* So far the calculations of the scattering probability (and hence also the weight multiplier) assumes the x-ray source to be unpolarized.
* Further the component does not yet account for absorption of x-rays. I.e. absorption is simply omitted.
*
* Example:
* COMPONENT air1 = Air(
* frac = 0.4, pressure = 50000, temperature = 270,
* xwidth = 0.5, yheight = 0.5, zdepth = 1.5, target_index = 1,
* focus_xw = 0.5, focus_yh = 0.5)
* AT (0, 0, 15) RELATIVE Origin
*
* %P
*
* xwidth: [m] Width of the air volume.
* yheight: [m] Height of the air volume.
* zdepth: [m] Depth of the air volume.
* radius: [m] Radius of spherical or cylindrical air volume.
* target_x:[m] X-coordinate of sampling window.
* target_y:[m] Y-coordinate of sampling window.
* target_z:[m] Z-coordinate of sampling window.
* target_index: [ ] index of target component putting sampling window on a subsequent component.
* focus_xw:[m] Width of the sampling window.
* focus_yh:[m] Height of the sampling window.
* focus_ah:[rad] Vertical (height) opening angle of sampling window.
* focus_aw:[rad] Horizontal (width) opening angle of sampling window.
* focus_r:[rad] Radius of circular sampling window.
* frac: fraction of rays to scatter from the air [ ]
* pressure: total pressure of the air gas [Pa]
* temperature: absolute temperature [K]
*
* %Link
* %E
*******************************************************************************/
DEFINE COMPONENT Air
DEFINITION PARAMETERS (R_gas=8.3144621, bond_N=1.0976, bond_O=1.2074, Nitrogen_part=0.781, Oxygen_part=0.21, Argon_part=0.009)
/* number-density (num_den) is declared and computed in TRACE !!! */
SETTING PARAMETERS (frac=0.3, pressure=101325, temperature=273.15+21.1,
xwidth=0, yheight=0, zdepth=0, radius=0,
target_x=0, target_y=0, target_z=6, int target_index=0,
focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, focus_r=0)
OUTPUT PARAMETERS (shape,an,bn,cn,ao,bo,co,aar,bar,car)
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
DECLARE
%{
double shape;
double an[4]={12.2126, 3.1322, 2.0125, 1.1663}, ao[4]={3.0485, 2.2868, 1.5463, 0.867}, aar[4]={7.4845, 6.7723, 0.6539, 1.6442};
double bn[4]={0.0057, 9.8933, 28.9975, 0.5826}, bo[4]={13.2771, 5.7011, 0.3239, 32.9089}, bar[4]={0.9072, 14.8407, 43.8983, 33.3929};
double cn=-11.529, co=0.2508, car=1.4445;
%}
INITIALIZE
%{
shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere */
if (xwidth && yheight && zdepth) shape=1; /* box */
else if (radius > 0 && yheight) shape=0; /* cylinder */
else if (radius > 0 && !yheight) shape=2; /* sphere */
if (shape < 0)
exit(fprintf(stderr,"Air: %s: sample has invalid dimensions.\n"
"ERROR Please check parameter values.\n", NAME_CURRENT_COMP));
/* now compute target coords if a component index is supplied */
if (!target_index && !target_x && !target_y && !target_z)
target_index=1;
if (target_index) {
Coords ToTarget;
ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index), POS_A_CURRENT_COMP);
ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
coords_get(ToTarget, &target_x, &target_y, &target_z);
}
if (!(target_x || target_y || target_z)) {
printf("Air: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP);
target_z=1;
}
%}
TRACE
%{
double l0, l1, l_full, l, dl, ran_var;
double aim_x=0, aim_y=0, aim_z=1, axis_x, axis_y, axis_z;
double f0, solid_angle, sc_angle, kx_i, ky_i, kz_i, k, q, qx, qy, qz, polarization, Ptot_sc_unpolarized;
char intersect=0;
/* Intersection X_ray trajectory / sample (sample surface) */
if (shape == 0)
intersect = cylinder_intersect(&l0, &l1, x, y, z, kx, ky, kz, radius, yheight);
else if (shape == 1)
intersect = box_intersect(&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth);
else if (shape == 2)
intersect = sphere_intersect(&l0, &l1, x, y, z, kx, ky, kz, radius);
if (intersect) {
if (l0 < 0) {
fprintf(stderr,"photon already inside sample %s - absorbing\n",NAME_CURRENT_COMP);
ABSORB;
}
/* Xray enters at l=l0. */
l_full = (l1 - l0); /* Length of full path through the air */
k = sqrt(kx*kx + ky*ky + kz*kz);
dl = rand01()*(l1 - l0) + l0; /* Time of scattering */
PROP_DL(dl); /* Point of scattering */
l = (dl-l0); /* Penetration in the air */
/* The total probability for the (unpolarized) X-ray to interact with the air as a function of k is
approximated with a 9th degree polynomium to the logarithm of the total probability. This is done
because the real expression had to be solved numerically. The exp-polynomial has units of squared
AAngstrom and to become the actual total scattering probabillity it has to be multiplied by
l_full*P*NA/(R*T). The factor of P*NA/(R*T) does however have units of reciprocal cubed meters so we have
to scale by 10^-30 [m^2/AA^2] to match the polynomium. NB: The factor of (n_i/n)*P*NA/(R*T) constitutes
the particle number density; while the factor of (n_i/n) is contained in the polynomium, the factor of
P*NA/(R*T) has to be multiplied below. The variable l_full has units of meters so it too has to be
rescaled to AAngstroms - hence we multiply by 10^10 [AA/m].*/
Ptot_sc_unpolarized = exp(-13.239209120968501922 + 7.7523973158902582507e-1*k - 1.9753647234751300397*k*k + 1.1395382484098412215*k*k*k - 3.5117800398913380701e-1*k*k*k*k + 6.493474609841591379e-2*k*k*k*k*k - 7.40740935645497008e-3*k*k*k*k*k*k + 5.1039561005291247e-4*k*k*k*k*k*k*k - 1.947399876590436e-5*k*k*k*k*k*k*k*k + 3.1587814532171995885e-7*k*k*k*k*k*k*k*k*k)*l_full*pressure*NA/(R_gas*temperature)*1e-30*1e10;
/* The mc-choice to scatter from the air or not is governed by the variable 'frac'*/
/* NB: The variable 'Ptot_sc' has to be evaluated BEFORE we make the mc-choice because we also need to be able to weight the non-interacting rays.*/
/* This non-interaction weight is p *= (1 - Ptot_sc_unpolarized)/(1 - frac)*/
ran_var = rand01();
if (ran_var < frac) {
kx_i=kx;
ky_i=ky;
kz_i=kz;
/* Vector pointing at target (anal./det.) */
if ((target_x || target_y || target_z)) {
aim_x = target_x-x;
aim_y = target_y-y;
aim_z = target_z-z;
}
if (focus_aw && focus_ah)
randvec_target_rect_angular(&kx, &ky, &kz, &solid_angle, aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
else if (focus_xw && focus_yh)
randvec_target_rect(&kx, &ky, &kz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP);
else
randvec_target_circle(&kx, &ky, &kz, &solid_angle, aim_x, aim_y, aim_z, focus_r);
NORM(kx, ky, kz);
kx *= k;
ky *= k;
kz *= k;
qx = (kx_i-kx);
qy = (ky_i-ky);
qz = (kz_i-kz);
q = sqrt(qx*qx+qy*qy+qz*qz);
sc_angle = acos( (kx_i*kx + ky_i*ky + kz_i*kz)/(k*k) );
polarization = (1 + cos(sc_angle)*cos(sc_angle))/2;
/* Now we pick out which type of particle in the gas to scatter from assuming that air behaves as an ideal gas.
This approximation is valid only when 0 < q < 25 [AA^-1] and hence 0