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

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

Corrected number of quadrature points in 20-pt smearing

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