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

Last change on this file since 801 was 801, checked in by srkline, 11 years ago

Smear_2D procedures to plot resolution function at each pixel have been improved, and to integrate the resolution to check that it normalizes to 1.

+ other minor tweaks to convert VAX time and to clean up list behavior.

File size: 16.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4
5
6// there are utlity functions here (commented out) for plotting the 2D resolution function
7// at any point on the detector
8
9
10// I can only calculate (nord) points in an AAO fashion, so threading the function calculation
11// doesn't help. Even if I can rewrite this to calculate nord*nord AAO, that will typically be 10*10=100
12// and that is not enough points to thread with any benefit
13
14// but the calculation of each q value is independent, so I can split the total number of q-values between processors
15//
16// -- but I must be careful to pass this a function that is not already threaded!
17// --- BUT - I can't pass either a function reference, OR a structure to a thread!
18// ---- So sadly, each 2D resolution calculation must be threaded by hand.
19//
20// See the Sphere_2D function for an example of how to thread the smearing
21//
22//
23// SRK May 2010
24
25//
26// NOTE: there is a definition of FindTheta = -1* FindPhi() that is a duplicate of what is in RawWindowHook.ipf
27// and should eventually be in a common location for analysis and reduction packages
28//
29
30
31//// this is the completely generic 2D smearing, not threaded, but takes FUNCREF and STRUCT parameters
32//
33// uses resolution ellipse defined perpendicular "Y" and parallel "X",
34// rotating the ellipse into its proper orientaiton based on qx,qy
35//
36// 5 gauss points is not enough - it gives artifacts as a function of phi
37// 10 gauss points is minimally sufficient
38// 20 gauss points are needed if lots of oscillations (just like in 1D)
39// even more may be necessary for highly peaked functions
40//
41//
42Function Smear_2DModel_PP(fcn,s,nord)
43        FUNCREF SANS_2D_ModelAAO_proto fcn
44        Struct ResSmear_2D_AAOStruct &s
45        Variable nord
46       
47        String weightStr,zStr
48
49        Variable ii,jj,kk,num
50        Variable qx,qy,qz,qval,fs
51        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
52       
53        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3
54       
55        switch(nord)   
56                case 5:         
57                        weightStr="gauss5wt"
58                        zStr="gauss5z"
59                        if (WaveExists($weightStr) == 0)
60                                Make/O/D/N=(nord) $weightStr,$zStr
61                                Make5GaussPoints($weightStr,$zStr)     
62                        endif
63                        break                           
64                case 10:               
65                        weightStr="gauss10wt"
66                        zStr="gauss10z"
67                        if (WaveExists($weightStr) == 0)
68                                Make/O/D/N=(nord) $weightStr,$zStr
69                                Make10GaussPoints($weightStr,$zStr)     
70                        endif
71                        break                           
72                case 20:               
73                        weightStr="gauss20wt"
74                        zStr="gauss20z"
75                        if (WaveExists($weightStr) == 0)
76                                Make/O/D/N=(nord) $weightStr,$zStr
77                                Make20GaussPoints($weightStr,$zStr)     
78                        endif
79                        break
80                default:                                                       
81                        Abort "Smear_2DModel_PP called with invalid nord value"                                 
82        endswitch
83       
84        Wave/Z wt = $weightStr
85        Wave/Z xi = $zStr
86       
87/// keep these waves local
88//      Make/O/D/N=1 yPtW
89        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
90       
91        // now just loop over the points as specified
92        num=numpnts(s.xw[0])
93       
94        answer=0
95       
96        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
97        Variable qperp_pt,phi_prime,q_prime
98
99        Variable t1=StopMSTimer(-2)
100
101        //loop over q-values
102        for(ii=0;ii<num;ii+=1)
103       
104//              if(mod(ii, 1000 ) == 0)
105//                      Print "ii= ",ii
106//              endif
107               
108                qx = s.xw[0][ii]
109                qy = s.xw[1][ii]
110                qz = s.qz[ii]
111                qval = sqrt(qx^2+qy^2+qz^2)
112                spl = s.sQpl[ii]
113                spp = s.sQpp[ii]
114                fs = s.fs[ii]
115               
116                normFactor = 2*pi*spl*spp
117               
118                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1
119               
120                apl = -numStdDev*spl + qval             //parallel = q integration limits
121                bpl = numStdDev*spl + qval
122///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG)
123///             bpp = numStdDev*spp + phi
124                app = -numStdDev*spp + 0                //q_perp = 0
125                bpp = numStdDev*spp + 0
126               
127                //make sure the limits are reasonable.
128                if(apl < 0)
129                        apl = 0
130                endif
131                // do I need to specially handle limits when phi ~ 0?
132       
133               
134                sumOut = 0
135                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
136///                     phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
137                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
138
139                        sumIn=0
140                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
141
142                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
143                               
144///                             FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
145
146                                // find QxQy given Qpl and Qperp on the grid
147                                //
148                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
149                                phi_prime = phi + qperp_pt/qpl_pt
150                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
151                               
152                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
153                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
154
155                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
156///                             res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
157                                res_tot[kk] /= normFactor
158//                              res_tot[kk] *= fs
159
160                        endfor
161                       
162                        fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord pts at a time
163                       
164                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
165                        fcnRet *= wt[jj]*wt*res_tot
166                        //
167                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
168                endfor
169
170                answer *= (bpp-app)/2.0
171                s.zw[ii] = answer
172        endfor
173       
174        Variable elap = (StopMSTimer(-2) - t1)/1e6
175        Print "elapsed time = ",elap
176       
177        return(0)
178end
179
180
181// this is generic, but I need to declare the Cylinder2D function threadsafe
182// and this calculation is significantly slower than the manually threaded calculation
183// if the function is fast to calculate. Once the function has polydispersity on 2 or more parameters
184// then this AAO calculation and the manual threading are both painfully slow, and more similar in execution time
185//
186// For 128x128 data, and 10x10 smearing, using 25 points for each polydispersity (using DANSE code)
187//      Successively making more things polydisperse gives the following timing (in seconds):
188//
189//                              monodisp                sigR            sigTheta                sigPhi
190// manual THR           1.1             3.9                     74                      1844
191//      this AAO                        8.6             11.4                    104             1930
192//
193// and using 5 points for each polydispersity: (I see no visual difference)
194// manual THR           1.1             1.6                     3.9             16.1
195//
196// so clearly -- use 5 points for the polydispersities, unless there's a good reason not to - and
197// certainly for a survey, it's the way to go.
198//
199Function Smear_2DModel_PP_AAO(fcn,s,nord)
200        FUNCREF SANS_2D_ModelAAO_proto fcn
201        Struct ResSmear_2D_AAOStruct &s
202        Variable nord
203       
204        String weightStr,zStr
205
206        Variable ii,jj,kk,num
207        Variable qx,qy,qz,qval,fs
208        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
209       
210        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3
211       
212        switch(nord)   
213                case 5:         
214                        weightStr="gauss5wt"
215                        zStr="gauss5z"
216                        if (WaveExists($weightStr) == 0)
217                                Make/O/D/N=(nord) $weightStr,$zStr
218                                Make5GaussPoints($weightStr,$zStr)     
219                        endif
220                        break                           
221                case 10:               
222                        weightStr="gauss10wt"
223                        zStr="gauss10z"
224                        if (WaveExists($weightStr) == 0)
225                                Make/O/D/N=(nord) $weightStr,$zStr
226                                Make10GaussPoints($weightStr,$zStr)     
227                        endif
228                        break                           
229                case 20:               
230                        weightStr="gauss20wt"
231                        zStr="gauss20z"
232                        if (WaveExists($weightStr) == 0)
233                                Make/O/D/N=(nord) $weightStr,$zStr
234                                Make20GaussPoints($weightStr,$zStr)     
235                        endif
236                        break
237                default:                                                       
238                        Abort "Smear_2DModel_PP called with invalid nord value"                                 
239        endswitch
240       
241        Wave/Z wt = $weightStr
242        Wave/Z xi = $zStr
243       
244/// keep these waves local
245//      Make/O/D/N=1 yPtW
246        Make/O/D/N=(nord*nord) fcnRet,xptW,res_tot,yptW,wts
247        Make/O/D/N=(nord) phi_pt,qpl_pt,qperp_pt
248               
249        // now just loop over the points as specified
250        num=numpnts(s.xw[0])
251       
252        answer=0
253       
254        Variable spl,spp,apl,app,bpl,bpp
255        Variable phi_prime,q_prime
256
257        Variable t1=StopMSTimer(-2)
258
259        //loop over q-values
260        for(ii=0;ii<num;ii+=1)
261       
262//              if(mod(ii, 1000 ) == 0)
263//                      Print "ii= ",ii
264//              endif
265               
266                qx = s.xw[0][ii]
267                qy = s.xw[1][ii]
268                qz = s.qz[ii]
269                qval = sqrt(qx^2+qy^2+qz^2)
270                spl = s.sQpl[ii]
271                spp = s.sQpp[ii]
272                fs = s.fs[ii]
273               
274                normFactor = 2*pi*spl*spp
275               
276                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1
277               
278                apl = -numStdDev*spl + qval             //parallel = q integration limits
279                bpl = numStdDev*spl + qval
280///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG)
281///             bpp = numStdDev*spp + phi
282                app = -numStdDev*spp + 0                //q_perp = 0
283                bpp = numStdDev*spp + 0
284               
285                //make sure the limits are reasonable.
286                if(apl < 0)
287                        apl = 0
288                endif
289                // do I need to specially handle limits when phi ~ 0?
290       
291               
292//              sumOut = 0
293                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
294//                      phi_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2
295                        qperp_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2             //this is now q_perp
296                       
297//                      sumIn=0
298                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
299
300                                qpl_pt[kk] = (xi[kk]*(bpl-apl)+apl+bpl)/2
301                               
302///                             FindQxQy(qpl_pt[kk],phi_pt[jj],qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
303                               
304                                // find QxQy given Qpl and Qperp on the grid
305                                //
306                                q_prime = sqrt(qpl_pt[kk]^2+qperp_pt[jj]^2)
307                                phi_prime = phi + qperp_pt[jj]/qpl_pt[kk]
308                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
309                                                               
310                                yPtw[nord*jj+kk] = qy_pt                                        //phi is the same in this loop, but qy is not
311                                xPtW[nord*jj+kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
312                               
313                                res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (qperp_pt[jj])^2/spp/spp ) )
314//                              res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (phi_pt[jj]-phi)^2/spp/spp ) )
315                                res_tot[nord*jj+kk] /= normFactor
316//                              res_tot[kk] *= fs
317
318                                //weighting
319                                wts[nord*jj+kk] = wt[jj]*wt[kk]
320                        endfor
321                       
322                endfor
323               
324                fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord*nord pts at a time
325               
326                fcnRet *= wts*res_tot
327                //
328                answer = (bpl-apl)/2.0*sum(fcnRet)              // get the sum, normalize to parallel direction
329                answer *= (bpp-app)/2.0                                         // and normalize to perpendicular direction
330               
331                s.zw[ii] = answer
332        endfor
333       
334        Variable elap = (StopMSTimer(-2) - t1)/1e6
335        Print "elapsed time = ",elap
336       
337        return(0)
338end
339
340
341
342
343
344
345//phi is defined from +x axis, proceeding CCW around [0,2Pi]
346//rotate the resolution function by theta,  = -phi
347//
348// this is only different by (-1) from FindPhi
349// I'd just call FindPhi, but it's awkward to include
350//
351Threadsafe Function FindTheta(vx,vy)
352        variable vx,vy
353
354       
355        return(-1 * FindPhi(vx,vy))
356end
357
358
359
360       
361//////////////////////////////////////////////////////////
362//////////////////////////////////////////////////////////
363/// utility functions to test the calculation of the resolution function and
364//       plotting of it for visualization
365//
366// -- these need to have the reduction package included for it to compile
367//
368// this works with the raw 2D data, calculating resolution in place
369// type is "RAW"
370// xx,yy are in detector coordinates
371//
372// call as:
373// PlotResolution_atPixel("RAW",pcsr(A),qcsr(A))
374//
375// then:
376// Display;AppendImage res;AppendMatrixContour res;ModifyContour res labels=0,autoLevels={*,*,3}
377//
378
379Function PlotResolution_atPixel(type,xx,yy)
380        String type
381        Variable xx,yy
382       
383///from QxQyExport
384        String destStr="",typeStr=""
385        Variable step=1,refnum
386        destStr = "root:Packages:NIST:"+type
387       
388        //must select the linear_data to export
389        NVAR isLog = $(destStr+":gIsLogScale")
390        if(isLog==1)
391                typeStr = ":linear_data"
392        else
393                typeStr = ":data"
394        endif
395       
396        NVAR pixelsX = root:myGlobals:gNPixelsX
397        NVAR pixelsY = root:myGlobals:gNPixelsY
398       
399        Wave data=$(destStr+typeStr)
400        WAVE intw=$(destStr + ":integersRead")
401        WAVE rw=$(destStr + ":realsRead")
402        WAVE/T textw=$(destStr + ":textRead")
403       
404//      Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist
405        Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist
406
407        Variable xctr,yctr,sdd,lambda,pixSize
408        xctr = rw[16]
409        yctr = rw[17]
410        sdd = rw[18]
411        lambda = rw[26]
412        pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels)
413       
414//      qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system
415//      qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
416       
417        qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)                //+1 converts to detector coordinate system
418        qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
419       
420//      Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val
421
422///************
423// do everything to write out the resolution too
424        // un-comment these if you want to write out qz_val and qval too, then use the proper save command
425        qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
426        qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
427        phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr))           //(dx,dy)
428        r_dist = sqrt(  (pixSize*((xx+1)-xctr))^2 +  (pixSize*((yy+1)-yctr))^2 )                //radial distance from ctr to pt
429//      Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist
430        //everything in 1D now
431//      Duplicate/O qval SigmaQX,SigmaQY,fsubS
432        Variable SigmaQX,SigmaQY,fsubS
433
434        Variable L2 = rw[18]
435        Variable BS = rw[21]
436        Variable S1 = rw[23]
437        Variable S2 = rw[24]
438        Variable L1 = rw[25]
439        Variable lambdaWidth = rw[27]   
440        Variable usingLenses = rw[28]           //new 2007
441
442        //Two parameters DDET and APOFF are instrument dependent.  Determine
443        //these from the instrument name in the header.
444        //From conversation with JB on 01.06.99 these are the current good values
445        Variable DDet
446        NVAR apOff = root:myGlobals:apOff               //in cm
447        DDet = rw[10]/10                        // header value (X) is in mm, want cm here
448
449        Variable ret1,ret2,ret3,del_r
450        del_r = rw[10]
451       
452        get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3)
453        SigmaQX = ret1 
454        SigmaQY = ret2 
455        fsubs = ret3   
456/////
457
458//      Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0
459        Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor
460        Variable qx_ret,qy_ret
461       
462//      theta = FindPhi(qx_val,qy_val)
463// need to rotate properly - theta is defined a starting from +y axis, moving CW
464// we define phi starting from +x and moving CCW
465        theta = -phi                    //seems to give the right behavior...
466       
467       
468        Print qx_val,qy_val,qval
469        Print "phi, theta",phi,theta
470       
471        FindQxQy(qval,phi,qx_ret,qy_ret)
472       
473        sx = SigmaQx
474        sy = sigmaQy
475        x0 = qx_val
476        y0 = qy_val
477       
478        a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy)
479        b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy)
480        c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy)
481       
482        normFactor = pi/sqrt(a*c-b*b)
483
484        Make/O/D/N=(num,num) res
485        // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle
486        maxSig = max(sx,sy)
487        Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res
488        Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res
489//      Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res
490//      Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res
491       
492        Variable xPt,yPt,delx,dely,offx,offy
493        delx = DimDelta(res,0)
494        dely = DimDelta(res,1)
495        offx = DimOffset(res,0)
496        offy = DimOffset(res,1)
497
498        Print "sx,sy = ",sx,sy
499        for(ii=0;ii<num;ii+=1)
500                xPt = offx + ii*delx
501                for(jj=0;jj<num;jj+=1)
502                        yPt = offy + jj*dely
503                        res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2))
504                endfor
505        endfor 
506        res /= normFactor
507       
508        //Print sum(res,-inf,inf)*delx*dely
509        if(WaveExists($"coef")==0)
510                Make/O/D/N=6 coef
511        endif
512        Wave coef=coef
513        coef[0] = 1
514        coef[1] = qx_val
515        coef[2] = qy_val
516        coef[3] = sx
517        coef[4] = sy
518        coef[5] = theta
519
520//      Variable t1=StopMSTimer(-2)
521
522//
523        do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0)
524//
525
526//      Variable elap = (StopMSTimer(-2) - t1)/1e6
527//      Print "elapsed time = ",elap
528//      Print "time for 16384 = (minutes)",16384*elap/60
529        return(0)
530End
531
532
533// this is called each time to integrate the gaussian
534Function do2dIntegrationGauss(xMin,xMax,yMin,yMax)
535        Variable xMin,xMax,yMin,yMax
536       
537        Variable/G globalXmin=xMin
538        Variable/G globalXmax=xMax
539        Variable/G globalY
540                       
541        Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)       
542        KillVariables/z globalXmax,globalXmin,globalY
543        print "integration of 2D = ",result
544End
545
546Function Gauss2DFuncOuter(inY)
547        Variable inY
548       
549        NVAR globalXmin,globalXmax,globalY
550        globalY=inY
551       
552        return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)         
553End
554
555Function Gauss2DFuncInner(inX)
556        Variable inX
557       
558        NVAR globalY
559        Wave coef=coef
560       
561        return Gauss2D_theta(coef,inX,GlobalY)
562End
563
564Function Gauss2D_theta(w,x,y)
565        Wave w
566        Variable x,y
567       
568        Variable val,a,b,c
569        Variable scale,x0,y0,sx,sy,theta,normFactor
570       
571        scale = w[0]
572        x0 = w[1]
573        y0 = w[2]
574        sx = w[3]
575        sy = w[4]
576        theta = w[5]
577       
578        a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy)
579        b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy)
580        c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy)
581       
582        val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2))
583       
584        normFactor = pi/sqrt(a*c-b*b)
585       
586        return(scale*val/normFactor)
587end
588
589////////////////////////////////////////////////////////////
590////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.