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

Last change on this file since 912 was 828, checked in by srkline, 12 years ago

adding some loaders to the macros menu - strictly for testing so that some of the newest procedures can be easily tested.
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//       now correctly renormalize the integral to 1. Hence the smeared calculation in the past was 0.27% low.
777//       Confimation of the integral 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.