source: sans/Dev/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf @ 819

Last change on this file since 819 was 800, checked in by srkline, 12 years ago

Changed the 2D resolution calculation to include gravity. This seems to work, but will still require some testing. It requires no modifications to the smearing calculation, still using parallel and perpendicular directions.

Added a Fake Update feature to the RealTime? update. There are specific, separate instructions for how to set this up. Usef (possibly) for summer schools or demos.

Adjusted the 2D MonteCarlo? simulation to a corrected beam center. The Gauss Peak was not symmetric around the old beam center.

Corrected the AAO resolution smearing functions to work with cursors.
File size: 51.2 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=4.00
3#pragma IgorVersion=6.1
4
5// GaussUtils.ipf
6// 22dec97 SRK
7
8//// NEW
9//
10// added two utility functions to do numerical
11// integration of an arbitrary function
12// Gaussian quadrature of 20 or 76 points is used
13//
14// fcn must be defined as in the prototype function
15// which is similar to the normal fitting function definition
16//
17// 09 SEP 03 SRK
18//
19//
20// 04DEC03 - added routines for 5 and 10 Gauss points
21//
22// 13 DEC 04 BSG/SRK - Modified Smear_Model_(5)(20) to allow
23// smearing of USANS data. dQv is passed in as negative value
24// in all resolution columns. the dQv value is read from the SigQ columnn
25// (the 4th column) - and all values are read, so fill the whole column
26//
27// Adaptive trapezoidal integration is used, 76 Gauss pts
28// did not give sufficient accuracy.
29//
30////////////////////////////////////////////////
31// 23 JUL 2007 SRK
32// major overhaul to use AAO functions in the smearing calculations
33//
34// - re-definition of function prototype in Smear_Model_N
35// - re-definition in trapezoidal routines too
36// - calls from model functions are somewhat different, so this will generate a lot of errors
37//   until all can be changed
38// - now is looking for a resolution matrix that contains the 3 resolution waves plus Q
39// = mat[num_Qvals][0] = sigQ
40// = mat[num_Qvals][1] = Qbar
41// = mat[num_Qvals][2] = fShad
42// = mat[num_Qvals][3] = qvals
43//
44// -- does not yet use the matrix calculation for USANS
45// -- SO THERE IS NO SWITCH YET IN Smear_Model_N()
46//
47//
48
49//maybe not the optimal set of parameters for the STRUCT
50//
51// resW is /N=(np,4) if SANS data
52// resW is /N=(np,np) if USANS data
53//
54// info may be useful later as a "KEY=value;" string to carry additional information
55Structure ResSmearAAOStruct
56        Wave coefW
57        Wave yW
58        Wave xW
59        Wave resW
60        String info
61EndStructure
62
63
64//// tentative pass at 2D resolution smearing
65////
66//Structure ResSmear_2D_AAOStruct
67//      Wave coefW
68//      Wave zw                 //answer
69//      Wave qy                 // q-value
70//      Wave qx
71//      Wave qz
72//      Wave sQpl               //resolution parallel to Q
73//      Wave sQpp               //resolution perpendicular to Q
74//      Wave fs
75//      String info
76//EndStructure
77
78// reformat the structure ?? WM Fit compatible
79// -- 2D resolution smearing
80//
81Structure ResSmear_2D_AAOStruct
82        Wave coefW
83        Wave zw                 //answer
84        Wave xw[2]              // qx-value is [0], qy is xw[1]
85        Wave qz
86        Wave sQpl               //resolution parallel to Q
87        Wave sQpp               //resolution perpendicular to Q
88        Wave fs
89        String info
90EndStructure
91
92
93
94Function Make5GaussPoints(w5,z5)
95        Wave w5,z5
96
97//      printf  "in make Gauss Pts\r"
98
99        z5[0] = -.906179845938664
100        z5[1] = -.538469310105683
101        z5[2] = -.0000000000000
102        z5[3] = .538469310105683
103        z5[4] = .906179845938664
104
105        w5[0] = .236926885056189
106        w5[1] = .478628670499366
107        w5[2] = .56888888888889
108        w5[3] = .478628670499366
109        w5[4] = .236926885056189
110
111//          printf "w[0],z[0] = %g %g\r", w5[0],z5[0]
112End
113
114Function Make10GaussPoints(w10,z10)
115        Wave w10,z10
116
117//      printf  "in make Gauss Pts\r"
118        z10[0] = -.973906528517172
119        z10[1] = -.865063366688985
120        z10[2] = -.679409568299024
121        z10[3] = -.433395394129247
122        z10[4] = -.148874338981631
123        z10[5] = .148874338981631
124        z10[6] = .433395394129247
125        z10[7] = .679409568299024
126        z10[8] = .865063366688985
127        z10[9] = .973906528517172
128       
129        w10[0] = .066671344308688
130        w10[1] = 0.149451349150581
131        w10[2] = 0.219086362515982
132        w10[3] = .269266719309996
133        w10[4] = 0.295524224714753
134        w10[5] = 0.295524224714753
135        w10[6] = .269266719309996
136        w10[7] = 0.219086362515982
137        w10[8] = 0.149451349150581
138        w10[9] = .066671344308688
139       
140//          printf "w[0],z[0] = %g %g\r", w10[0],z10[0]
141End
142
143Function Make20GaussPoints(w20,z20)
144        Wave w20,z20
145
146//      printf  "in make Gauss Pts\r"
147
148         z20[0] = -.993128599185095
149         z20[1] =  -.963971927277914
150         z20[2] =    -.912234428251326
151         z20[3] =    -.839116971822219
152         z20[4] =   -.746331906460151
153         z20[5] =   -.636053680726515
154         z20[6] =   -.510867001950827
155         z20[7] =     -.37370608871542
156         z20[8] =     -.227785851141645
157         z20[9] =     -.076526521133497
158         z20[10] =     .0765265211334973
159         z20[11] =     .227785851141645
160         z20[12] =     .37370608871542
161         z20[13] =     .510867001950827
162         z20[14] =     .636053680726515
163         z20[15] =     .746331906460151
164         z20[16] =     .839116971822219
165         z20[17] =     .912234428251326
166         z20[18] =    .963971927277914
167         z20[19] =  .993128599185095
168           
169        w20[0] =  .0176140071391521
170        w20[1] =  .0406014298003869
171        w20[2] =      .0626720483341091
172        w20[3] =      .0832767415767047
173        w20[4] =     .10193011981724
174        w20[5] =      .118194531961518
175        w20[6] =      .131688638449177
176        w20[7] =      .142096109318382
177        w20[8] =      .149172986472604
178        w20[9] =      .152753387130726
179        w20[10] =      .152753387130726
180        w20[11] =      .149172986472604
181        w20[12] =      .142096109318382
182        w20[13] =     .131688638449177
183        w20[14] =     .118194531961518
184        w20[15] =     .10193011981724
185        w20[16] =     .0832767415767047
186        w20[17] =     .0626720483341091
187        w20[18] =     .0406014298003869
188        w20[19] =     .0176140071391521
189//          printf "w[0],z[0] = %g %g\r", w20[0],z20[0]
190End
191
192Function Make76GaussPoints(w76,z76)
193        Wave w76,z76
194
195//      printf  "in make Gauss Pts\r"
196       
197                z76[0] = .999505948362153*(-1.0)
198            z76[75] = -.999505948362153*(-1.0)
199            z76[1] = .997397786355355*(-1.0)
200            z76[74] = -.997397786355355*(-1.0)
201            z76[2] = .993608772723527*(-1.0)
202            z76[73] = -.993608772723527*(-1.0)
203            z76[3] = .988144453359837*(-1.0)
204            z76[72] = -.988144453359837*(-1.0)
205            z76[4] = .981013938975656*(-1.0)
206            z76[71] = -.981013938975656*(-1.0)
207            z76[5] = .972229228520377*(-1.0)
208            z76[70] = -.972229228520377*(-1.0)
209            z76[6] = .961805126758768*(-1.0)
210            z76[69] = -.961805126758768*(-1.0)
211            z76[7] = .949759207710896*(-1.0)
212            z76[68] = -.949759207710896*(-1.0)
213            z76[8] = .936111781934811*(-1.0)
214            z76[67] = -.936111781934811*(-1.0)
215            z76[9] = .92088586125215*(-1.0)
216            z76[66] = -.92088586125215*(-1.0)
217            z76[10] = .904107119545567*(-1.0)
218            z76[65] = -.904107119545567*(-1.0)
219            z76[11] = .885803849292083*(-1.0)
220            z76[64] = -.885803849292083*(-1.0)
221            z76[12] = .866006913771982*(-1.0)
222            z76[63] = -.866006913771982*(-1.0)
223            z76[13] = .844749694983342*(-1.0)
224            z76[62] = -.844749694983342*(-1.0)
225            z76[14] = .822068037328975*(-1.0)
226            z76[61] = -.822068037328975*(-1.0)
227            z76[15] = .7980001871612*(-1.0)
228            z76[60] = -.7980001871612*(-1.0)
229            z76[16] = .77258672828181*(-1.0)
230            z76[59] = -.77258672828181*(-1.0)
231            z76[17] = .74587051350361*(-1.0)
232            z76[58] = -.74587051350361*(-1.0)
233            z76[18] = .717896592387704*(-1.0)
234            z76[57] = -.717896592387704*(-1.0)
235            z76[19] = .688712135277641*(-1.0)
236            z76[56] = -.688712135277641*(-1.0)
237            z76[20] = .658366353758143*(-1.0)
238            z76[55] = -.658366353758143*(-1.0)
239            z76[21] = .626910417672267*(-1.0)
240            z76[54] = -.626910417672267*(-1.0)
241            z76[22] = .594397368836793*(-1.0)
242            z76[53] = -.594397368836793*(-1.0)
243            z76[23] = .560882031601237*(-1.0)
244            z76[52] = -.560882031601237*(-1.0)
245            z76[24] = .526420920401243*(-1.0)
246            z76[51] = -.526420920401243*(-1.0)
247            z76[25] = .491072144462194*(-1.0)
248            z76[50] = -.491072144462194*(-1.0)
249            z76[26] = .454895309813726*(-1.0)
250            z76[49] = -.454895309813726*(-1.0)
251            z76[27] = .417951418780327*(-1.0)
252            z76[48] = -.417951418780327*(-1.0)
253            z76[28] = .380302767117504*(-1.0)
254            z76[47] = -.380302767117504*(-1.0)
255            z76[29] = .342012838966962*(-1.0)
256            z76[46] = -.342012838966962*(-1.0)
257            z76[30] = .303146199807908*(-1.0)
258            z76[45] = -.303146199807908*(-1.0)
259            z76[31] = .263768387584994*(-1.0)
260            z76[44] = -.263768387584994*(-1.0)
261            z76[32] = .223945802196474*(-1.0)
262            z76[43] = -.223945802196474*(-1.0)
263            z76[33] = .183745593528914*(-1.0)
264            z76[42] = -.183745593528914*(-1.0)
265            z76[34] = .143235548227268*(-1.0)
266            z76[41] = -.143235548227268*(-1.0)
267            z76[35] = .102483975391227*(-1.0)
268            z76[40] = -.102483975391227*(-1.0)
269            z76[36] = .0615595913906112*(-1.0)
270            z76[39] = -.0615595913906112*(-1.0)
271            z76[37] = .0205314039939986*(-1.0)
272            z76[38] = -.0205314039939986*(-1.0)
273           
274                w76[0] =  .00126779163408536
275                w76[75] = .00126779163408536
276                w76[1] =  .00294910295364247
277            w76[74] = .00294910295364247
278            w76[2] = .00462793522803742
279            w76[73] =  .00462793522803742
280            w76[3] = .00629918049732845
281            w76[72] = .00629918049732845
282            w76[4] = .00795984747723973
283            w76[71] = .00795984747723973
284            w76[5] = .00960710541471375
285            w76[70] =  .00960710541471375
286            w76[6] = .0112381685696677
287            w76[69] = .0112381685696677
288            w76[7] =  .0128502838475101
289            w76[68] = .0128502838475101
290            w76[8] = .0144407317482767
291            w76[67] =  .0144407317482767
292            w76[9] = .0160068299122486
293            w76[66] = .0160068299122486
294            w76[10] = .0175459372914742
295            w76[65] = .0175459372914742
296            w76[11] = .0190554584671906
297            w76[64] = .0190554584671906
298            w76[12] = .020532847967908
299            w76[63] = .020532847967908
300            w76[13] = .0219756145344162
301            w76[62] = .0219756145344162
302            w76[14] = .0233813253070112
303            w76[61] = .0233813253070112
304            w76[15] = .0247476099206597
305            w76[60] = .0247476099206597
306            w76[16] = .026072164497986
307            w76[59] = .026072164497986
308            w76[17] = .0273527555318275
309            w76[58] = .0273527555318275
310            w76[18] = .028587223650054
311            w76[57] = .028587223650054
312            w76[19] = .029773487255905
313            w76[56] = .029773487255905
314            w76[20] = .0309095460374916
315            w76[55] = .0309095460374916
316            w76[21] = .0319934843404216
317            w76[54] = .0319934843404216
318            w76[22] = .0330234743977917
319            w76[53] = .0330234743977917
320            w76[23] = .0339977794120564
321            w76[52] = .0339977794120564
322            w76[24] = .0349147564835508
323            w76[51] = .0349147564835508
324            w76[25] = .0357728593807139
325            w76[50] = .0357728593807139
326            w76[26] = .0365706411473296
327            w76[49] = .0365706411473296
328            w76[27] = .0373067565423816
329            w76[48] = .0373067565423816
330            w76[28] = .0379799643084053
331            w76[47] = .0379799643084053
332            w76[29] = .0385891292645067
333            w76[46] = .0385891292645067
334            w76[30] = .0391332242205184
335            w76[45] = .0391332242205184
336            w76[31] = .0396113317090621
337            w76[44] = .0396113317090621
338            w76[32] = .0400226455325968
339            w76[43] = .0400226455325968
340            w76[33] = .040366472122844
341            w76[42] = .040366472122844
342            w76[34] = .0406422317102947
343            w76[41] = .0406422317102947
344            w76[35] = .0408494593018285
345            w76[40] = .0408494593018285
346            w76[36] = .040987805464794
347            w76[39] = .040987805464794
348            w76[37] = .0410570369162294
349            w76[38] = .0410570369162294
350
351End             //Make76GaussPoints()
352
353// !!!!! reduces the length of wt and zi by one !!!!!
354//
355Function Make_N_GaussPoints(wt,zi)
356        Wave wt,zi
357       
358        Variable num
359        num = numpnts(wt) - 1
360       
361        gauleg(-1,1,zi,wt,num)
362       
363        DeletePoints 0,1,wt,zi
364       
365        return(0)
366End
367
368/// gauleg subroutine from NR to calculate weights and abscissae for
369// Gauss-Legendre quadrature
370//
371//
372// arrays are indexed from 1
373//
374Function gauleg( x1, x2, x, w, n)
375        Variable x1, x2
376        Wave x, w
377        Variable n
378       
379        variable m,j,i
380        variable z1,z,xm,xl,pp,p3,p2,p1
381        Variable eps = 3e-11
382
383        m=(n+1)/2
384        xm=0.5*(x2+x1)
385        xl=0.5*(x2-x1)
386        for (i=1;i<=m;i+=1)
387                z=cos(pi*(i-0.25)/(n+0.5))
388                do
389                        p1=1.0
390                        p2=0.0
391                        for (j=1;j<=n;j+=1)
392                                p3=p2
393                                p2=p1
394                                p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
395                        endfor
396                        pp=n*(z*p1-p2)/(z*z-1.0)
397                        z1=z
398                        z=z1-p1/pp
399                while (abs(z-z1) > EPS)
400                x[i]=xm-xl*z
401                x[n+1-i]=xm+xl*z
402                w[i]=2.0*xl/((1.0-z*z)*pp*pp)
403                w[n+1-i]=w[i]
404        Endfor
405End
406
407
408/// uses a user-supplied number of Gauss points, and generates them on-the-fly as needed
409// using a Numerical Recipes routine
410//
411// - note that this has an extra input parameter, nord
412//
413////////////
414Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord)                               
415        FUNCREF GenericQuadrature_proto fcn
416        Variable loLim,upLim    //limits of integration
417        Wave w                  //coefficients of function fcn(w,x)
418        Variable x      //x-value (q) for the calculation
419        Variable nord                   //number of quadrature points to used
420
421
422// special case of integral limits that are identical
423        if(upLim == loLim)
424                return( fcn(w,x, loLim))
425        endif
426       
427// local variables
428        Variable ii,va,vb,summ,yyy,zi
429        Variable answer,dum
430        String weightStr,zStr
431       
432        weightStr = "gauss"+num2iStr(nord)+"wt"
433        zStr = "gauss"+num2istr(nord)+"z"
434               
435        if (WaveExists($weightStr) == 0) // wave reference is not valid,
436                Make/D/N=(nord+1) $weightStr,$zStr
437                Wave wt = $weightStr
438                Wave xx = $zStr         // wave references to pass
439                Make_N_GaussPoints(wt,xx)                               //generates the gauss points and removes the extra point
440        else
441                if(exists(weightStr) > 1)
442                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
443                endif
444                Wave wt = $weightStr
445                Wave xx = $zStr         // create the wave references
446        endif
447
448        //limits of integration are input to function
449        va = loLim
450        vb = upLim
451        // Using 5 Gauss points             
452        // remember to index from 0,size-1
453       
454        summ = 0.0              // initialize integral
455        ii=0                    // loop counter
456        do
457                // calculate Gauss points on integration interval (q-value for evaluation)
458                zi = ( xx[ii]*(vb-va) + vb + va )/2.0
459                //calculate partial sum for the passed-in model function       
460                yyy = wt[ii] *  fcn(w,x,zi)                                             
461                summ += yyy             //add to the running total of the quadrature
462       ii+=1           
463        while (ii<nord)                         // end of loop over quadrature points
464   
465        // calculate value of integral to return
466        answer = (vb-va)/2.0*summ
467
468        Return (answer)
469End
470
471
472
473////////////
474Function IntegrateFn5(fcn,loLim,upLim,w,x)                             
475        FUNCREF GenericQuadrature_proto fcn
476        Variable loLim,upLim    //limits of integration
477        Wave w                  //coefficients of function fcn(w,x)
478        Variable x      //x-value (q) for the calculation
479
480
481// special case of integral limits that are identical
482        if(upLim == loLim)
483                return( fcn(w,x, loLim))
484        endif
485       
486// local variables
487        Variable nord,ii,va,vb,summ,yyy,zi
488        Variable answer,dum
489        String weightStr,zStr
490       
491        weightStr = "gauss5wt"
492        zStr = "gauss5z"
493        nord = 5
494               
495        if (WaveExists($weightStr) == 0) // wave reference is not valid,
496                Make/D/N=5 $weightStr,$zStr
497                Wave w5 = $weightStr
498                Wave z5 = $zStr         // wave references to pass
499                Make5GaussPoints(w5,z5)
500        else
501                if(exists(weightStr) > 1)
502                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
503                endif
504                Wave w5 = $weightStr
505                Wave z5 = $zStr         // create the wave references
506        endif
507
508        //limits of integration are input to function
509        va = loLim
510        vb = upLim
511        // Using 5 Gauss points             
512        // remember to index from 0,size-1
513
514        summ = 0.0              // initialize integral
515        ii=0                    // loop counter
516        do
517                // calculate Gauss points on integration interval (q-value for evaluation)
518                zi = ( z5[ii]*(vb-va) + vb + va )/2.0
519                //calculate partial sum for the passed-in model function       
520                yyy = w5[ii] *  fcn(w,x,zi)                                             
521                summ += yyy             //add to the running total of the quadrature
522        ii+=1           
523        while (ii<nord)                         // end of loop over quadrature points
524   
525        // calculate value of integral to return
526        answer = (vb-va)/2.0*summ
527
528        Return (answer)
529End
530///////////////////////////////////////////////////////////////
531
532////////////
533Function IntegrateFn10(fcn,loLim,upLim,w,x)                             
534        FUNCREF GenericQuadrature_proto fcn
535        Variable loLim,upLim    //limits of integration
536        Wave w                  //coefficients of function fcn(w,x)
537        Variable x      //x-value (q) for the calculation
538
539// special case of integral limits that are identical
540        if(upLim == loLim)
541                return( fcn(w,x, loLim))
542        endif
543       
544// local variables
545        Variable nord,ii,va,vb,summ,yyy,zi
546        Variable answer,dum
547        String weightStr,zStr
548       
549        weightStr = "gauss10wt"
550        zStr = "gauss10z"
551        nord = 10
552               
553        if (WaveExists($weightStr) == 0) // wave reference is not valid,
554                Make/D/N=10 $weightStr,$zStr
555                Wave w10 = $weightStr
556                Wave z10 = $zStr                // wave references to pass
557                Make10GaussPoints(w10,z10)     
558        else
559                if(exists(weightStr) > 1)
560                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
561                endif
562                Wave w10 = $weightStr
563                Wave z10 = $zStr                // create the wave references
564        endif
565
566        //limits of integration are input to function
567        va = loLim
568        vb = upLim
569        // Using 10 Gauss points                   
570        // remember to index from 0,size-1
571
572        summ = 0.0              // initialize integral
573        ii=0                    // loop counter
574        do
575                // calculate Gauss points on integration interval (q-value for evaluation)
576                zi = ( z10[ii]*(vb-va) + vb + va )/2.0
577                //calculate partial sum for the passed-in model function       
578                yyy = w10[ii] *  fcn(w,x,zi)                                           
579                summ += yyy             //add to the running total of the quadrature
580        ii+=1           
581        while (ii<nord)                         // end of loop over quadrature points
582   
583        // calculate value of integral to return
584        answer = (vb-va)/2.0*summ
585
586        Return (answer)
587End
588///////////////////////////////////////////////////////////////
589
590////////////
591Function IntegrateFn20(fcn,loLim,upLim,w,x)                             
592        FUNCREF GenericQuadrature_proto fcn
593        Variable loLim,upLim    //limits of integration
594        Wave w                  //coefficients of function fcn(w,x)
595        Variable x      //x-value (q) for the calculation
596
597// special case of integral limits that are identical
598        if(upLim == loLim)
599                return( fcn(w,x, loLim))
600        endif
601       
602// local variables
603        Variable nord,ii,va,vb,summ,yyy,zi
604        Variable answer,dum
605        String weightStr,zStr
606       
607        weightStr = "gauss20wt"
608        zStr = "gauss20z"
609        nord = 20
610               
611        if (WaveExists($weightStr) == 0) // wave reference is not valid,
612                Make/D/N=20 $weightStr,$zStr
613                Wave w20 = $weightStr
614                Wave z20 = $zStr                // wave references to pass
615                Make20GaussPoints(w20,z20)     
616        else
617                if(exists(weightStr) > 1)
618                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
619                endif
620                Wave w20 = $weightStr
621                Wave z20 = $zStr                // create the wave references
622        endif
623
624        //limits of integration are input to function
625        va = loLim
626        vb = upLim
627        // Using 20 Gauss points                   
628        // remember to index from 0,size-1
629
630        summ = 0.0              // initialize integral
631        ii=0                    // loop counter
632        do
633                // calculate Gauss points on integration interval (q-value for evaluation)
634                zi = ( z20[ii]*(vb-va) + vb + va )/2.0
635                //calculate partial sum for the passed-in model function       
636                yyy = w20[ii] *  fcn(w,x,zi)                                           
637                summ += yyy             //add to the running total of the quadrature
638        ii+=1           
639        while (ii<nord)                         // end of loop over quadrature points
640   
641        // calculate value of integral to return
642        answer = (vb-va)/2.0*summ
643
644        Return (answer)
645End
646///////////////////////////////////////////////////////////////
647
648Function IntegrateFn76(fcn,loLim,upLim,w,x)                             
649        FUNCREF GenericQuadrature_proto fcn
650        Variable loLim,upLim    //limits of integration
651        Wave w                  //coefficients of function fcn(w,x)
652        Variable x      //x-value (q) for the calculation
653
654//**** The coefficient wave is passed into this function and straight through to the unsmeared model function
655// special case of integral limits that are identical
656        if(upLim == loLim)
657                return( fcn(w,x, loLim))
658        endif
659       
660// local variables
661        Variable nord,ii,va,vb,summ,yyy,zi
662        Variable answer,dum
663        String weightStr,zStr
664       
665        weightStr = "gauss76wt"
666        zStr = "gauss76z"
667        nord = 76
668               
669        if (WaveExists($weightStr) == 0) // wave reference is not valid,
670                Make/D/N=76 $weightStr,$zStr
671                Wave w76 = $weightStr
672                Wave z76 = $zStr                // wave references to pass
673                Make76GaussPoints(w76,z76)     
674        else
675                if(exists(weightStr) > 1)
676                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
677                endif
678                Wave w76 = $weightStr
679                Wave z76 = $zStr                // create the wave references
680        endif
681
682        //limits of integration are input to function
683        va = loLim
684        vb = upLim
685        // Using 76 Gauss points                   
686        // remember to index from 0,size-1
687
688        summ = 0.0              // initialize integral
689        ii=0                    // loop counter
690        do
691                // calculate Gauss points on integration interval (q-value for evaluation)
692                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
693                //calculate partial sum for the passed-in model function       
694                yyy = w76[ii] *  fcn(w,x,zi)                                           
695                summ += yyy             //add to the running total of the quadrature
696     ii+=1     
697        while (ii<nord)                         // end of loop over quadrature points
698   
699        // calculate value of integral to return
700        answer = (vb-va)/2.0*summ
701
702        Return (answer)
703End
704///////////////////////////////////////////////////////////////
705
706//////Resolution Smearing Utilities
707
708// To check for the existence of all waves needed for smearing
709// returns 1 if any waves are missing, 0 if all is OK
710Function ResolutionWavesMissing()
711
712        SVAR/Z sq = gSig_Q
713        SVAR/Z qb = gQ_bar
714        SVAR/Z sh = gShadow
715        SVAR/Z gQ = gQVals
716       
717        Variable err=0
718        if(!SVAR_Exists(sq) || !SVAR_Exists(qb) || !SVAR_Exists(sh) || !SVAR_Exists(gQ))
719                DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with LoadOneDData macro"
720                return(1)
721        endif
722       
723        if(WaveExists($sq) == 0)        //wave ref does not exist
724                err=1
725        endif
726        if(WaveExists($qb) == 0)        //wave ref does not exist
727                err=1
728        endif
729        if(WaveExists($sh) == 0)        //wave ref does not exist
730                err=1
731        endif
732        if(WaveExists($gQ) == 0)        //wave ref does not exist
733                err=1
734        endif
735       
736        if(err)
737                DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with 'Load SANS or USANS Data' macro"
738        endif
739       
740        return(err)
741end
742
743//////Resolution Smearing Utilities
744
745// To check for the existence of all waves needed for smearing
746// returns 1 if any waves are missing, 0 if all is OK
747// str passed in is the data folder containing the data
748//
749// 19 JUN07 using new data folder structure for loading
750// and resolution matrix
751Function ResolutionWavesMissingDF(str)
752        String str
753       
754        String DF="root:"+str+":"
755
756        WAVE/Z res = $(DF+str+"_res")
757       
758        if(!WaveExists(res))
759                DoAlert 0,"The resolution matrix is missing. Re-load experimental data with the 'Load SANS or USANS Data' macro"
760                return(1)
761        endif
762       
763        return(0)
764end
765
766///////////////////////////////////////////////////////////////
767
768// "backwards" wrapped to reduce redundant code
769// there are only 4 choices of N (5,10,20,76) for smearing
770//
771//
772// 4 MAR 2011
773// Note: In John's paper, he integrated the Gaussian to +/- 3 sigma and then renormalized
774//       to an integral of 1. This "truncated" gaussian was a somewhat better approximation
775//       to the triangular resolution function. Here, I integrate to +/- 3 sigma and
776//       do not renormalize the integral to 1. Hence the smeared calculation is 0.27% low.
777//       This is easily seen by smearing a constant value.
778//
779// Using 5 quadrature points is not recommended, as it doesn't normalize properly using .9973
780//  -- instead, it normalizes to 1.0084,
781//
782Function Smear_Model_N(fcn,w,x,resW,wi,zi,nord)
783        FUNCREF SANSModelAAO_proto fcn
784        Wave w                  //coefficients of function fcn(w,x)
785        Variable x      //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE
786        Wave resW               // Nx4 or NxN matrix of resolution
787        Wave wi         //weight wave
788        Wave zi         //abscissa wave
789        Variable nord           //order of integration
790
791        NVAR dQv = root:Packages:NIST:USANS_dQv
792
793// local variables
794        Variable ii,va,vb
795        Variable answer,i_shad,i_qbar,i_sigq,normalize=1
796
797        // current x point is the q-value for evaluation
798        //
799        // * for the input x, the resolution function waves are interpolated to get the correct values for
800        //  sigq, qbar and shad - since the model x-spacing may not be the same as
801        // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is
802        // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different
803        // from experimental data.
804        // **note** if the (x) passed in is the experimental q-values, these values are
805        // returned from the interpolation (as expected)
806
807        Make/O/D/N=(DimSize(resW, 0)) sigQ,qbar,shad,qvals
808        sigq = resW[p][0]               //std dev of resolution fn
809        qbar = resW[p][1]               //mean q-value
810        shad = resW[p][2]               //beamstop shadow factor
811        qvals = resW[p][3]      //q-values where R(q) is known
812
813        i_shad = interp(x,qvals,shad)
814        i_qbar = interp(x,qvals,qbar)
815        i_sigq = interp(x,qvals,sigq)
816                       
817// set up the integration
818// number of Gauss Quadrature points
819               
820        if (isSANSResolution(i_sigq))
821               
822                // end points of integration
823                // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
824                // +/- 3 sigq catches 99.73% of distrubution
825                // change limits (and spacing of zi) at each evaluation based on R()
826                //integration from va to vb
827       
828                // for +/- 3 sigma ONLY
829                if(nord == 5)
830                        normalize = 1.0057              //empirical correction, N=5 shouldn't be any different
831                else
832                        normalize = 0.9973
833                endif
834               
835                va = -3*i_sigq + i_qbar
836                if (va<0)
837                        va=0            //to avoid numerical error when  va<0 (-ve q-value)
838//                      Print "truncated Gaussian at nominal q = ",x
839                endif
840                vb = 3*i_sigq + i_qbar
841               
842               
843                // Using 20 Gauss points                   
844                ii=0                    // loop counter
845                // do the calculation as a single pass w/AAO function
846                Make/O/D/N=(nord) Resoln,yyy,xGauss
847                do
848                        // calculate Gauss points on integration interval (q-value for evaluation)
849                        xGauss[ii] = ( zi[ii]*(vb-va) + vb + va )/2.0
850                        // calculate resolution function at input q-value (use the interpolated values and zi)
851                        Resoln[ii] = i_shad/sqrt(2*pi*i_sigq*i_sigq)
852                        Resoln[ii] *= exp((-1*(xGauss[ii] - i_qbar)^2)/(2*i_sigq*i_sigq))
853                        ii+=1           
854                while (ii<nord)                         // end of loop over quadrature points
855               
856                fcn(w,yyy,xGauss)               //yyy is the return value as a wave
857               
858                yyy *= wi *Resoln               //multiply function by resolution and weights
859                // calculate value of integral to return
860                answer = (vb-va)/2.0*sum(yyy)
861                // all scaling, background addition... etc. is done in the model calculation
862       
863                        // renormalize to 1
864                        answer /= normalize
865        else
866                //smear with the USANS routine
867                // Make global string and local variables
868                // now data folder aware, necessary for GlobalFit = FULL path to wave   
869                String/G gTrap_coefStr = GetWavesDataFolder(w, 2 )     
870                Variable maxiter=20, tol=1e-4
871               
872                // set up limits for the integration
873                va=0
874                vb=abs(dQv)
875               
876                Variable/G gEvalQval = x
877               
878                // call qtrap to do actual work
879                answer = qtrap_USANS(fcn,va,vb,tol,maxiter)
880                answer /= vb
881               
882        endif
883       
884        //killing these waves is cleaner, but MUCH SLOWER
885//      Killwaves/Z Resoln,yyy,xGauss
886//      Killwaves/Z sigQ,qbar,shad,qvals
887        Return (answer)
888       
889End
890
891//resolution smearing, using only 5 Gauss points
892//
893//
894Function Smear_Model_5(fcn,w,x,answer,resW)                             
895        FUNCREF SANSModelAAO_proto fcn
896        Wave w                  //coefficients of function fcn(w,x)
897        Wave x  //x-value (q) for the calculation
898        Wave answer // ywave for calculation result
899        Wave resW               // Nx4 or NxN matrix of resolution
900        NVAR useTrap = root:Packages:NIST:USANSUseTrap
901
902        String weightStr,zStr
903        Variable nord=5
904       
905        if (dimsize(resW,1) > 4 && useTrap !=1)
906                if(dimsize(resW,1) != dimsize(answer,0) )
907                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
908                endif
909                //USANS Weighting matrix is present.
910                fcn(w,answer,x)
911       
912                MatrixOP/O  answer = resW x answer
913                //Duplicate/O answer,tmpMat
914                //MatrixOP/O answer = resW x tmpMat
915                Return(0)
916        else
917                weightStr = "gauss5wt"
918                zStr = "gauss5z"
919       
920        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
921                if (WaveExists($weightStr) == 0) // wave reference is not valid,
922                        Make/D/N=(nord) $weightStr,$zStr
923                        Wave weightW = $weightStr
924                        Wave abscissW = $zStr           // wave references to pass
925                        Make5GaussPoints(weightW,abscissW)     
926                else
927                        if(exists(weightStr) > 1)
928                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
929                        endif
930                        Wave weightW = $weightStr
931                        Wave abscissW = $zStr           // create the wave references
932                endif
933       
934//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
935                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer)
936               
937                Return (0)
938        endif
939       
940End
941
942//resolution smearing, using only 10 Gauss points
943//
944//
945Function Smear_Model_10(fcn,w,x,answer,resW)                           
946        FUNCREF SANSModelAAO_proto fcn
947        Wave w                  //coefficients of function fcn(w,x)
948        Wave x  //x-value (q) for the calculation
949        Wave answer // ywave for calculation result
950        Wave resW               // Nx4 or NxN matrix of resolution
951
952        String weightStr,zStr
953        Variable nord=10
954       
955        if (dimsize(resW,1) > 4)
956                if(dimsize(resW,1) != dimsize(answer,0) )
957                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
958                endif
959                //USANS Weighting matrix is present.
960                fcn(w,answer,x)
961       
962                MatrixOP/O  answer = resW x answer
963                //Duplicate/O answer,tmpMat
964                //MatrixOP/O answer = resW x tmpMat
965                Return(0)
966        else
967                weightStr = "gauss10wt"
968                zStr = "gauss10z"
969       
970        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
971                if (WaveExists($weightStr) == 0) // wave reference is not valid,
972                        Make/D/N=(nord) $weightStr,$zStr
973                        Wave weightW = $weightStr
974                        Wave abscissW = $zStr           // wave references to pass
975                        Make10GaussPoints(weightW,abscissW)     
976                else
977                        if(exists(weightStr) > 1)
978                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
979                        endif
980                        Wave weightW = $weightStr
981                        Wave abscissW = $zStr           // create the wave references
982                endif
983       
984//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
985                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer)
986
987                Return (0)
988        endif
989       
990End
991
992//
993//Smear_Model_20(SphereForm,s.coefW,s.yW,s.xW,s.resW)
994//
995//      Wave sigq               //std dev of resolution fn
996//      Wave qbar               //mean q-value
997//      Wave shad               //beamstop shadow factor
998//      Wave qvals      //q-values where R(q) is known
999//
1000Function Smear_Model_20(fcn,w,x,answer,resW)                           
1001        FUNCREF SANSModelAAO_proto fcn
1002        Wave w                  //coefficients of function fcn(w,x)
1003        Wave x  //x-value (q) for the calculation
1004        Wave answer // ywave for calculation result
1005        Wave resW               // Nx4 or NxN matrix of resolution
1006        NVAR useTrap = root:Packages:NIST:USANSUseTrap
1007
1008        String weightStr,zStr
1009        Variable nord=20
1010
1011       
1012        if (dimsize(resW,1) > 4 && useTrap != 1)
1013                if(dimsize(resW,1) != dimsize(answer,0) )
1014                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
1015                endif
1016                //USANS Weighting matrix is present.
1017                fcn(w,answer,x)
1018       
1019                MatrixOP/O  answer = resW x answer
1020                //Duplicate/O answer,tmpMat
1021                //MatrixOP/O answer = resW x tmpMat
1022                Return(0)
1023        else
1024                weightStr = "gauss20wt"
1025                zStr = "gauss20z"
1026       
1027        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
1028                if (WaveExists($weightStr) == 0) // wave reference is not valid,
1029                        Make/D/N=(nord) $weightStr,$zStr
1030                        Wave weightW = $weightStr
1031                        Wave abscissW = $zStr           // wave references to pass
1032                        Make20GaussPoints(weightW,abscissW)     
1033                else
1034                        if(exists(weightStr) > 1)
1035                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
1036                        endif
1037                        Wave weightW = $weightStr
1038                        Wave abscissW = $zStr           // create the wave references
1039                endif
1040       
1041//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
1042                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer)
1043
1044                Return (0)
1045        endif
1046       
1047End
1048///////////////////////////////////////////////////////////////
1049Function Smear_Model_76(fcn,w,x,answer,resW)                           
1050        FUNCREF SANSModelAAO_proto fcn
1051        Wave w                  //coefficients of function fcn(w,x)
1052        Wave x  //x-value (q) for the calculation
1053        Wave answer // ywave for calculation result
1054        Wave resW               // Nx4 or NxN matrix of resolution
1055        NVAR useTrap = root:Packages:NIST:USANSUseTrap
1056
1057        String weightStr,zStr
1058        Variable nord=76
1059       
1060        if (dimsize(resW,1) > 4  && useTrap != 1)
1061                if(dimsize(resW,1) != dimsize(answer,0) )
1062                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
1063                endif
1064                //USANS Weighting matrix is present.
1065                fcn(w,answer,x)
1066       
1067                MatrixOP/O  answer = resW x answer
1068                //Duplicate/O answer,tmpMat
1069                //MatrixOP/O answer = resW x tmpMat
1070                Return(0)
1071        else
1072                weightStr = "gauss76wt"
1073                zStr = "gauss76z"
1074       
1075        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
1076                if (WaveExists($weightStr) == 0) // wave reference is not valid,
1077                        Make/D/N=(nord) $weightStr,$zStr
1078                        Wave weightW = $weightStr
1079                        Wave abscissW = $zStr           // wave references to pass
1080                        Make76GaussPoints(weightW,abscissW)     
1081                else
1082                        if(exists(weightStr) > 1)
1083                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
1084                        endif
1085                        Wave weightW = $weightStr
1086                        Wave abscissW = $zStr           // create the wave references
1087                endif
1088       
1089//              answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
1090                Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer)
1091
1092                Return (0)
1093        endif
1094       
1095End
1096
1097
1098///////////////////////////////////////////////////////////////
1099
1100//typically, the first point (or any point) of sigQ is passed
1101// if negative, it's USANS data...
1102Function isSANSResolution(val)
1103        Variable val
1104        if(val>=0)
1105                return(1)               //true, SANS data
1106        else
1107                return(0)               //false, USANS
1108        endif
1109End
1110
1111Function GenericQuadrature_proto(w,x,dum)
1112        Wave w
1113        Variable x,dum
1114       
1115        Print "in GenericQuadrature_proto function"
1116        return(1)
1117end
1118 
1119// prototype function for smearing routines, Smear_Model_N
1120// and trapzd_USANS() and qtrap_USANS()
1121// it intentionally does nothing
1122Function SANSModel_proto(w,x)
1123        Wave w
1124        Variable x
1125       
1126        Print "in SANSModel_proto function"
1127        return(1)
1128end
1129
1130// prototype function for smearing routines, Smear_Model_N
1131// and trapzd_USANS() and qtrap_USANS()
1132// it intentionally does nothing
1133Function SANSModelAAO_proto(w,yw,xw)
1134        Wave w,yw,xw
1135       
1136        Print "in SANSModelAAO_proto function"
1137        return(1)
1138end
1139
1140// prototype function for fit wrapper
1141// it intentionally does nothing
1142Function SANSModelSTRUCT_proto(s)
1143        Struct ResSmearAAOStruct &s     
1144
1145        Print "in SANSModelSTRUCT_proto function"
1146        return(1)
1147end
1148
1149// prototype function for 2D smearing routine
1150ThreadSafe Function SANS_2D_ModelAAO_proto(w,zw,xw,yw)
1151        Wave w,zw,xw,yw
1152       
1153        Print "in SANSModelAAO_proto function"
1154        return(1)
1155end
1156
1157// prototype function for fit wrapper using 2D smearing
1158// not used (yet)
1159Function SANS_2D_ModelSTRUCT_proto(s)
1160        Struct ResSmear_2D_AAOStruct &s
1161
1162        Print "in SANSModelSTRUCT_proto function"
1163        return(1)
1164end
1165
1166//Numerical Recipes routine to calculate the nn(th) stage
1167//refinement of a trapezoid integration
1168//
1169//must be called sequentially from nn=1...n from qtrap()
1170// to cumulatively refine the integration value
1171//
1172// in the conversion:
1173// -- s was replaced with sVal and declared global (rather than static)
1174//  so that the nn-1 value would be available during the nn(th) call
1175//
1176// -- the specific coefficient wave for func() is passed in as a
1177//  global string (then converted to a wave reference), since
1178//  func() will eventually call sphereForm()
1179//
1180Function trapzd_USANS(fcn,aa,bb,nn)
1181        FUNCREF SANSModelAAO_proto fcn
1182        Variable aa,bb,nn
1183       
1184        Variable xx,tnm,summ,del
1185        Variable it,jj,arg1,arg2
1186        NVAR sVal=sVal          //calling function must initialize this global
1187        NVAR qval = gEvalQval
1188        SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
1189        Wave cw=$cwStr                 
1190        Variable temp=0
1191        Make/D/O/N=2 tmp_xw,tmp_yw
1192        if(nn==1)
1193                tmp_xw[0] = sqrt(qval^2 + aa^2)
1194                tmp_xw[1] = sqrt(qval^2 + bb^2)
1195                fcn(cw,tmp_yw,tmp_xw)
1196                temp = 0.5*(bb-aa)*(tmp_yw[0] + tmp_yw[1])
1197                sval = temp
1198                return(sVal)
1199        else
1200                it=1
1201                it= 2^(nn-2)  //done in NR with a bit shift <<=
1202                tnm = it
1203                del = (bb - aa)/tnm             //this is the spacing of points to add
1204                xx = aa+0.5*del
1205                summ=0
1206                for(jj=1;jj<=it;jj+=1)
1207                        tmp_xw = sqrt(qval^2 + xx^2)
1208                        fcn(cw,tmp_yw,tmp_xw)           //not the most efficient... but replaced by the matrix method
1209                        summ += tmp_yw[0]
1210                        xx += del
1211                endfor
1212                sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1213                return (sval)
1214        endif
1215        //KillWaves/Z tmp_xw,tmp_yw
1216End
1217
1218////Numerical Recipes routine to calculate the nn(th) stage
1219////refinement of a trapezoid integration
1220////
1221////must be called sequentially from nn=1...n from qtrap()
1222//// to cumulatively refine the integration value
1223////
1224//// in the conversion:
1225//// -- s was replaced with sVal and declared global (rather than static)
1226////  so that the nn-1 value would be available during the nn(th) call
1227////
1228//// -- the specific coefficient wave for func() is passed in as a
1229////  global string (then converted to a wave reference), since
1230////  func() will eventually call sphereForm()
1231////
1232//Function trapzd_USANS_point(fcn,aa,bb,nn)
1233//      FUNCREF SANSModel_proto fcn
1234//      Variable aa,bb,nn
1235//     
1236//      Variable xx,tnm,summ,del
1237//      Variable it,jj,arg1,arg2
1238//      NVAR sVal=sVal          //calling function must initialize this global
1239//      NVAR qval = gEvalQval
1240//      SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
1241//      Wave cw=$cwStr                 
1242//      Variable temp=0
1243//      if(nn==1)
1244//              arg1 = qval^2 + aa^2
1245//              arg2 = qval^2 + bb^2
1246//              temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2)))
1247//              sval = temp
1248//              return(sVal)
1249//      else
1250//              it=1
1251//              it= 2^(nn-2)  //done in NR with a bit shift <<=
1252//              tnm = it
1253//              del = (bb - aa)/tnm             //this is the spacing of points to add
1254//              xx = aa+0.5*del
1255//              summ=0
1256//              for(jj=1;jj<=it;jj+=1)
1257//                      arg1 = qval^2 + xx^2
1258//                      summ += fcn(cw,sqrt(arg1))
1259//                      xx += del
1260//              endfor
1261//              sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1262//              return (sval)
1263//      endif
1264//     
1265//End
1266
1267// Numerical Recipes routine to calculate the integral of a
1268// specified function, trapezoid rule is used to a user-specified
1269// level of refinement using sequential calls to trapzd()
1270//
1271// in NR, eps and maxIt were global, pass them in here...
1272// eps typically 1e-5
1273// maxit typically 20
1274Function qtrap_USANS(fcn,aa,bb,eps,maxIt)
1275        FUNCREF SANSModelAAO_proto fcn
1276        Variable aa,bb,eps,maxit
1277       
1278        Variable/G sVal=0               //create and initialize what trapzd will return
1279        Variable jj,ss,olds
1280       
1281        olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1282        for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1283                ss = trapzd_USANS(fcn,aa,bb,jj)
1284                if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1285                        return ss
1286                endif
1287                if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1288                        return ss
1289                endif
1290                olds = ss
1291        endfor
1292
1293        Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1294        return(ss)              //should never get here if function is well-behaved
1295       
1296End     
1297
1298//// Numerical Recipes routine to calculate the integral of a
1299//// specified function, trapezoid rule is used to a user-specified
1300//// level of refinement using sequential calls to trapzd()
1301////
1302//// in NR, eps and maxIt were global, pass them in here...
1303//// eps typically 1e-5
1304//// maxit typically 20
1305//Function qtrap_USANS_point(fcn,aa,bb,eps,maxIt)
1306//      FUNCREF SANSModel_proto fcn
1307//      Variable aa,bb,eps,maxit
1308//     
1309//      Variable/G sVal=0               //create and initialize what trapzd will return
1310//      Variable jj,ss,olds
1311//     
1312//      olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1313//      for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1314//              ss = trapzd_USANS_point(fcn,aa,bb,jj)
1315//              if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1316//                      return ss
1317//              endif
1318//              if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1319//                      return ss
1320//              endif
1321//              olds = ss
1322//      endfor
1323//
1324//      Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1325//      return(ss)              //should never get here if function is well-behaved
1326//     
1327//End   
1328
1329Proc RRW()
1330        ResetResolutionWaves()
1331End
1332
1333//utility procedures that are currently untied to any actions, although useful...
1334Proc ResetResolutionWaves(str)
1335        String Str
1336        Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","")
1337
1338
1339        Abort "This function is not data floder aware and does nothing..."
1340               
1341        String/G root:gQvals = str[0,strlen(str)-3]+"_q"
1342        String/G root:gSig_Q = str[0,strlen(str)-3]+"sq"
1343        String/G root:gQ_bar = str[0,strlen(str)-3]+"qb"
1344        String/G root:gShadow = str[0,strlen(str)-3]+"fs"
1345       
1346        //touch everything to make sure that the dependencies are
1347        //properly updated - especially the $gQvals reference in the
1348        // dependency assignment
1349        fKillDependencies("Smear*")
1350       
1351        //replace the q-values and intensity (so they're the right length)
1352        fResetSmearedModels("Smear*",root:gQvals)
1353       
1354        fRestoreDependencies("Smear*")
1355End
1356
1357// pass "*" as the matchString to do ALL dependent waves
1358// or "abc*" to get just the matching waves
1359//
1360Function fKillDependencies(matchStr)
1361        String matchStr
1362
1363        String str=WaveList(matchStr, ";", "" ),item,formula
1364        Variable ii
1365       
1366        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1367                item = StringFromList(ii, str ,";")
1368                formula = GetFormula($item)
1369                if(cmpstr("",formula)!=0)
1370                        Printf "wave %s had the formula %s removed\r",item,formula
1371                        Note $item, "FORMULA:"+formula
1372                        SetFormula $item, ""                    //clears the formula
1373                endif
1374        endfor
1375        return(0)
1376end
1377
1378// pass "*" as the matchString to do ALL dependent waves
1379// or "abc*" to get just the matching waves
1380//
1381Function fRestoreDependencies(matchStr)
1382        String matchStr
1383
1384        String str=WaveList(matchStr, ";", "" ),item,formula
1385        Variable ii
1386       
1387        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1388                item = StringFromList(ii, str ,";")
1389                formula = StringByKey("FORMULA", note($item),":",";")
1390                if(cmpstr("",formula)!=0)
1391                        Printf "wave %s had the formula %s restored\r",item,formula
1392                        Note/K $item
1393                        SetFormula $item, formula                       //restores the formula
1394                endif
1395        endfor
1396        return(0)
1397end
1398
1399Function fResetSmearedModels(matchStr,qStr)
1400        String matchStr,qStr
1401
1402        Duplicate/O $qStr root:smeared_qvals   
1403       
1404        String str=WaveList(matchStr, ";", "" ),item,formula
1405        Variable ii
1406       
1407        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1408                item = StringFromList(ii, str ,";")
1409                formula = StringByKey("FORMULA", note($item),":",";")
1410                if(cmpstr("",formula)!=0)
1411                        Printf "wave %s has been duplicated to gQvals\r",item
1412                        Duplicate/O $qStr $item
1413                        Note $item, "FORMULA:"+formula          //be sure to keep the formula note
1414                        Print "   and the formula is",formula
1415                endif
1416        endfor
1417        return(0)
1418end
1419
1420
1421////
1422//// moved from RawWindowHook to here - where the Q calculations are available to
1423//   reduction and analysis
1424//
1425
1426//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1427Threadsafe Function FindPhi(vx,vy)
1428        variable vx,vy
1429       
1430        variable phi
1431       
1432        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
1433       
1434        // special cases
1435        if(vx==0 && vy > 0)
1436                return(pi/2)
1437        endif
1438        if(vx==0 && vy < 0)
1439                return(3*pi/2)
1440        endif
1441        if(vx >= 0 && vy == 0)
1442                return(0)
1443        endif
1444        if(vx < 0 && vy == 0)
1445                return(pi)
1446        endif
1447       
1448       
1449        if(vx > 0 && vy > 0)
1450                return(phi)
1451        endif
1452        if(vx < 0 && vy > 0)
1453                return(phi + pi)
1454        endif
1455        if(vx < 0 && vy < 0)
1456                return(phi + pi)
1457        endif
1458        if( vx > 0 && vy < 0)
1459                return(phi + 2*pi)
1460        endif
1461       
1462        return(phi)
1463end
1464
1465       
1466//function to calculate the overall q-value, given all of the necesary trig inputs
1467//NOTE: detector locations passed in are pixels = 0.5cm real space on the detector
1468//and are in detector coordinates (1,128) rather than axis values
1469//the pixel locations need not be integers, reals are ok inputs
1470//sdd is in meters
1471//wavelength is in Angstroms
1472//
1473//returned magnitude of Q is in 1/Angstroms
1474//
1475Function CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1476        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize
1477       
1478        Variable dx,dy,qval,two_theta,dist
1479       
1480        Variable pixSizeX=pixSize
1481        Variable pixSizeY=pixSize
1482       
1483        sdd *=100               //convert to cm
1484        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
1485        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
1486        dist = sqrt(dx^2 + dy^2)
1487       
1488        two_theta = atan(dist/sdd)
1489
1490        qval = 4*Pi/lam*sin(two_theta/2)
1491       
1492        return qval
1493End
1494
1495//calculates just the q-value in the x-direction on the detector
1496//input/output is the same as CalcQval()
1497//ALL inputs are in detector coordinates
1498//
1499//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
1500//sdd is in meters
1501//wavelength is in Angstroms
1502//
1503// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
1504// now properly accounts for qz
1505//
1506Function CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1507        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize
1508
1509        Variable qx,qval,phi,dx,dy,dist,two_theta
1510       
1511        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1512       
1513        sdd *=100               //convert to cm
1514        dx = (xaxval - xctr)*pixSize            //delta x in cm
1515        dy = (yaxval - yctr)*pixSize            //delta y in cm
1516        phi = FindPhi(dx,dy)
1517       
1518        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
1519        dist = sqrt(dx^2 + dy^2)
1520        two_theta = atan(dist/sdd)
1521
1522        qx = qval*cos(two_theta/2)*cos(phi)
1523       
1524        return qx
1525End
1526
1527//calculates just the q-value in the y-direction on the detector
1528//input/output is the same as CalcQval()
1529//ALL inputs are in detector coordinates
1530//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
1531//sdd is in meters
1532//wavelength is in Angstroms
1533//
1534// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
1535// now properly accounts for qz
1536//
1537Function CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1538        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize
1539       
1540        Variable dy,qval,dx,phi,qy,dist,two_theta
1541       
1542        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1543       
1544        sdd *=100               //convert to cm
1545        dx = (xaxval - xctr)*pixSize            //delta x in cm
1546        dy = (yaxval - yctr)*pixSize            //delta y in cm
1547        phi = FindPhi(dx,dy)
1548       
1549        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
1550        dist = sqrt(dx^2 + dy^2)
1551        two_theta = atan(dist/sdd)
1552       
1553        qy = qval*cos(two_theta/2)*sin(phi)
1554       
1555        return qy
1556End
1557
1558//calculates just the z-component of the q-vector, not measured on the detector
1559//input/output is the same as CalcQval()
1560//ALL inputs are in detector coordinates
1561//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
1562//sdd is in meters
1563//wavelength is in Angstroms
1564//
1565// not actually used, but here for completeness if anyone asks
1566//
1567Function CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1568        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize
1569       
1570        Variable dy,qval,dx,phi,qz,dist,two_theta
1571       
1572        qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize)
1573       
1574        sdd *=100               //convert to cm
1575       
1576        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
1577        dx = (xaxval - xctr)*pixSize            //delta x in cm
1578        dy = (yaxval - yctr)*pixSize            //delta y in cm
1579        dist = sqrt(dx^2 + dy^2)
1580        two_theta = atan(dist/sdd)
1581       
1582        qz = qval*sin(two_theta/2)
1583       
1584        return qz
1585End
1586
1587//for command-line testing, replace the function declaration
1588//Function FindQxQy(qq,phi)
1589//      Variable qq,phi
1590//      Variable qx,qy
1591//
1592//
1593ThreadSafe Function FindQxQy(qq,phi,qx,qy)
1594        Variable qq,phi,&qx,&qy
1595
1596        qx = sqrt(qq^2/(1+tan(phi)*tan(phi)))
1597        qy = qx*tan(phi)
1598       
1599        if(phi >= 0 && phi <= pi/2)
1600                qx = abs(qx)
1601                qy = abs(qy)
1602        endif
1603       
1604        if(phi > pi/2 && phi <= pi)
1605                qx = -abs(qx)
1606                qy = abs(qy)
1607        endif
1608       
1609        if(phi > pi && phi <= pi*3/2)
1610                qx = -abs(qx)
1611                qy = -abs(qy)
1612        endif
1613       
1614        if(phi > pi*3/2 && phi < 2*pi)
1615                qx = abs(qx)
1616                qy = -abs(qy)
1617        endif   
1618       
1619       
1620//      Print "recalculated qx,qy,q = ",qx,qy,sqrt(qx*qx+qy*qy)
1621       
1622        return(0)
1623end
1624
1625
1626// 7 MAR 2011 SRK
1627//
1628// calculate the resolution smearing AAO
1629//
1630// - many of the form factor calculations are threaded, so they benefit
1631// from being passed large numbers of q-values at once, rather than suffering the
1632// overhead penalty of setting up threads.
1633//
1634// In general, single integral functions benefit from this, multiple integrals not so much.
1635// As an example, a fit using SmearedCylinderForm took 4.3s passing nord (=20) q-values
1636// at a time, but only 1.1s by passing all (Nq*nord) q-values at once. For Cyl_polyRad,
1637// the difference was not so large, 16.2s vs. 11.9s. This is due to CylPolyRad being a
1638// double integral and slow enough of a calculation that passing even 20 points at once
1639// provides some speedup.
1640//
1641//// APRIL 2011 *** this function is now cursor aware. The whole input x-wave is interpolated
1642//
1643//
1644// 4 MAR 2011 SRK
1645// Note: In John's paper, he integrated the Gaussian to +/- 3 sigma and then renormalized
1646//       to an integral of 1. This "truncated" gaussian was a somewhat better approximation
1647//       to the triangular resolution function. Here, I integrate to +/- 3 sigma and
1648//       do not renormalize the integral to 1. Hence the smeared calculation is 0.27% low.
1649//       This is easily seen by smearing a constant value.
1650//
1651// Using 5 quadrature points is not recommended, as it doesn't normalize properly using .9973
1652//  -- instead, it normalizes to 1.0084.
1653//
1654Function Smear_Model_N_AAO(fcn,w,x,resW,wi,zi,nord,sm_ans)
1655        FUNCREF SANSModelAAO_proto fcn
1656        Wave w                  //coefficients of function fcn(w,x)
1657        WAVE x                  //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE
1658        Wave resW               // Nx4 or NxN matrix of resolution
1659        Wave wi         //weight wave
1660        Wave zi         //abscissa wave
1661        Variable nord           //order of integration
1662        Wave sm_ans             // wave returned with the smeared model
1663
1664        NVAR dQv = root:Packages:NIST:USANS_dQv
1665
1666// local variables
1667        Variable ii,jj
1668        Variable normalize=1
1669        Variable nTot,num,block_sum
1670
1671
1672        // current x point is the q-value for evaluation
1673        //
1674        // * for the input x, the resolution function waves are interpolated to get the correct values for
1675        //  sigq, qbar and shad - since the model x-spacing may not be the same as
1676        // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is
1677        // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different
1678        // from experimental data.
1679        // **note** if the (x) passed in is the experimental q-values, these values are
1680        // returned from the interpolation (as expected)
1681
1682        Make/O/D/N=(numpnts(x)) sigQ,qbar,shad,qvals,va,vb
1683        Make/O/D/N=(DimSize(resW, 0)) tmpsigQ,tmpqbar,tmpshad,tmpqvals
1684        tmpsigq = resW[p][0]            //std dev of resolution fn
1685        tmpqbar = resW[p][1]            //mean q-value
1686        tmpshad = resW[p][2]            //beamstop shadow factor
1687        tmpqvals = resW[p][3]   //q-values where R(q) is known
1688
1689        //interpolate the whole input x-wave to make sure that the resolution and input x are in sync if cursors are used
1690        shad = interp(x,tmpqvals,tmpshad)
1691        qbar = interp(x,tmpqvals,tmpqbar)
1692        sigq = interp(x,tmpqvals,tmpsigq)
1693       
1694        // if USANS data, handle separately
1695        // -- but this would only ever be used if the calculation was forced to use trapezoid integration
1696        if ( ! isSANSResolution(sigq[0]) )
1697                        //smear with the USANS routine
1698                // Make global string and local variables
1699                // now data folder aware, necessary for GlobalFit = FULL path to wave   
1700                String/G gTrap_coefStr = GetWavesDataFolder(w, 2 )     
1701                Variable maxiter=20, tol=1e-4,uva,uvb
1702               
1703                num=numpnts(x)
1704                // set up limits for the integration
1705                uva=0
1706                uvb=abs(dQv)
1707               
1708                //loop over the q-values
1709                for(jj=0;jj<num;jj+=1)
1710                        Variable/G gEvalQval = x[jj]
1711
1712                        // call qtrap to do actual work
1713                        sm_ans[jj] = qtrap_USANS(fcn,uva,uvb,tol,maxiter)
1714                        sm_ans[jj] /= (uvb - uva)
1715                endfor
1716               
1717                return(0)       
1718        endif
1719
1720
1721// now the idea is to calculate a long vector of all of the zi's (Nq * nord)
1722// and pass these AAO to the AAO function, to make the most use of the threading
1723// passing repeated short lengths of q to the function can actually be slower
1724// due to the overhead.
1725       
1726        num = numpnts(x)
1727        nTot = nord*num
1728       
1729        Make/O/D/N=(nTot) Resoln,yyy,xGauss,wts
1730
1731        //loop over q
1732        for(jj=0;jj<num;jj+=1)
1733       
1734                //for each q, set up the integration range
1735                // end points of integration limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
1736                // +/- 3 sigq catches 99.73% of distrubution
1737                // change limits (and spacing of zi) at each evaluation based on R()
1738                //integration from va to vb
1739
1740                va[jj] = -3*sigq[jj] + qbar[jj]
1741                if (va[jj]<0)
1742                        va[jj]=0                //to avoid numerical error when  va<0 (-ve q-value)
1743//                      Print "truncated Gaussian at nominal q = ",x
1744                endif
1745                vb[jj] = 3*sigq[jj] + qbar[jj]
1746       
1747                // loop over the Gauss points
1748                for(ii=0;ii<nord;ii+=1)
1749                        // calculate Gauss points on integration interval (q-value for evaluation)
1750                        xGauss[nord*jj+ii] = ( zi[ii]*(vb[jj]-va[jj]) + vb[jj] + va[jj] )/2.0
1751                        // calculate resolution function at input q-value (use the interpolated values and zi)
1752                        Resoln[nord*jj+ii] = shad[jj]/sqrt(2*pi*sigq[jj]*sigq[jj])
1753                        Resoln[nord*jj+ii] *= exp((-1*(xGauss[nord*jj+ii] - qbar[jj])^2)/(2*sigq[jj]*sigq[jj]))
1754//                      Resoln[nord*jj+ii] *= exp((-1*(xGauss[nord*jj+ii] - qvals[jj])^2)/(2*sigq[jj]*sigq[jj]))                //WRONG, but just for testing
1755                        // carry a copy of the weights
1756                        wts[nord*jj+ii] = wi[ii]
1757                endfor          // end of loop over quadrature points
1758       
1759        endfor          //loop over q
1760       
1761        //calculate AAO
1762        yyy = 0
1763        fcn(w,yyy,xGauss)               //yyy is the return value as a wave
1764
1765        //multiply by weights
1766        yyy *= wts*Resoln               //multiply function by resolution and weights
1767       
1768        //sum up blockwise to get the final answer
1769        for(jj=0;jj<num;jj+=1)
1770                block_sum = 0
1771                for(ii=0;ii<nord;ii+=1)
1772                        block_sum += yyy[nord*jj+ii]
1773                endfor
1774                sm_ans[jj] = (vb[jj]-va[jj])/2.0*block_sum
1775        endfor
1776       
1777       
1778        // then normalize for +/- 3 sigma ONLY
1779        if(nord == 5)
1780                normalize = 1.0057              //empirical correction, N=5 shouldn't be any different
1781        else
1782                normalize = 0.9973
1783        endif
1784       
1785        sm_ans /= normalize
1786       
1787        KillWaves/Z tmpsigQ,tmpqbar,tmpshad,tmpqvals
1788       
1789        return(0)
1790       
1791End
1792
Note: See TracBrowser for help on using the repository browser.