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

Last change on this file since 524 was 524, checked in by srkline, 14 years ago

first pass at utilities for 2D smearing. they're all wrong now, and need more work before they are anywhere near useful.

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                Make5GaussPoints(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                // 5-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.