/******************************************************************************* * * McXtrace, x-ray tracing package * Copyright, All rights reserved * DTU Physics, Kgs. Lyngby, Denmark * Synchrotron SOLEIL, Saint-Aubin, France * * Component: Air * * %Identification * Written by: M. B. Nielsen * Date: 05.02.2015 * Origin: DTU Fysik, NBI * Release: McXtrace 1.4 * * Component simulating atmospheric air. * * %Description * 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, but it may OR may NOT be negligible. * I have not yet looked into this last question, so I can't say if the lack of absorption is a bad thing or if it is allowable. * * 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 * * %Parameters * * 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: [0-1] Fraction of rays to scatter from the air * pressure: [Pa] Total pressure of the air gas * temperature:[K] Absolute temperature * * %End *******************************************************************************/ DEFINE COMPONENT Air /* number-density (num_den) is declared and computed in TRACE !!! */ SETTING PARAMETERS (frac=0.3, pressure=101325, temperature=273.15+21.1, R_gas=8.3144621, bond_N=1.0976, bond_O=1.2074, Nitrogen_part=0.781, Oxygen_part=0.21, Argon_part=0.009, 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) DECLARE %{ int shape; %} 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) { fprintf (stderr, "ERROR: %s: Air volume has invalid dimensions, please check parameter values.\n", NAME_CURRENT_COMP); exit (-1); } /* 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 ("INFO %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; /*material constants for the constituent gases*/ double an[4] = { 12.2126, 3.1322, 2.0125, 1.1663 }; double ao[4] = { 3.0485, 2.2868, 1.5463, 0.867 }; double aar[4] = { 7.4845, 6.7723, 0.6539, 1.6442 }; double bn[4] = { 0.0057, 9.8933, 28.9975, 0.5826 }; double bo[4] = { 13.2771, 5.7011, 0.3239, 32.9089 }; double bar[4] = { 0.9072, 14.8407, 43.8983, 33.3929 }; double cn = -11.529; double co = 0.2508; double car = 1.4445; /* 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-polynomium has units of squared Angstrom 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) constitues 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 Angstrom - 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) { p *= Ptot_sc_unpolarized / 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; aim_y = target_y; aim_z = target_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_real (&kx, &ky, &kz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP, x, y, z, 2); } 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