source: sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf @ 601

Last change on this file since 601 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.
File size: 8.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4Function 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       
30end
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//
39Function 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
124END
125
126/// --- either way, same to machine precision
127Function 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)
146END
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//
176Function 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)
253end
254
255Function 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)
332end
Note: See TracBrowser for help on using the repository browser.