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
222
223
224
225
226
227
228
229
230
231
|
/*******************************************************************************
*
* McXtrace, X-ray tracing package
* Copyright, All rights reserved
* DTU Physics, Kgs. Lyngby, Denmark
* Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Source_flat
*
* %Identification
* Written by:Erik Knudsen
* Date: September 25, 2009
* Origin: Risoe
* Release: McXtrace 0.1_alpha
*
* A flat rectangular or circular surface emitting x-rays
*
* %Description
* A circular or rectangular xray source. Spectrum may be either gaussian or uniform around a central wavelength/energy
* or read from a datafile. Xrays are considered emitted uniformly into 4pi, but a square target retricts the beam to
* that window and scales the beam intensity accordingly.
* If an input spectrum datafile (spectrum_file) is not specified, the beam is restricted to emit photons between E0+-dE keV, or lambda0+-dlambda AA, whichever is given.
* The input spectrum file should be formatted such that x-ray energy/wavelength is in the first column and the intensity in the second. Any preceding
* lines starting with # are considered part of the file header. If a datafile is given, a nonzero E0 value indicates that is is parametrized by energy (in keV)
* as opposed to wavelength (in AA). Wavelength is the default.
* Flux is set in the unit photons/s
*
* Example: Source_flat(xwidth=1e-3,yheight=1e-3, focus_xw=0.5e-2, focus_yh=0.45e-2,dist=1, E0=E0, dE=DE)
*
* %Parameters
* radius: [m] Radius of circle in (x,y,0) plane where x-rays are generated.
* yheight: [m] Height of rectangle in (x,y,0) plane where x-rays are generated.
* xwidth: [m] Width of rectangle in (x,y,0) plane where x-rays are generated. Overrides xmin and xmax.
* xmax: [m] Upper bound of x-interval where photons are generated.
* xmin: [m] Lower bound of x-interval where photons are generated.
* dist: [m] Distance to target along z axis.
* focus_xw: [m] Width of target
* focus_yh: [m] Height of target
* E0: [keV] Mean energy of xrays.
* dE: [keV] Energy half spread of x-rays (flat or gaussian sigma).
* lambda0: [AA] Mean wavelength of x-rays.
* dlambda: [AA] Wavelength half spread of x-rays.
* gauss: [1] Gaussian (1) or Flat (0) energy/wavelength distribution
* flux: [pht/s] Total flux radiated from the source
* randomphase: [ ] If nonzero, the phase of the emotted photon is random, i.e. source is fully incoherent. otherwise the value of phase is used.
* phase: [rad] Set phase to something given.
* spectrum_file: [string] Filename for optional spectrum-file
* %End
*******************************************************************************/
DEFINE COMPONENT Source_flat
SETTING PARAMETERS (radius=0, yheight=0, xwidth=0,xmin=0,xmax=0,
dist=0, focus_xw=.045, focus_yh=.12,
E0=0, dE=0, lambda0=0, dlambda=0, flux=0,gauss=0,randomphase=1,phase=0,
string spectrum_file="")
/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */
SHARE
%{
%include "read_table-lib"
%}
DECLARE
%{
double pmul;
double srcArea;
int square;
t_Table spectrum_T;
%}
INITIALIZE
%{
square = 0;
srcArea=0;
/* Determine source area */
if (xmin && xmax && !xwidth){
xwidth=xmax-xmin;
}
if (radius && !yheight && !xwidth ) {
square = 0;
srcArea = PI*radius*radius;
} else if(yheight && xwidth) {
square = 1;
srcArea = xwidth * yheight;
}
if (srcArea <= 0) {
printf("Source_flat: %s: Source area is <= 0 !\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(-1);
}
if (dist <= 0 || focus_xw <= 0 || focus_yh <= 0) {
printf("Source_flat: %s: Target area unmeaningful! (negative dist / focus_xw / focus_yh)\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(-1);
}
if (spectrum_file && strlen(spectrum_file)>0 && strcmp(spectrum_file,"NULL")!=0 ){
/*read spectrum from file*/
int status=0;
if ( (status=Table_Read(&(spectrum_T),spectrum_file,0))==-1){
fprintf(stderr,"Source_flat(%s) Error: Could not parse file \"%s\"\n",NAME_CURRENT_COMP,spectrum_file?spectrum_file:"");
exit(-1);
}
/*data is now in table spectrum_T*/
/*integrate to get total flux, assuming numbers have been corrected for measuring aperture*/
int i;
double pint=0;
t_Table *T=&(spectrum_T);
for (i=0;i<spectrum_T.rows-1;i++){
pint+=((T->data[i*T->columns+1]+T->data[(i+1)*T->columns+1])/2.0)*fabs(T->data[(i+1)*T->columns]-T->data[i*T->columns]);
}
printf("Source_flat(%s) Integrated intensity radiated is %g pht/s\n",NAME_CURRENT_COMP,pint);
if(E0) printf("Source_flat(%s) E0!=0 -> assuming intensity spectrum is parametrized by energy [keV]\n",NAME_CURRENT_COMP);
} else if (!E0 && !lambda0){
fprintf(stderr,"Source_flat(%s): Error: Must specify either wavelength or energy distribution\n",NAME_CURRENT_COMP);
exit(-1);
}
/*calculate the X-ray weight from the flux*/
if (flux){
pmul=flux;
}else{
pmul=1;
}
pmul*=1.0/((double) mcget_ncount());
%}
TRACE
%{
double chi,e,k,l,r, xf, yf, rf, dx, dy, pdir;
phi=0;
z=0;
if (square == 1) {
if (xmin && xmax){
x = xmin + xwidth*rand01();
}else{
x = xwidth * (rand01() - 0.5);
}
y = yheight * (rand01() - 0.5);
} else {
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
}
randvec_target_rect_real(&xf, &yf, &rf, &pdir,
0, 0, dist, focus_xw, focus_yh, ROT_A_CURRENT_COMP, x, y, z, 2);
dx = xf-x;
dy = yf-y;
rf = sqrt(dx*dx+dy*dy+dist*dist);
/*pdir contains the unnormalized solid angle weighting */
p = pmul*pdir/(4*M_PI);
if (spectrum_file && strlen(spectrum_file) && strcmp(spectrum_file,"NULL")!=0){
double pp=0;
//while (pp<=0){
l=spectrum_T.data[0]+ (spectrum_T.data[(spectrum_T.rows-1)*spectrum_T.columns] -spectrum_T.data[0])*rand01();
pp=Table_Value(spectrum_T,l,1);
//}
p*=pp;
/*if E0!=0 the tabled value is assumed to be energy in keV*/
if (E0!=0){
k=E2K*l;
}else{
k=(2*M_PI/l);
}
}else if (E0){
if(!dE){
e=E0;
}else if (gauss){
e=E0+dE*randnorm();
}else{
e=randpm1()*dE + E0;
}
k=E2K*e;
}else if (lambda0){
if (!dlambda){
l=lambda0;
}else if (gauss){
l=lambda0+dlambda*randnorm();
}else{
l=randpm1()*dlambda*0.5 + lambda0;
}
k=(2*M_PI/l);
}
ky=k*dy/rf;
kx=k*dx/rf;
kz=k*dist/rf;
/*randomly pick phase or set to something real*/
if (randomphase){
phi=rand01()*2*M_PI;
}else{
phi=phase;
}
/*set polarization vector*/
Ex=0;Ey=0;Ez=0;
%}
FINALLY
%{
Table_Free(&(spectrum_T));
%}
MCDISPLAY
%{
if (square == 1) {
rectangle("xy",0,0,0,xwidth,yheight);
} else {
circle("xy",0,0,0,radius);
}
if (dist) {
dashed_line(0,0,0, -focus_xw/2,-focus_yh/2,dist, 4);
dashed_line(0,0,0, focus_xw/2,-focus_yh/2,dist, 4);
dashed_line(0,0,0, focus_xw/2, focus_yh/2,dist, 4);
dashed_line(0,0,0, -focus_xw/2, focus_yh/2,dist, 4);
}
%}
END
|
