1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion=6.1 |
---|
3 | |
---|
4 | Function TestSmear_2D() |
---|
5 | |
---|
6 | Variable DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA |
---|
7 | DX = 0.5 |
---|
8 | num = 128 |
---|
9 | // x0 = 64 |
---|
10 | x0 = 114 |
---|
11 | y0 = 64 |
---|
12 | L1 = 300 //units of cm ?? |
---|
13 | L2 = 130 |
---|
14 | s1 = 5.0/2 |
---|
15 | s2 = 1.27/2 |
---|
16 | sig_det = 0.5 //not sure about this |
---|
17 | dlamb = 0.15 |
---|
18 | lambda = 6 |
---|
19 | |
---|
20 | Duplicate/O root:no_gravity_dat:no_gravity_dat_mat root:no_gravity_dat:John_mat |
---|
21 | Wave data=root:no_gravity_dat:John_mat |
---|
22 | |
---|
23 | SUB_SMEAR_2D(DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA,DATA) |
---|
24 | |
---|
25 | Duplicate/O root:no_gravity_dat:John_mat root:no_gravity_dat:John_mat_log |
---|
26 | Wave log_data=root:no_gravity_dat:John_mat_log |
---|
27 | |
---|
28 | log_data = log(data) |
---|
29 | |
---|
30 | end |
---|
31 | |
---|
32 | // I have changed the array indexing to [0,], so subtract 1 from the x0,Y0 center |
---|
33 | // to shift from detector coordinates to Igor array index |
---|
34 | // |
---|
35 | // |
---|
36 | // !! the wi values do not match what is written in John's notebook. Are these the |
---|
37 | // correct values for hermite integration?? |
---|
38 | // |
---|
39 | Function SUB_SMEAR_2D(DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA,DATA) //,Q_MODEL_NAME) |
---|
40 | Variable DX,NUM,X0,Y0,L1,L2,S1,S2,SIG_DET,DLAMB,LAMBDA |
---|
41 | Wave data |
---|
42 | |
---|
43 | Variable I,J,KI,KJ //integers |
---|
44 | Variable SUMm,THET0,Q0,R_PL,R_PD,Q0_PL,Q0_PD,LP,V_R,V_L |
---|
45 | Variable PHI,R0,SIGQ_R,SIGQ_A,Q_PL,Q_PD,DIF_PD_I |
---|
46 | Variable RES_I,RES_J,RES,DIF_PL_J,DIF_PD_J,DIF_PL_I |
---|
47 | // DIMENSION DATA(128,128),XI(3),WI(3) |
---|
48 | // EXTERNAL Q_MODEL_NAME |
---|
49 | // PARAMETER PI = 3.14159265 |
---|
50 | Variable N_QUAD = 3 |
---|
51 | Make/O/D xi_h = {.6167065887,1.889175877,3.324257432} |
---|
52 | Make/O/D wi_h = {.72462959522,.15706732032,.45300099055E-2} |
---|
53 | |
---|
54 | //C DATA XI/.4360774119,1.3358490740,2.3506049736/ |
---|
55 | // DATA XI/.6167065887,1.889175877,3.324257432/ |
---|
56 | // DATA WI/.72462959522,.15706732032,.45300099055E-2/ |
---|
57 | //C DX : PIXEL SIZE, CM |
---|
58 | //C NUM: NUMBER OF PIXELS ACROSS DETECTOR. (128) |
---|
59 | //C X0,Y0: BEAM CENTER, IN UNITS OF PIXELS. |
---|
60 | //C L1: SOURCE TO SAMPLE DISTANCE. |
---|
61 | //C L2: SAMPLE TO DETECTOR DISTANCE. |
---|
62 | //C S1: SOURCE APERTURE RADIUS. |
---|
63 | //C S2: SAMPLE APERTURE RADIUS. |
---|
64 | //C SIG_DET:STANDARD DEVIATION OF DETECTOR SPATIAL RESOLUTION. |
---|
65 | //C DLAMB: FWHM WAVLENGTH RESOLUTION. |
---|
66 | //C LAMBDA: MEAN WAVELENGTH. |
---|
67 | //C DATA: OUTPUT SMEARED ARRAY (NUM,NUM) |
---|
68 | |
---|
69 | Make/O/D/N=(128,128) sigQR, sigQA |
---|
70 | |
---|
71 | |
---|
72 | LP = 1 / ( 1/L1 + 1/L2 ) |
---|
73 | V_R = 0.25*(S1/L1)^2 + 0.25*(S2/LP)^2 + (SIG_DET/L2)^2 |
---|
74 | V_L = DLAMB^2/6. |
---|
75 | for(i=0;i<num;i+=1) |
---|
76 | R_PL = DX*(I-X0) |
---|
77 | for(j=0;j<num;j+=1) |
---|
78 | R_PD = DX*(J-Y0) |
---|
79 | PHI = ATAN(R_PD/R_PL) //do I need atan2 here? |
---|
80 | R0 = SQRT(R_PL^2+R_PD^2) |
---|
81 | THET0 = ATAN(R0/L2) |
---|
82 | Q0 = 4*PI*SIN(0.5*THET0)/LAMBDA |
---|
83 | //C DETERMINE Q VECTOR, CARTESIAN REPRESENTATION. |
---|
84 | Q0_PL = Q0*COS(PHI) |
---|
85 | Q0_PD = Q0*SIN(PHI) |
---|
86 | //C DETERMINE SIGMA'S FOR RESOLUTION FUNCTION, RADIALLY, AZIMUTHAL |
---|
87 | SIGQ_R = Q0*SQRT(V_R+V_L) |
---|
88 | SIGQ_A = Q0*SQRT(V_R) |
---|
89 | |
---|
90 | sigQR[i][j] = sigq_R |
---|
91 | sigQA[i][j] = sigq_A |
---|
92 | |
---|
93 | SUMm = 0.0 |
---|
94 | for(KI=0;ki<N_quad;ki+=1) |
---|
95 | DIF_PL_I = SIGQ_R*COS(PHI)*xi_h[ki] |
---|
96 | DIF_PD_I = SIGQ_R*SIN(PHI)*xi_h[ki] |
---|
97 | for( KJ=0;kj<N_QUAD;kj+=1) |
---|
98 | DIF_PL_J = SIGQ_A*SIN(PHI)*xi_h[kj] |
---|
99 | DIF_PD_J = SIGQ_A*COS(PHI)*xi_h[kj] |
---|
100 | //C -,- |
---|
101 | Q_PL = Q0_PL - DIF_PL_I - DIF_PL_J |
---|
102 | Q_PD = Q0_PD - DIF_PD_I - DIF_PD_J |
---|
103 | SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD) //,Q_MODEL_NAME) |
---|
104 | //C -,+ |
---|
105 | Q_PL = Q0_PL - DIF_PL_I + DIF_PL_J |
---|
106 | Q_PD = Q0_PD - DIF_PD_I + DIF_PD_J |
---|
107 | SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD) //,Q_MODEL_NAME) |
---|
108 | //C +,- |
---|
109 | Q_PL = Q0_PL + DIF_PL_I - DIF_PL_J |
---|
110 | Q_PD = Q0_PD + DIF_PD_I - DIF_PD_J |
---|
111 | SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD) //,Q_MODEL_NAME) |
---|
112 | //C +,+ |
---|
113 | Q_PL = Q0_PL + DIF_PL_I + DIF_PL_J |
---|
114 | Q_PD = Q0_PD + DIF_PD_I + DIF_PD_J |
---|
115 | SUMm = SUMm + wi_h[ki]*wi_h[kj]*I_MACRO(Q_PL,Q_PD) //,Q_MODEL_NAME) |
---|
116 | endfor // KJ |
---|
117 | endfor // KI |
---|
118 | DATA[i][j] = SUMm / PI |
---|
119 | endfor // J |
---|
120 | endfor // I |
---|
121 | |
---|
122 | RETURN(0) |
---|
123 | |
---|
124 | END |
---|
125 | |
---|
126 | /// --- either way, same to machine precision |
---|
127 | Function I_MACRO(Q_PL,Q_PD) //,Q_MODEL_NAME) |
---|
128 | Variable Q_PL,Q_PD |
---|
129 | |
---|
130 | Variable I_MACRO,Q,PHI,PHI_MODEL,NU |
---|
131 | |
---|
132 | //Q_MODEL_NAME |
---|
133 | //eccentricity factor for ellipse in John's code... |
---|
134 | NU = 1 |
---|
135 | |
---|
136 | //C PHI = ATAN(Q_PD/Q_PL) |
---|
137 | |
---|
138 | Q = SQRT((NU*Q_PD)^2+Q_PL^2) |
---|
139 | |
---|
140 | WAVE cw = $"root:coef_Peak_Gauss" |
---|
141 | |
---|
142 | I_MACRO = Peak_Gauss_modelX(cw,Q) |
---|
143 | // I_MACRO = Q_MODEL_NAME(Q) |
---|
144 | |
---|
145 | RETURN(I_MACRO) |
---|
146 | END |
---|
147 | |
---|
148 | //Function I_MACRO(Q_PL,Q_PD) //,Q_MODEL_NAME) |
---|
149 | // Variable Q_PL,Q_PD |
---|
150 | // |
---|
151 | // Variable I_MACRO |
---|
152 | // |
---|
153 | // Make/O/D/N=1 fcnRet,xptW,yPtw |
---|
154 | // xptw[0] = q_pl |
---|
155 | // yptw[0] = q_pd |
---|
156 | // |
---|
157 | // WAVE cw = $"root:coef_sf" |
---|
158 | // |
---|
159 | // I_MACRO = Sphere_2DX(cw,xptw,yptw) |
---|
160 | // |
---|
161 | // RETURN(I_MACRO) |
---|
162 | //END |
---|
163 | |
---|
164 | ////Structure ResSmear_2D_AAOStruct |
---|
165 | //// Wave coefW |
---|
166 | //// Wave zw //answer |
---|
167 | //// Wave qy // q-value |
---|
168 | //// Wave qx |
---|
169 | //// Wave qz |
---|
170 | //// Wave sigQx //resolution |
---|
171 | //// Wave sigQy |
---|
172 | //// Wave fs |
---|
173 | //// String info |
---|
174 | ////EndStructure |
---|
175 | // |
---|
176 | Function Smear_2DModel_5(fcn,s) |
---|
177 | FUNCREF SANS_2D_ModelAAO_proto fcn |
---|
178 | Struct ResSmear_2D_AAOStruct &s |
---|
179 | |
---|
180 | String weightStr="gauss5wt",zStr="gauss5z" |
---|
181 | Variable nord=5 |
---|
182 | |
---|
183 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
184 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
185 | Make/D/N=(nord) $weightStr,$zStr |
---|
186 | Wave wt = $weightStr |
---|
187 | Wave xi = $zStr // wave references to pass |
---|
188 | Make5GaussPoints(wt,xi) |
---|
189 | else |
---|
190 | if(exists(weightStr) > 1) |
---|
191 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
192 | endif |
---|
193 | Wave wt = $weightStr |
---|
194 | Wave xi = $zStr // create the wave references |
---|
195 | endif |
---|
196 | |
---|
197 | Variable ii,jj,kk,ax,bx,ay,by,num |
---|
198 | Variable qx,qy,qz,qval,sqx,sqy,fs |
---|
199 | Variable qy_pt,qx_pt,res_x,res_y,res_tot,answer,sumIn,sumOut |
---|
200 | num=numpnts(s.qx) |
---|
201 | |
---|
202 | // end points of integration |
---|
203 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
204 | // +/- 3 sigq catches 99.73% of distrubution |
---|
205 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
206 | //integration from va to vb |
---|
207 | Make/O/D/N=1 fcnRet,xptW,yPtw |
---|
208 | |
---|
209 | answer=0 |
---|
210 | //loop over q-values |
---|
211 | for(ii=0;ii<num;ii+=1) |
---|
212 | qx = s.qx[ii] |
---|
213 | qy = s.qy[ii] |
---|
214 | qz = s.qz[ii] |
---|
215 | qval = sqrt(qx^2+qy^2+qz^2) |
---|
216 | sqx = s.sigQx[ii] |
---|
217 | sqy = s.sigQy[ii] |
---|
218 | fs = s.fs[ii] |
---|
219 | |
---|
220 | ax = -3*sqx + qx //qx integration limits |
---|
221 | bx = 3*sqx + qx |
---|
222 | ay = -3*sqy + qy //qy integration limits |
---|
223 | by = 3*sqy + qy |
---|
224 | |
---|
225 | // 5-pt quadrature loops |
---|
226 | sumOut = 0 |
---|
227 | for(jj=0;jj<nord;jj+=1) // call qy the "outer' |
---|
228 | qy_pt = (xi[jj]*(by-ay)+ay+by)/2 |
---|
229 | res_y = exp((-1*(qy - qy_pt)^2)/(2*sqy*sqy)) |
---|
230 | |
---|
231 | sumIn=0 |
---|
232 | for(kk=0;kk<nord;kk+=1) |
---|
233 | |
---|
234 | qx_pt = (xi[kk]*(bx-ax)+ax+bx)/2 |
---|
235 | res_x = exp((-1*(qx - qx_pt)^2)/(2*sqx*sqx)) |
---|
236 | |
---|
237 | res_tot = res_x*res_y/(2*pi*sqx*sqy) |
---|
238 | xptw[0] = qx_pt |
---|
239 | yptw[0] = qy_pt |
---|
240 | fcn(s.coefW,fcnRet,xptw,yptw) //fcn passed in is an AAO |
---|
241 | sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] |
---|
242 | endfor |
---|
243 | answer += (bx-ax)/2.0*sumIn //this is NOT the right normalization |
---|
244 | endfor |
---|
245 | |
---|
246 | answer *= (by-ay)/2.0 |
---|
247 | s.zw[ii] = answer |
---|
248 | // s.zw[ii] = sumIn |
---|
249 | endfor |
---|
250 | |
---|
251 | |
---|
252 | return(0) |
---|
253 | end |
---|
254 | |
---|
255 | Function Smear_2DModel_20(fcn,s) |
---|
256 | FUNCREF SANS_2D_ModelAAO_proto fcn |
---|
257 | Struct ResSmear_2D_AAOStruct &s |
---|
258 | |
---|
259 | String weightStr="gauss20wt",zStr="gauss20z" |
---|
260 | Variable nord=20 |
---|
261 | |
---|
262 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
263 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
264 | Make/D/N=(nord) $weightStr,$zStr |
---|
265 | Wave wt = $weightStr |
---|
266 | Wave xi = $zStr // wave references to pass |
---|
267 | Make20GaussPoints(wt,xi) |
---|
268 | else |
---|
269 | if(exists(weightStr) > 1) |
---|
270 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
271 | endif |
---|
272 | Wave wt = $weightStr |
---|
273 | Wave xi = $zStr // create the wave references |
---|
274 | endif |
---|
275 | |
---|
276 | Variable ii,jj,kk,ax,bx,ay,by,num |
---|
277 | Variable qx,qy,qz,qval,sqx,sqy,fs |
---|
278 | Variable qy_pt,qx_pt,res_x,res_y,res_tot,answer,sumIn,sumOut |
---|
279 | num=numpnts(s.qx) |
---|
280 | |
---|
281 | // end points of integration |
---|
282 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
283 | // +/- 3 sigq catches 99.73% of distrubution |
---|
284 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
285 | //integration from va to vb |
---|
286 | Make/O/D/N=1 fcnRet,xptW,yPtw |
---|
287 | |
---|
288 | answer=0 |
---|
289 | //loop over q-values |
---|
290 | for(ii=0;ii<num;ii+=1) |
---|
291 | qx = s.qx[ii] |
---|
292 | qy = s.qy[ii] |
---|
293 | qz = s.qz[ii] |
---|
294 | qval = sqrt(qx^2+qy^2+qz^2) |
---|
295 | sqx = s.sigQx[ii] |
---|
296 | sqy = s.sigQy[ii] |
---|
297 | fs = s.fs[ii] |
---|
298 | |
---|
299 | ax = -3*sqx + qx //qx integration limits |
---|
300 | bx = 3*sqx + qx |
---|
301 | ay = -3*sqy + qy //qy integration limits |
---|
302 | by = 3*sqy + qy |
---|
303 | |
---|
304 | // 20-pt quadrature loops |
---|
305 | sumOut = 0 |
---|
306 | for(jj=0;jj<nord;jj+=1) // call qy the "outer' |
---|
307 | qy_pt = (xi[jj]*(by-ay)+ay+by)/2 |
---|
308 | res_y = exp((-1*(qy - qy_pt)^2)/(2*sqy*sqy)) |
---|
309 | |
---|
310 | sumIn=0 |
---|
311 | for(kk=0;kk<nord;kk+=1) |
---|
312 | |
---|
313 | qx_pt = (xi[kk]*(bx-ax)+ax+bx)/2 |
---|
314 | res_x = exp((-1*(qx - qx_pt)^2)/(2*sqx*sqx)) |
---|
315 | |
---|
316 | res_tot = res_x*res_y/(2*pi*sqx*sqy) |
---|
317 | xptw[0] = qx_pt |
---|
318 | yptw[0] = qy_pt |
---|
319 | fcn(s.coefW,fcnRet,xptw,yptw) //fcn passed in is an AAO |
---|
320 | sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] |
---|
321 | endfor |
---|
322 | answer += (bx-ax)/2.0*sumIn //this is NOT the right normalization |
---|
323 | endfor |
---|
324 | |
---|
325 | answer *= (by-ay)/2.0 |
---|
326 | s.zw[ii] = answer |
---|
327 | // s.zw[ii] = sumIn |
---|
328 | endfor |
---|
329 | |
---|
330 | |
---|
331 | return(0) |
---|
332 | end |
---|