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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
|
/*******************************************************************************
*
* 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
|
