source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v4.00/GaussUtils_v40.ipf @ 305

Last change on this file since 305 was 305, checked in by ajj, 14 years ago

Adding option to force use of Trapezoidal integration in USANS smearing for debug purposes.

Setting global variable root:NIST:Packages:USANSUseTrap = 1 will turn on the trapezoidal routine.

File size: 35.2 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=4.00
3#pragma IgorVersion=6.0
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// utility used in the "PlotSmeared...() macros to get a list of data folders
64//
65//1:    Waves.
66//2:    Numeric variables.
67//3:    String variables.
68//4:    Data folders.
69Function/S GetAList(type)
70        Variable type
71       
72        SetDataFolder root:
73       
74        String objName,str=""
75        Variable index = 0
76        do
77                objName = GetIndexedObjName(":", type, index)
78                if (strlen(objName) == 0)
79                        break
80                endif
81                //Print objName
82                str += objName + ";"
83                index += 1
84        while(1)
85       
86        // remove myGlobals, Packages from the folder list
87        if(type==4)
88                str = RemoveFromList("myGlobals", str , ";" )
89                str = RemoveFromList("Packages", str, ";")
90                str = RemoveFromList("AutoFit", str, ";")
91        endif
92       
93        return(str)
94End
95
96
97Function Make5GaussPoints(w5,z5)
98        Wave w5,z5
99
100//      printf  "in make Gauss Pts\r"
101
102        z5[0] = -.906179845938664
103        z5[1] = -.538469310105683
104        z5[2] = -.0000000000000
105        z5[3] = .538469310105683
106        z5[4] = .906179845938664
107
108        w5[0] = .236926885056189
109        w5[1] = .478628670499366
110        w5[2] = .56888888888889
111        w5[3] = .478628670499366
112        w5[4] = .236926885056189
113
114//          printf "w[0],z[0] = %g %g\r", w5[0],z5[0]
115End
116
117Function Make10GaussPoints(w10,z10)
118        Wave w10,z10
119
120//      printf  "in make Gauss Pts\r"
121        z10[0] = -.973906528517172
122        z10[1] = -.865063366688985
123        z10[2] = -.679409568299024
124        z10[3] = -.433395394129247
125        z10[4] = -.148874338981631
126        z10[5] = .148874338981631
127        z10[6] = .433395394129247
128        z10[7] = .679409568299024
129        z10[8] = .865063366688985
130        z10[9] = .973906528517172
131       
132        w10[0] = .066671344308688
133        w10[1] = 0.149451349150581
134        w10[2] = 0.219086362515982
135        w10[3] = .269266719309996
136        w10[4] = 0.295524224714753
137        w10[5] = 0.295524224714753
138        w10[6] = .269266719309996
139        w10[7] = 0.219086362515982
140        w10[8] = 0.149451349150581
141        w10[9] = .066671344308688
142       
143//          printf "w[0],z[0] = %g %g\r", w10[0],z10[0]
144End
145
146Function Make20GaussPoints(w20,z20)
147        Wave w20,z20
148
149//      printf  "in make Gauss Pts\r"
150
151         z20[0] = -.993128599185095
152         z20[1] =  -.963971927277914
153         z20[2] =    -.912234428251326
154         z20[3] =    -.839116971822219
155         z20[4] =   -.746331906460151
156         z20[5] =   -.636053680726515
157         z20[6] =   -.510867001950827
158         z20[7] =     -.37370608871542
159         z20[8] =     -.227785851141645
160         z20[9] =     -.076526521133497
161         z20[10] =     .0765265211334973
162         z20[11] =     .227785851141645
163         z20[12] =     .37370608871542
164         z20[13] =     .510867001950827
165         z20[14] =     .636053680726515
166         z20[15] =     .746331906460151
167         z20[16] =     .839116971822219
168         z20[17] =     .912234428251326
169         z20[18] =    .963971927277914
170         z20[19] =  .993128599185095
171           
172        w20[0] =  .0176140071391521
173        w20[1] =  .0406014298003869
174        w20[2] =      .0626720483341091
175        w20[3] =      .0832767415767047
176        w20[4] =     .10193011981724
177        w20[5] =      .118194531961518
178        w20[6] =      .131688638449177
179        w20[7] =      .142096109318382
180        w20[8] =      .149172986472604
181        w20[9] =      .152753387130726
182        w20[10] =      .152753387130726
183        w20[11] =      .149172986472604
184        w20[12] =      .142096109318382
185        w20[13] =     .131688638449177
186        w20[14] =     .118194531961518
187        w20[15] =     .10193011981724
188        w20[16] =     .0832767415767047
189        w20[17] =     .0626720483341091
190        w20[18] =     .0406014298003869
191        w20[19] =     .0176140071391521
192//          printf "w[0],z[0] = %g %g\r", w20[0],z20[0]
193End
194
195Function Make76GaussPoints(w76,z76)
196        Wave w76,z76
197
198//      printf  "in make Gauss Pts\r"
199       
200                z76[0] = .999505948362153*(-1.0)
201            z76[75] = -.999505948362153*(-1.0)
202            z76[1] = .997397786355355*(-1.0)
203            z76[74] = -.997397786355355*(-1.0)
204            z76[2] = .993608772723527*(-1.0)
205            z76[73] = -.993608772723527*(-1.0)
206            z76[3] = .988144453359837*(-1.0)
207            z76[72] = -.988144453359837*(-1.0)
208            z76[4] = .981013938975656*(-1.0)
209            z76[71] = -.981013938975656*(-1.0)
210            z76[5] = .972229228520377*(-1.0)
211            z76[70] = -.972229228520377*(-1.0)
212            z76[6] = .961805126758768*(-1.0)
213            z76[69] = -.961805126758768*(-1.0)
214            z76[7] = .949759207710896*(-1.0)
215            z76[68] = -.949759207710896*(-1.0)
216            z76[8] = .936111781934811*(-1.0)
217            z76[67] = -.936111781934811*(-1.0)
218            z76[9] = .92088586125215*(-1.0)
219            z76[66] = -.92088586125215*(-1.0)
220            z76[10] = .904107119545567*(-1.0)
221            z76[65] = -.904107119545567*(-1.0)
222            z76[11] = .885803849292083*(-1.0)
223            z76[64] = -.885803849292083*(-1.0)
224            z76[12] = .866006913771982*(-1.0)
225            z76[63] = -.866006913771982*(-1.0)
226            z76[13] = .844749694983342*(-1.0)
227            z76[62] = -.844749694983342*(-1.0)
228            z76[14] = .822068037328975*(-1.0)
229            z76[61] = -.822068037328975*(-1.0)
230            z76[15] = .7980001871612*(-1.0)
231            z76[60] = -.7980001871612*(-1.0)
232            z76[16] = .77258672828181*(-1.0)
233            z76[59] = -.77258672828181*(-1.0)
234            z76[17] = .74587051350361*(-1.0)
235            z76[58] = -.74587051350361*(-1.0)
236            z76[18] = .717896592387704*(-1.0)
237            z76[57] = -.717896592387704*(-1.0)
238            z76[19] = .688712135277641*(-1.0)
239            z76[56] = -.688712135277641*(-1.0)
240            z76[20] = .658366353758143*(-1.0)
241            z76[55] = -.658366353758143*(-1.0)
242            z76[21] = .626910417672267*(-1.0)
243            z76[54] = -.626910417672267*(-1.0)
244            z76[22] = .594397368836793*(-1.0)
245            z76[53] = -.594397368836793*(-1.0)
246            z76[23] = .560882031601237*(-1.0)
247            z76[52] = -.560882031601237*(-1.0)
248            z76[24] = .526420920401243*(-1.0)
249            z76[51] = -.526420920401243*(-1.0)
250            z76[25] = .491072144462194*(-1.0)
251            z76[50] = -.491072144462194*(-1.0)
252            z76[26] = .454895309813726*(-1.0)
253            z76[49] = -.454895309813726*(-1.0)
254            z76[27] = .417951418780327*(-1.0)
255            z76[48] = -.417951418780327*(-1.0)
256            z76[28] = .380302767117504*(-1.0)
257            z76[47] = -.380302767117504*(-1.0)
258            z76[29] = .342012838966962*(-1.0)
259            z76[46] = -.342012838966962*(-1.0)
260            z76[30] = .303146199807908*(-1.0)
261            z76[45] = -.303146199807908*(-1.0)
262            z76[31] = .263768387584994*(-1.0)
263            z76[44] = -.263768387584994*(-1.0)
264            z76[32] = .223945802196474*(-1.0)
265            z76[43] = -.223945802196474*(-1.0)
266            z76[33] = .183745593528914*(-1.0)
267            z76[42] = -.183745593528914*(-1.0)
268            z76[34] = .143235548227268*(-1.0)
269            z76[41] = -.143235548227268*(-1.0)
270            z76[35] = .102483975391227*(-1.0)
271            z76[40] = -.102483975391227*(-1.0)
272            z76[36] = .0615595913906112*(-1.0)
273            z76[39] = -.0615595913906112*(-1.0)
274            z76[37] = .0205314039939986*(-1.0)
275            z76[38] = -.0205314039939986*(-1.0)
276           
277                w76[0] =  .00126779163408536
278                w76[75] = .00126779163408536
279                w76[1] =  .00294910295364247
280            w76[74] = .00294910295364247
281            w76[2] = .00462793522803742
282            w76[73] =  .00462793522803742
283            w76[3] = .00629918049732845
284            w76[72] = .00629918049732845
285            w76[4] = .00795984747723973
286            w76[71] = .00795984747723973
287            w76[5] = .00960710541471375
288            w76[70] =  .00960710541471375
289            w76[6] = .0112381685696677
290            w76[69] = .0112381685696677
291            w76[7] =  .0128502838475101
292            w76[68] = .0128502838475101
293            w76[8] = .0144407317482767
294            w76[67] =  .0144407317482767
295            w76[9] = .0160068299122486
296            w76[66] = .0160068299122486
297            w76[10] = .0175459372914742
298            w76[65] = .0175459372914742
299            w76[11] = .0190554584671906
300            w76[64] = .0190554584671906
301            w76[12] = .020532847967908
302            w76[63] = .020532847967908
303            w76[13] = .0219756145344162
304            w76[62] = .0219756145344162
305            w76[14] = .0233813253070112
306            w76[61] = .0233813253070112
307            w76[15] = .0247476099206597
308            w76[60] = .0247476099206597
309            w76[16] = .026072164497986
310            w76[59] = .026072164497986
311            w76[17] = .0273527555318275
312            w76[58] = .0273527555318275
313            w76[18] = .028587223650054
314            w76[57] = .028587223650054
315            w76[19] = .029773487255905
316            w76[56] = .029773487255905
317            w76[20] = .0309095460374916
318            w76[55] = .0309095460374916
319            w76[21] = .0319934843404216
320            w76[54] = .0319934843404216
321            w76[22] = .0330234743977917
322            w76[53] = .0330234743977917
323            w76[23] = .0339977794120564
324            w76[52] = .0339977794120564
325            w76[24] = .0349147564835508
326            w76[51] = .0349147564835508
327            w76[25] = .0357728593807139
328            w76[50] = .0357728593807139
329            w76[26] = .0365706411473296
330            w76[49] = .0365706411473296
331            w76[27] = .0373067565423816
332            w76[48] = .0373067565423816
333            w76[28] = .0379799643084053
334            w76[47] = .0379799643084053
335            w76[29] = .0385891292645067
336            w76[46] = .0385891292645067
337            w76[30] = .0391332242205184
338            w76[45] = .0391332242205184
339            w76[31] = .0396113317090621
340            w76[44] = .0396113317090621
341            w76[32] = .0400226455325968
342            w76[43] = .0400226455325968
343            w76[33] = .040366472122844
344            w76[42] = .040366472122844
345            w76[34] = .0406422317102947
346            w76[41] = .0406422317102947
347            w76[35] = .0408494593018285
348            w76[40] = .0408494593018285
349            w76[36] = .040987805464794
350            w76[39] = .040987805464794
351            w76[37] = .0410570369162294
352            w76[38] = .0410570369162294
353
354End             //Make76GaussPoints()
355
356////////////
357Function IntegrateFn5(fcn,loLim,upLim,w,x)                             
358        FUNCREF GenericQuadrature_proto fcn
359        Variable loLim,upLim    //limits of integration
360        Wave w                  //coefficients of function fcn(w,x)
361        Variable x      //x-value (q) for the calculation
362
363// local variables
364        Variable nord,ii,va,vb,summ,yyy,zi
365        Variable answer,dum
366        String weightStr,zStr
367       
368        weightStr = "gauss5wt"
369        zStr = "gauss5z"
370        nord = 5
371               
372        if (WaveExists($weightStr) == 0) // wave reference is not valid,
373                Make/D/N=5 $weightStr,$zStr
374                Wave w5 = $weightStr
375                Wave z5 = $zStr         // wave references to pass
376                Make5GaussPoints(w5,z5)
377        else
378                if(exists(weightStr) > 1)
379                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
380                endif
381                Wave w5 = $weightStr
382                Wave z5 = $zStr         // create the wave references
383        endif
384
385        //limits of integration are input to function
386        va = loLim
387        vb = upLim
388        // Using 5 Gauss points             
389        // remember to index from 0,size-1
390
391        summ = 0.0              // initialize integral
392        ii=0                    // loop counter
393        do
394                // calculate Gauss points on integration interval (q-value for evaluation)
395                zi = ( z5[ii]*(vb-va) + vb + va )/2.0
396                //calculate partial sum for the passed-in model function       
397                yyy = w5[ii] *  fcn(w,x,zi)                                             
398                summ += yyy             //add to the running total of the quadrature
399        ii+=1           
400        while (ii<nord)                         // end of loop over quadrature points
401   
402        // calculate value of integral to return
403        answer = (vb-va)/2.0*summ
404
405        Return (answer)
406End
407///////////////////////////////////////////////////////////////
408
409////////////
410Function IntegrateFn10(fcn,loLim,upLim,w,x)                             
411        FUNCREF GenericQuadrature_proto fcn
412        Variable loLim,upLim    //limits of integration
413        Wave w                  //coefficients of function fcn(w,x)
414        Variable x      //x-value (q) for the calculation
415
416// local variables
417        Variable nord,ii,va,vb,summ,yyy,zi
418        Variable answer,dum
419        String weightStr,zStr
420       
421        weightStr = "gauss10wt"
422        zStr = "gauss10z"
423        nord = 10
424               
425        if (WaveExists($weightStr) == 0) // wave reference is not valid,
426                Make/D/N=10 $weightStr,$zStr
427                Wave w10 = $weightStr
428                Wave z10 = $zStr                // wave references to pass
429                Make10GaussPoints(w10,z10)     
430        else
431                if(exists(weightStr) > 1)
432                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
433                endif
434                Wave w10 = $weightStr
435                Wave z10 = $zStr                // create the wave references
436        endif
437
438        //limits of integration are input to function
439        va = loLim
440        vb = upLim
441        // Using 10 Gauss points                   
442        // remember to index from 0,size-1
443
444        summ = 0.0              // initialize integral
445        ii=0                    // loop counter
446        do
447                // calculate Gauss points on integration interval (q-value for evaluation)
448                zi = ( z10[ii]*(vb-va) + vb + va )/2.0
449                //calculate partial sum for the passed-in model function       
450                yyy = w10[ii] *  fcn(w,x,zi)                                           
451                summ += yyy             //add to the running total of the quadrature
452        ii+=1           
453        while (ii<nord)                         // end of loop over quadrature points
454   
455        // calculate value of integral to return
456        answer = (vb-va)/2.0*summ
457
458        Return (answer)
459End
460///////////////////////////////////////////////////////////////
461
462////////////
463Function IntegrateFn20(fcn,loLim,upLim,w,x)                             
464        FUNCREF GenericQuadrature_proto fcn
465        Variable loLim,upLim    //limits of integration
466        Wave w                  //coefficients of function fcn(w,x)
467        Variable x      //x-value (q) for the calculation
468
469// local variables
470        Variable nord,ii,va,vb,summ,yyy,zi
471        Variable answer,dum
472        String weightStr,zStr
473       
474        weightStr = "gauss20wt"
475        zStr = "gauss20z"
476        nord = 20
477               
478        if (WaveExists($weightStr) == 0) // wave reference is not valid,
479                Make/D/N=20 $weightStr,$zStr
480                Wave w20 = $weightStr
481                Wave z20 = $zStr                // wave references to pass
482                Make20GaussPoints(w20,z20)     
483        else
484                if(exists(weightStr) > 1)
485                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
486                endif
487                Wave w20 = $weightStr
488                Wave z20 = $zStr                // create the wave references
489        endif
490
491        //limits of integration are input to function
492        va = loLim
493        vb = upLim
494        // Using 20 Gauss points                   
495        // remember to index from 0,size-1
496
497        summ = 0.0              // initialize integral
498        ii=0                    // loop counter
499        do
500                // calculate Gauss points on integration interval (q-value for evaluation)
501                zi = ( z20[ii]*(vb-va) + vb + va )/2.0
502                //calculate partial sum for the passed-in model function       
503                yyy = w20[ii] *  fcn(w,x,zi)                                           
504                summ += yyy             //add to the running total of the quadrature
505        ii+=1           
506        while (ii<nord)                         // end of loop over quadrature points
507   
508        // calculate value of integral to return
509        answer = (vb-va)/2.0*summ
510
511        Return (answer)
512End
513///////////////////////////////////////////////////////////////
514
515Function IntegrateFn76(fcn,loLim,upLim,w,x)                             
516        FUNCREF GenericQuadrature_proto fcn
517        Variable loLim,upLim    //limits of integration
518        Wave w                  //coefficients of function fcn(w,x)
519        Variable x      //x-value (q) for the calculation
520
521//**** The coefficient wave is passed into this function and straight through to the unsmeared model function
522
523// local variables
524        Variable nord,ii,va,vb,summ,yyy,zi
525        Variable answer,dum
526        String weightStr,zStr
527       
528        weightStr = "gauss76wt"
529        zStr = "gauss76z"
530        nord = 76
531               
532        if (WaveExists($weightStr) == 0) // wave reference is not valid,
533                Make/D/N=76 $weightStr,$zStr
534                Wave w76 = $weightStr
535                Wave z76 = $zStr                // wave references to pass
536                Make76GaussPoints(w76,z76)     
537        else
538                if(exists(weightStr) > 1)
539                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
540                endif
541                Wave w76 = $weightStr
542                Wave z76 = $zStr                // create the wave references
543        endif
544
545        //limits of integration are input to function
546        va = loLim
547        vb = upLim
548        // Using 76 Gauss points                   
549        // remember to index from 0,size-1
550
551        summ = 0.0              // initialize integral
552        ii=0                    // loop counter
553        do
554                // calculate Gauss points on integration interval (q-value for evaluation)
555                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
556                //calculate partial sum for the passed-in model function       
557                yyy = w76[ii] *  fcn(w,x,zi)                                           
558                summ += yyy             //add to the running total of the quadrature
559     ii+=1     
560        while (ii<nord)                         // end of loop over quadrature points
561   
562        // calculate value of integral to return
563        answer = (vb-va)/2.0*summ
564
565        Return (answer)
566End
567///////////////////////////////////////////////////////////////
568
569//////Resolution Smearing Utilities
570
571// To check for the existence of all waves needed for smearing
572// returns 1 if any waves are missing, 0 if all is OK
573Function ResolutionWavesMissing()
574
575        SVAR/Z sq = gSig_Q
576        SVAR/Z qb = gQ_bar
577        SVAR/Z sh = gShadow
578        SVAR/Z gQ = gQVals
579       
580        Variable err=0
581        if(!SVAR_Exists(sq) || !SVAR_Exists(qb) || !SVAR_Exists(sh) || !SVAR_Exists(gQ))
582                DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with LoadOneDData macro"
583                return(1)
584        endif
585       
586        if(WaveExists($sq) == 0)        //wave ref does not exist
587                err=1
588        endif
589        if(WaveExists($qb) == 0)        //wave ref does not exist
590                err=1
591        endif
592        if(WaveExists($sh) == 0)        //wave ref does not exist
593                err=1
594        endif
595        if(WaveExists($gQ) == 0)        //wave ref does not exist
596                err=1
597        endif
598       
599        if(err)
600                DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with 'Load SANS or USANS Data' macro"
601        endif
602       
603        return(err)
604end
605
606//////Resolution Smearing Utilities
607
608// To check for the existence of all waves needed for smearing
609// returns 1 if any waves are missing, 0 if all is OK
610// str passed in is the data folder containing the data
611//
612// 19 JUN07 using new data folder structure for loading
613// and resolution matrix
614Function ResolutionWavesMissingDF(str)
615        String str
616       
617        String DF="root:"+str+":"
618
619        WAVE/Z res = $(DF+str+"_res")
620       
621        if(!WaveExists(res))
622                DoAlert 0,"The resolution matrix is missing. Re-load experimental data with the 'Load SANS or USANS Data' macro"
623                return(1)
624        endif
625       
626        return(0)
627end
628
629///////////////////////////////////////////////////////////////
630
631// "backwards" wrapped to reduce redundant code
632// there are only 4 choices of N (5,10,20,76) for smearing
633Function Smear_Model_N(fcn,w,x,resW,wi,zi,nord)
634        FUNCREF SANSModelAAO_proto fcn
635        Wave w                  //coefficients of function fcn(w,x)
636        Variable x      //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE
637        Wave resW               // Nx4 or NxN matrix of resolution
638        Wave wi         //weight wave
639        Wave zi         //abscissa wave
640        Variable nord           //order of integration
641
642// local variables
643        Variable ii,va,vb
644        Variable answer,i_shad,i_qbar,i_sigq
645
646        // current x point is the q-value for evaluation
647        //
648        // * for the input x, the resolution function waves are interpolated to get the correct values for
649        //  sigq, qbar and shad - since the model x-spacing may not be the same as
650        // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is
651        // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different
652        // from experimental data.
653        // **note** if the (x) passed in is the experimental q-values, these values are
654        // returned from the interpolation (as expected)
655
656        Make/O/D/N=(DimSize(resW, 0)) sigQ,qbar,shad,qvals
657        sigq = resW[p][0]               //std dev of resolution fn
658        qbar = resW[p][1]               //mean q-value
659        shad = resW[p][2]               //beamstop shadow factor
660        qvals = resW[p][3]      //q-values where R(q) is known
661
662        i_shad = interp(x,qvals,shad)
663        i_qbar = interp(x,qvals,qbar)
664        i_sigq = interp(x,qvals,sigq)
665                       
666// set up the integration
667// number of Gauss Quadrature points
668               
669        if (isSANSResolution(i_sigq))
670               
671                // end points of integration
672                // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
673                // +/- 3 sigq catches 99.73% of distrubution
674                // change limits (and spacing of zi) at each evaluation based on R()
675                //integration from va to vb
676       
677                va = -3*i_sigq + i_qbar
678                if (va<0)
679                        va=0            //to avoid numerical error when  va<0 (-ve q-value)
680//                      Print "truncated Gaussian at nominal q = ",x
681                endif
682                vb = 3*i_sigq + i_qbar
683               
684                // Using 20 Gauss points                   
685                ii=0                    // loop counter
686                // do the calculation as a single pass w/AAO function
687                Make/O/D/N=(nord) Resoln,yyy,xGauss
688                do
689                        // calculate Gauss points on integration interval (q-value for evaluation)
690                        xGauss[ii] = ( zi[ii]*(vb-va) + vb + va )/2.0
691                        // calculate resolution function at input q-value (use the interpolated values and zi)
692                        Resoln[ii] = i_shad/sqrt(2*pi*i_sigq*i_sigq)
693                        Resoln[ii] *= exp((-1*(xGauss[ii] - i_qbar)^2)/(2*i_sigq*i_sigq))
694                        ii+=1           
695                while (ii<nord)                         // end of loop over quadrature points
696               
697                fcn(w,yyy,xGauss)               //yyy is the return value as a wave
698               
699                yyy *= wi *Resoln               //multiply function by resolution and weights
700                // calculate value of integral to return
701                answer = (vb-va)/2.0*sum(yyy)
702                // all scaling, background addition... etc. is done in the model calculation
703       
704        else
705                //smear with the USANS routine
706                // Make global string and local variables
707                // now data folder aware, necessary for GlobalFit = FULL path to wave   
708                String/G gTrap_coefStr = GetWavesDataFolder(w, 2 )     
709                Variable maxiter=20, tol=1e-4
710               
711                // set up limits for the integration
712                va=0
713                vb=abs(i_sigq)
714               
715                Variable/G gEvalQval = x
716               
717                // call qtrap to do actual work
718                answer = qtrap_USANS(fcn,va,vb,tol,maxiter)
719                answer /= vb
720               
721        endif
722       
723        //killing these waves is cleaner, but MUCH SLOWER
724//      Killwaves/Z Resoln,yyy,xGauss
725//      Killwaves/Z sigQ,qbar,shad,qvals
726        Return (answer)
727       
728End
729
730//resolution smearing, using only 5 Gauss points
731//
732//
733Function Smear_Model_5(fcn,w,x,answer,resW)                             
734        FUNCREF SANSModelAAO_proto fcn
735        Wave w                  //coefficients of function fcn(w,x)
736        Wave x  //x-value (q) for the calculation
737        Wave answer // ywave for calculation result
738        Wave resW               // Nx4 or NxN matrix of resolution
739        NVAR useTrap = root:Packages:NIST:USANSUseTrap
740
741        String weightStr,zStr
742        Variable nord=5
743       
744        if (dimsize(resW,1) > 4 && useTrap !=1)
745                //USANS Weighting matrix is present.
746                fcn(w,answer,x)
747       
748                MatrixOP/O  answer = resW x answer
749                //Duplicate/O answer,tmpMat
750                //MatrixOP/O answer = resW x tmpMat
751                Return(0)
752        else
753                weightStr = "gauss20wt"
754                zStr = "gauss20z"
755       
756        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
757                if (WaveExists($weightStr) == 0) // wave reference is not valid,
758                        Make/D/N=(nord) $weightStr,$zStr
759                        Wave weightW = $weightStr
760                        Wave abscissW = $zStr           // wave references to pass
761                        Make20GaussPoints(weightW,abscissW)     
762                else
763                        if(exists(weightStr) > 1)
764                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
765                        endif
766                        Wave weightW = $weightStr
767                        Wave abscissW = $zStr           // create the wave references
768                endif
769       
770                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
771                Return (0)
772        endif
773       
774End
775
776//resolution smearing, using only 10 Gauss points
777//
778//
779Function Smear_Model_10(fcn,w,x,answer,resW)                           
780        FUNCREF SANSModelAAO_proto fcn
781        Wave w                  //coefficients of function fcn(w,x)
782        Wave x  //x-value (q) for the calculation
783        Wave answer // ywave for calculation result
784        Wave resW               // Nx4 or NxN matrix of resolution
785
786        String weightStr,zStr
787        Variable nord=10
788       
789        if (dimsize(resW,1) > 4)
790                //USANS Weighting matrix is present.
791                fcn(w,answer,x)
792       
793                MatrixOP/O  answer = resW x answer
794                //Duplicate/O answer,tmpMat
795                //MatrixOP/O answer = resW x tmpMat
796                Return(0)
797        else
798                weightStr = "gauss20wt"
799                zStr = "gauss20z"
800       
801        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
802                if (WaveExists($weightStr) == 0) // wave reference is not valid,
803                        Make/D/N=(nord) $weightStr,$zStr
804                        Wave weightW = $weightStr
805                        Wave abscissW = $zStr           // wave references to pass
806                        Make20GaussPoints(weightW,abscissW)     
807                else
808                        if(exists(weightStr) > 1)
809                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
810                        endif
811                        Wave weightW = $weightStr
812                        Wave abscissW = $zStr           // create the wave references
813                endif
814       
815                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
816                Return (0)
817        endif
818       
819End
820
821//
822//Smear_Model_20(SphereForm,s.coefW,s.yW,s.xW,s.resW)
823//
824//      Wave sigq               //std dev of resolution fn
825//      Wave qbar               //mean q-value
826//      Wave shad               //beamstop shadow factor
827//      Wave qvals      //q-values where R(q) is known
828//
829Function Smear_Model_20(fcn,w,x,answer,resW)                           
830        FUNCREF SANSModelAAO_proto fcn
831        Wave w                  //coefficients of function fcn(w,x)
832        Wave x  //x-value (q) for the calculation
833        Wave answer // ywave for calculation result
834        Wave resW               // Nx4 or NxN matrix of resolution
835        NVAR useTrap = root:Packages:NIST:USANSUseTrap
836
837        String weightStr,zStr
838        Variable nord=20
839       
840        if (dimsize(resW,1) > 4 && useTrap != 1)
841                //USANS Weighting matrix is present.
842                fcn(w,answer,x)
843       
844                MatrixOP/O  answer = resW x answer
845                //Duplicate/O answer,tmpMat
846                //MatrixOP/O answer = resW x tmpMat
847                Return(0)
848        else
849                weightStr = "gauss20wt"
850                zStr = "gauss20z"
851       
852        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
853                if (WaveExists($weightStr) == 0) // wave reference is not valid,
854                        Make/D/N=(nord) $weightStr,$zStr
855                        Wave weightW = $weightStr
856                        Wave abscissW = $zStr           // wave references to pass
857                        Make20GaussPoints(weightW,abscissW)     
858                else
859                        if(exists(weightStr) > 1)
860                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
861                        endif
862                        Wave weightW = $weightStr
863                        Wave abscissW = $zStr           // create the wave references
864                endif
865       
866                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
867                Return (0)
868        endif
869       
870End
871///////////////////////////////////////////////////////////////
872Function Smear_Model_76(fcn,w,x,answer,resW)                           
873        FUNCREF SANSModelAAO_proto fcn
874        Wave w                  //coefficients of function fcn(w,x)
875        Wave x  //x-value (q) for the calculation
876        Wave answer // ywave for calculation result
877        Wave resW               // Nx4 or NxN matrix of resolution
878        NVAR useTrap = root:Packages:NIST:USANSUseTrap
879
880        String weightStr,zStr
881        Variable nord=76
882       
883        if (dimsize(resW,1) > 4  && useTrap != 1)
884                //USANS Weighting matrix is present.
885                fcn(w,answer,x)
886       
887                MatrixOP/O  answer = resW x answer
888                //Duplicate/O answer,tmpMat
889                //MatrixOP/O answer = resW x tmpMat
890                Return(0)
891        else
892                weightStr = "gauss20wt"
893                zStr = "gauss20z"
894       
895        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
896                if (WaveExists($weightStr) == 0) // wave reference is not valid,
897                        Make/D/N=(nord) $weightStr,$zStr
898                        Wave weightW = $weightStr
899                        Wave abscissW = $zStr           // wave references to pass
900                        Make20GaussPoints(weightW,abscissW)     
901                else
902                        if(exists(weightStr) > 1)
903                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
904                        endif
905                        Wave weightW = $weightStr
906                        Wave abscissW = $zStr           // create the wave references
907                endif
908       
909                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
910                Return (0)
911        endif
912       
913End
914
915
916///////////////////////////////////////////////////////////////
917
918//typically, the first point (or any point) of sigQ is passed
919// if negative, it's USANS data...
920Function isSANSResolution(val)
921        Variable val
922        if(val>=0)
923                return(1)               //true, SANS data
924        else
925                return(0)               //false, USANS
926        endif
927End
928
929Function GenericQuadrature_proto(w,x,dum)
930        Wave w
931        Variable x,dum
932       
933        Print "in GenericQuadrature_proto function"
934        return(1)
935end
936 
937// prototype function for smearing routines, Smear_Model_N
938// and trapzd_USANS() and qtrap_USANS()
939// it intentionally does nothing
940Function SANSModel_proto(w,x)
941        Wave w
942        Variable x
943       
944        Print "in SANSModel_proto function"
945        return(1)
946end
947
948// prototype function for smearing routines, Smear_Model_N
949// and trapzd_USANS() and qtrap_USANS()
950// it intentionally does nothing
951Function SANSModelAAO_proto(w,yw,xw)
952        Wave w,yw,xw
953       
954        Print "in SANSModelAAO_proto function"
955        return(1)
956end
957
958// prototype function for fit wrapper
959// it intentionally does nothing
960Function SANSModelSTRUCT_proto(s)
961        Struct ResSmearAAOStruct &s     
962
963        Print "in SANSModelSTRUCT_proto function"
964        return(1)
965end
966
967//Numerical Recipes routine to calculate the nn(th) stage
968//refinement of a trapezoid integration
969//
970//must be called sequentially from nn=1...n from qtrap()
971// to cumulatively refine the integration value
972//
973// in the conversion:
974// -- s was replaced with sVal and declared global (rather than static)
975//  so that the nn-1 value would be available during the nn(th) call
976//
977// -- the specific coefficient wave for func() is passed in as a
978//  global string (then converted to a wave reference), since
979//  func() will eventually call sphereForm()
980//
981Function trapzd_USANS(fcn,aa,bb,nn)
982        FUNCREF SANSModelAAO_proto fcn
983        Variable aa,bb,nn
984       
985        Variable xx,tnm,summ,del
986        Variable it,jj,arg1,arg2
987        NVAR sVal=sVal          //calling function must initialize this global
988        NVAR qval = gEvalQval
989        SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
990        Wave cw=$cwStr                 
991        Variable temp=0
992        Make/D/O/N=2 tmp_xw,tmp_yw
993        if(nn==1)
994                tmp_xw[0] = sqrt(qval^2 + aa^2)
995                tmp_xw[1] = sqrt(qval^2 + bb^2)
996                fcn(cw,tmp_yw,tmp_xw)
997                temp = 0.5*(bb-aa)*(tmp_yw[0] + tmp_yw[1])
998                sval = temp
999                return(sVal)
1000        else
1001                it=1
1002                it= 2^(nn-2)  //done in NR with a bit shift <<=
1003                tnm = it
1004                del = (bb - aa)/tnm             //this is the spacing of points to add
1005                xx = aa+0.5*del
1006                summ=0
1007                for(jj=1;jj<=it;jj+=1)
1008                        tmp_xw = sqrt(qval^2 + xx^2)
1009                        fcn(cw,tmp_yw,tmp_xw)           //not the most efficient... but replaced by the matrix method
1010                        summ += tmp_yw[0]
1011                        xx += del
1012                endfor
1013                sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1014                return (sval)
1015        endif
1016        //KillWaves/Z tmp_xw,tmp_yw
1017End
1018
1019////Numerical Recipes routine to calculate the nn(th) stage
1020////refinement of a trapezoid integration
1021////
1022////must be called sequentially from nn=1...n from qtrap()
1023//// to cumulatively refine the integration value
1024////
1025//// in the conversion:
1026//// -- s was replaced with sVal and declared global (rather than static)
1027////  so that the nn-1 value would be available during the nn(th) call
1028////
1029//// -- the specific coefficient wave for func() is passed in as a
1030////  global string (then converted to a wave reference), since
1031////  func() will eventually call sphereForm()
1032////
1033//Function trapzd_USANS_point(fcn,aa,bb,nn)
1034//      FUNCREF SANSModel_proto fcn
1035//      Variable aa,bb,nn
1036//     
1037//      Variable xx,tnm,summ,del
1038//      Variable it,jj,arg1,arg2
1039//      NVAR sVal=sVal          //calling function must initialize this global
1040//      NVAR qval = gEvalQval
1041//      SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
1042//      Wave cw=$cwStr                 
1043//      Variable temp=0
1044//      if(nn==1)
1045//              arg1 = qval^2 + aa^2
1046//              arg2 = qval^2 + bb^2
1047//              temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2)))
1048//              sval = temp
1049//              return(sVal)
1050//      else
1051//              it=1
1052//              it= 2^(nn-2)  //done in NR with a bit shift <<=
1053//              tnm = it
1054//              del = (bb - aa)/tnm             //this is the spacing of points to add
1055//              xx = aa+0.5*del
1056//              summ=0
1057//              for(jj=1;jj<=it;jj+=1)
1058//                      arg1 = qval^2 + xx^2
1059//                      summ += fcn(cw,sqrt(arg1))
1060//                      xx += del
1061//              endfor
1062//              sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1063//              return (sval)
1064//      endif
1065//     
1066//End
1067
1068// Numerical Recipes routine to calculate the integral of a
1069// specified function, trapezoid rule is used to a user-specified
1070// level of refinement using sequential calls to trapzd()
1071//
1072// in NR, eps and maxIt were global, pass them in here...
1073// eps typically 1e-5
1074// maxit typically 20
1075Function qtrap_USANS(fcn,aa,bb,eps,maxIt)
1076        FUNCREF SANSModelAAO_proto fcn
1077        Variable aa,bb,eps,maxit
1078       
1079        Variable/G sVal=0               //create and initialize what trapzd will return
1080        Variable jj,ss,olds
1081       
1082        olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1083        for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1084                ss = trapzd_USANS(fcn,aa,bb,jj)
1085                if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1086                        return ss
1087                endif
1088                if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1089                        return ss
1090                endif
1091                olds = ss
1092        endfor
1093
1094        Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1095        return(ss)              //should never get here if function is well-behaved
1096       
1097End     
1098
1099//// Numerical Recipes routine to calculate the integral of a
1100//// specified function, trapezoid rule is used to a user-specified
1101//// level of refinement using sequential calls to trapzd()
1102////
1103//// in NR, eps and maxIt were global, pass them in here...
1104//// eps typically 1e-5
1105//// maxit typically 20
1106//Function qtrap_USANS_point(fcn,aa,bb,eps,maxIt)
1107//      FUNCREF SANSModel_proto fcn
1108//      Variable aa,bb,eps,maxit
1109//     
1110//      Variable/G sVal=0               //create and initialize what trapzd will return
1111//      Variable jj,ss,olds
1112//     
1113//      olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1114//      for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1115//              ss = trapzd_USANS_point(fcn,aa,bb,jj)
1116//              if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1117//                      return ss
1118//              endif
1119//              if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1120//                      return ss
1121//              endif
1122//              olds = ss
1123//      endfor
1124//
1125//      Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1126//      return(ss)              //should never get here if function is well-behaved
1127//     
1128//End   
1129
1130Proc RRW()
1131        ResetResolutionWaves()
1132End
1133
1134//utility procedures that are currently untied to any actions, although useful...
1135Proc ResetResolutionWaves(str)
1136        String Str
1137        Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","")
1138
1139
1140        Abort "This function is not data floder aware and does nothing..."
1141               
1142        String/G root:gQvals = str[0,strlen(str)-3]+"_q"
1143        String/G root:gSig_Q = str[0,strlen(str)-3]+"sq"
1144        String/G root:gQ_bar = str[0,strlen(str)-3]+"qb"
1145        String/G root:gShadow = str[0,strlen(str)-3]+"fs"
1146       
1147        //touch everything to make sure that the dependencies are
1148        //properly updated - especially the $gQvals reference in the
1149        // dependency assignment
1150        fKillDependencies("Smear*")
1151       
1152        //replace the q-values and intensity (so they're the right length)
1153        fResetSmearedModels("Smear*",root:gQvals)
1154       
1155        fRestoreDependencies("Smear*")
1156End
1157
1158// pass "*" as the matchString to do ALL dependent waves
1159// or "abc*" to get just the matching waves
1160//
1161Function fKillDependencies(matchStr)
1162        String matchStr
1163
1164        String str=WaveList(matchStr, ";", "" ),item,formula
1165        Variable ii
1166       
1167        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1168                item = StringFromList(ii, str ,";")
1169                formula = GetFormula($item)
1170                if(cmpstr("",formula)!=0)
1171                        Printf "wave %s had the formula %s removed\r",item,formula
1172                        Note $item, "FORMULA:"+formula
1173                        SetFormula $item, ""                    //clears the formula
1174                endif
1175        endfor
1176        return(0)
1177end
1178
1179// pass "*" as the matchString to do ALL dependent waves
1180// or "abc*" to get just the matching waves
1181//
1182Function fRestoreDependencies(matchStr)
1183        String matchStr
1184
1185        String str=WaveList(matchStr, ";", "" ),item,formula
1186        Variable ii
1187       
1188        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1189                item = StringFromList(ii, str ,";")
1190                formula = StringByKey("FORMULA", note($item),":",";")
1191                if(cmpstr("",formula)!=0)
1192                        Printf "wave %s had the formula %s restored\r",item,formula
1193                        Note/K $item
1194                        SetFormula $item, formula                       //restores the formula
1195                endif
1196        endfor
1197        return(0)
1198end
1199
1200Function fResetSmearedModels(matchStr,qStr)
1201        String matchStr,qStr
1202
1203        Duplicate/O $qStr root:smeared_qvals   
1204       
1205        String str=WaveList(matchStr, ";", "" ),item,formula
1206        Variable ii
1207       
1208        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1209                item = StringFromList(ii, str ,";")
1210                formula = StringByKey("FORMULA", note($item),":",";")
1211                if(cmpstr("",formula)!=0)
1212                        Printf "wave %s has been duplicated to gQvals\r",item
1213                        Duplicate/O $qStr $item
1214                        Note $item, "FORMULA:"+formula          //be sure to keep the formula note
1215                        Print "   and the formula is",formula
1216                endif
1217        endfor
1218        return(0)
1219end
Note: See TracBrowser for help on using the repository browser.