/******************************************************************************* * * McXtrace, x-ray tracing package * Copyright, All rights reserved * DTU Physics, Kgs. Lyngby, Denmark * Synchrotron SOLEIL, Saint-Aubin, France * * Component: Grating_trans * * %Identification * * Written by: Erik B Knudsen (erkn@fysik.dtu.dk) * Date: December 2016 * Version: 1.0 * Release: McXtrace 1.4 * Origin: DTU Physics * * Transmission grating * * %Description * Model of a 1D rectangular transmission grating based on the theory developed in Schnopper et. al., Applied Optics, 1977. * The grating lines are assumed to be vertical. Within each period a fraction gamma is the "open" fraction. * (I.e. 1 is completely open). At present only absorption in the substrate (modelled by the thickness sdepth) * is included. * * This component is currently undergoing validation. * * Example: Grating_trans( * xwidth=25e-3, yheight=25e-3, gamma=0.4, period=2000e-10, zdepth=5100e-10, max_order=3, material="Au.txt") * * %Parameters * Input parameters: * period: [m] Distance between grating grooves. * zdepth: [m] Depth of grooves. * gamma: [0-1] Ratio between groove and period aka duty cycle. 1 means fully open. * sdepth: [m] Thickness of substrate. The default is to have no substrate - i.e. rods. * xwidth: [m] Width of the grating. Defines how many lines there are in total. * yheight: [m] Height of the grating. * material: [str] Data file containing the material from which the grating is made. * substrate: [str] Data file containing material data for the substrate. * fixed_delta:[0/1] Set delta to the given constant. Useful for debugging. * max_order: [1] Maximum order to diffract * * %End *******************************************************************************/ DEFINE COMPONENT Grating_trans SETTING PARAMETERS (xwidth=1e-3, yheight=1e-3, period=1e-6,gamma=0.5,zdepth=1e-6,sdepth=0, string material="Au.txt", string substrate="",int max_order=2,fixed_delta=0) SHARE %{ %include "read_table-lib" %} DECLARE %{ int Z; int mu_c; double Ar; double rho; double delta_prefactor; t_Table table; double srho; int smu_c; t_Table stable; int order; %} INITIALIZE %{ int status; if (xwidth == 0 || yheight == 0) { fprintf (stderr, "Error (%s): Grating has 0 area.\n", NAME_CURRENT_COMP); exit (-1); } if (zdepth == 0) { fprintf (stderr, "Error (%s): Grating has no grooves (zdepth==0).\n", NAME_CURRENT_COMP); exit (-1); } /*check if material datafiles are present - if so load them*/ if ((status = Table_Read (&(table), material, 0)) == -1) { fprintf (stderr, "Error: Could not parse file \"%s\" in COMP %s\n", material, NAME_CURRENT_COMP); exit (-1); } char** header_parsed; header_parsed = Table_ParseHeader (table.header, "Z", "A[r]", "rho", "Z/A", "sigma[a]", NULL); if (header_parsed[0]) { Z = strtol (header_parsed[0], NULL, 10); } if (header_parsed[1]) { Ar = strtod (header_parsed[1], NULL); } if (header_parsed[2]) { rho = strtod (header_parsed[2], NULL); } else { fprintf (stderr, "Warning(%s): %s not found in header of %s, set to 1\n", NAME_CURRENT_COMP, "rho", material); rho = 1; } /*which columns holds the mus*/ mu_c = 5; if (table.columns == 3) mu_c = 1; delta_prefactor = NA * (rho * 1e-24) / Ar * 2.0 * M_PI * RE; /*read substrate datafile if present*/ if (substrate && strlen (substrate)) { if ((status = Table_Read (&(stable), substrate, 0)) == -1) { fprintf (stderr, "Error(%s): Could not parse file \"%s\".\n", NAME_CURRENT_COMP, material); exit (-1); } header_parsed = Table_ParseHeader (stable.header, "Z", "A[r]", "rho", "Z/A", "sigma[a]", NULL); if (header_parsed[2]) { srho = strtod (header_parsed[2], NULL); } smu_c = 5; if (stable.columns == 3) smu_c = 1; } %} TRACE %{ double e, k, f, delta, beta, mu0, smu0, pmul, order_factor, theta, order_likelihood; double thx, thy, period_eff, zdepth_eff; int M; // order; PROP_Z0; if (x < -xwidth * 0.5 || x > xwidth * 0.5 || y < -yheight * 0.5 || y > yheight * 0.5) { RESTORE_XRAY (INDEX_CURRENT_COMP, x, y, z, kx, ky, kz, phi, t, Ex, Ey, Ez, p); } else { /*get the index of refraction*/ k = sqrt (kx * kx + ky * ky + kz * kz); e = k * K2E; /*Material's Number Density of Electrons [e/A^3] incl f' scattering length correction*/ /*We have reparametrized e as log(e) for a more practical constant step table.*/ f = Table_Value (table, e, 1); delta = (fixed_delta ? fixed_delta : f / (k * k) * delta_prefactor); mu0 = Table_Value (table, e, mu_c) * rho * 1e2; beta = mu0 / (2 * k * 1e10); /*pick an order according the relative intensity as given by Schnopper*/ /* correct the period for horizontal angle*/ thx = atan (fabs (kx / kz)); period_eff = cos (thx) * period; M = xwidth / period_eff; order = floor (rand01 () * (max_order * 2 + 1)) - max_order; /* n_m(q)/N_0(q)=([sin(Mmpi)/Msin(mpi)]^2 [sin( (a/d)mpi)/mpi]^2 [1+exp(-2qzk)-2exp(-qzk)cos(qzdelta)]*/ /* M=number of wires in the grating, * d grating spacing (period) * a width of grating opening (duty-cycle) * z the thickness of grating wire * k imaginary part of refr index * n=1-delta real part of refr index * q wavenumber. */ /* correct the period for horizontal angle*/ thx = atan (fabs (kx / kz)); period_eff = cos (thx) * period; M = xwidth / period_eff; double tot = 0; int i; for (i = 0; i < max_order; i++) { if (i) { order_factor = pow (sin (gamma * i * M_PI) / (i * M_PI), 2.0); order_likelihood = order_factor * (1 + exp (-2.0 * k * 1e10 * zdepth * beta) - 2.0 * exp (-k * 1e10 * zdepth * beta) * cos (k * 1e10 * zdepth * delta)); } else { order_likelihood = gamma * gamma + (1 - gamma) * (1 - gamma) * exp (-2.0 * k * 1e10 * zdepth * beta) - 2.0 * gamma * (1 - gamma) * exp (-k * 1e10 * zdepth * beta) * cos (k * 1e10 * zdepth * delta); } tot += order_likelihood; } /* correct zdepth for vertical angle*/ thy = atan (fabs (ky / kz)); zdepth_eff = zdepth / cos (thy); if (order) { order_factor = pow (sin (gamma * order * M_PI) / (order * M_PI), 2.0); order_likelihood = order_factor * (1 + exp (-2.0 * k * 1e10 * zdepth_eff * beta) - 2.0 * exp (-k * 1e10 * zdepth_eff * beta) * cos (k * 1e10 * zdepth_eff * delta)); } else { order_likelihood = gamma * gamma + (1 - gamma) * (1 - gamma) * exp (-2.0 * k * 1e10 * zdepth_eff * beta) - 2.0 * gamma * (1 - gamma) * exp (-k * 1e10 * zdepth_eff * beta) * cos (k * 1e10 * zdepth_eff * delta); } double fmc; fmc = 1.0 / (max_order * 2 + 1); pmul = order_likelihood / fmc; /*Pick an angle centered on the order*/ theta = order * 2 * M_PI / (k * 1e10 * period_eff); /*the 1e10 factor because k is in AA^-1 and period is in m*/ /* we should rotate around an axis perpendicular to k and in yz-plane*/ /*i.e. (ay * ky + az*kz ) = 0 && ay>0 => Without loss of gen. ay=1 => az= -ky/kz*/ double az; az = -ky / kz; /*change direction*/ rotate (kx, ky, kz, kx, ky, kz, theta, 0, 1, az); /*do absorption according to mean absorption coefficient in the "slab", remembering we can be at an angle. Just use the angle theta.*/ smu0 = Table_Value (stable, e, smu_c) * srho * 1e2; p *= pmul * exp (-smu0 * (sdepth) / cos (theta)) * exp (-mu0 * (zdepth * (1 - gamma))); SCATTER; } %} MCDISPLAY %{ int i; rectangle ("xy", 0, 0, -(sdepth) / 2.0, xwidth, yheight); rectangle ("xy", 0, 0, (sdepth) / 2.0, xwidth, yheight); for (i = -1; i < 2; i++) { box (period / 4.0 + i * period, 0, sdepth / 2.0 + zdepth / 2.0, period / 2.0, yheight, zdepth, 0, 0, 1, 0); box (-3 * period / 4.0 + i * period, 0, sdepth / 2.0 + zdepth / 2.0, period / 2.0, yheight, zdepth, 0, 0, 1, 0); } %} END