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

Last change on this file since 309 was 309, checked in by srkline, 15 years ago

Changes to how USANS resolution matrix is handled in fitting, both regular and global.
(1) - resolution matirx is re-calculated if the data is masked (with cursors)
(2) - if the matrix is appropriate, as stored in the wave note, it is not recalculated
(3) - global fit does the same, being sure to recalculate the matrix before the Global fit has started
(4) - the global fit report now has pararameter names rather than "coef_n". will need to be tweaked in future to accomodate different fit functions, and to improve the formatting.

File size: 35.9 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        NVAR dQv = root:Packages:NIST:USANS_dQv
643
644// local variables
645        Variable ii,va,vb
646        Variable answer,i_shad,i_qbar,i_sigq
647
648        // current x point is the q-value for evaluation
649        //
650        // * for the input x, the resolution function waves are interpolated to get the correct values for
651        //  sigq, qbar and shad - since the model x-spacing may not be the same as
652        // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is
653        // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different
654        // from experimental data.
655        // **note** if the (x) passed in is the experimental q-values, these values are
656        // returned from the interpolation (as expected)
657
658        Make/O/D/N=(DimSize(resW, 0)) sigQ,qbar,shad,qvals
659        sigq = resW[p][0]               //std dev of resolution fn
660        qbar = resW[p][1]               //mean q-value
661        shad = resW[p][2]               //beamstop shadow factor
662        qvals = resW[p][3]      //q-values where R(q) is known
663
664        i_shad = interp(x,qvals,shad)
665        i_qbar = interp(x,qvals,qbar)
666        i_sigq = interp(x,qvals,sigq)
667                       
668// set up the integration
669// number of Gauss Quadrature points
670               
671        if (isSANSResolution(i_sigq))
672               
673                // end points of integration
674                // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
675                // +/- 3 sigq catches 99.73% of distrubution
676                // change limits (and spacing of zi) at each evaluation based on R()
677                //integration from va to vb
678       
679                va = -3*i_sigq + i_qbar
680                if (va<0)
681                        va=0            //to avoid numerical error when  va<0 (-ve q-value)
682//                      Print "truncated Gaussian at nominal q = ",x
683                endif
684                vb = 3*i_sigq + i_qbar
685               
686                // Using 20 Gauss points                   
687                ii=0                    // loop counter
688                // do the calculation as a single pass w/AAO function
689                Make/O/D/N=(nord) Resoln,yyy,xGauss
690                do
691                        // calculate Gauss points on integration interval (q-value for evaluation)
692                        xGauss[ii] = ( zi[ii]*(vb-va) + vb + va )/2.0
693                        // calculate resolution function at input q-value (use the interpolated values and zi)
694                        Resoln[ii] = i_shad/sqrt(2*pi*i_sigq*i_sigq)
695                        Resoln[ii] *= exp((-1*(xGauss[ii] - i_qbar)^2)/(2*i_sigq*i_sigq))
696                        ii+=1           
697                while (ii<nord)                         // end of loop over quadrature points
698               
699                fcn(w,yyy,xGauss)               //yyy is the return value as a wave
700               
701                yyy *= wi *Resoln               //multiply function by resolution and weights
702                // calculate value of integral to return
703                answer = (vb-va)/2.0*sum(yyy)
704                // all scaling, background addition... etc. is done in the model calculation
705       
706        else
707                //smear with the USANS routine
708                // Make global string and local variables
709                // now data folder aware, necessary for GlobalFit = FULL path to wave   
710                String/G gTrap_coefStr = GetWavesDataFolder(w, 2 )     
711                Variable maxiter=20, tol=1e-4
712               
713                // set up limits for the integration
714                va=0
715                vb=abs(dQv)
716               
717                Variable/G gEvalQval = x
718               
719                // call qtrap to do actual work
720                answer = qtrap_USANS(fcn,va,vb,tol,maxiter)
721                answer /= vb
722               
723        endif
724       
725        //killing these waves is cleaner, but MUCH SLOWER
726//      Killwaves/Z Resoln,yyy,xGauss
727//      Killwaves/Z sigQ,qbar,shad,qvals
728        Return (answer)
729       
730End
731
732//resolution smearing, using only 5 Gauss points
733//
734//
735Function Smear_Model_5(fcn,w,x,answer,resW)                             
736        FUNCREF SANSModelAAO_proto fcn
737        Wave w                  //coefficients of function fcn(w,x)
738        Wave x  //x-value (q) for the calculation
739        Wave answer // ywave for calculation result
740        Wave resW               // Nx4 or NxN matrix of resolution
741        NVAR useTrap = root:Packages:NIST:USANSUseTrap
742
743        String weightStr,zStr
744        Variable nord=5
745       
746        if (dimsize(resW,1) > 4 && useTrap !=1)
747                if(dimsize(resW,1) != dimsize(answer,0) )
748                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
749                endif
750                //USANS Weighting matrix is present.
751                fcn(w,answer,x)
752       
753                MatrixOP/O  answer = resW x answer
754                //Duplicate/O answer,tmpMat
755                //MatrixOP/O answer = resW x tmpMat
756                Return(0)
757        else
758                weightStr = "gauss20wt"
759                zStr = "gauss20z"
760       
761        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
762                if (WaveExists($weightStr) == 0) // wave reference is not valid,
763                        Make/D/N=(nord) $weightStr,$zStr
764                        Wave weightW = $weightStr
765                        Wave abscissW = $zStr           // wave references to pass
766                        Make20GaussPoints(weightW,abscissW)     
767                else
768                        if(exists(weightStr) > 1)
769                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
770                        endif
771                        Wave weightW = $weightStr
772                        Wave abscissW = $zStr           // create the wave references
773                endif
774       
775                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
776                Return (0)
777        endif
778       
779End
780
781//resolution smearing, using only 10 Gauss points
782//
783//
784Function Smear_Model_10(fcn,w,x,answer,resW)                           
785        FUNCREF SANSModelAAO_proto fcn
786        Wave w                  //coefficients of function fcn(w,x)
787        Wave x  //x-value (q) for the calculation
788        Wave answer // ywave for calculation result
789        Wave resW               // Nx4 or NxN matrix of resolution
790
791        String weightStr,zStr
792        Variable nord=10
793       
794        if (dimsize(resW,1) > 4)
795                if(dimsize(resW,1) != dimsize(answer,0) )
796                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
797                endif
798                //USANS Weighting matrix is present.
799                fcn(w,answer,x)
800       
801                MatrixOP/O  answer = resW x answer
802                //Duplicate/O answer,tmpMat
803                //MatrixOP/O answer = resW x tmpMat
804                Return(0)
805        else
806                weightStr = "gauss20wt"
807                zStr = "gauss20z"
808       
809        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
810                if (WaveExists($weightStr) == 0) // wave reference is not valid,
811                        Make/D/N=(nord) $weightStr,$zStr
812                        Wave weightW = $weightStr
813                        Wave abscissW = $zStr           // wave references to pass
814                        Make20GaussPoints(weightW,abscissW)     
815                else
816                        if(exists(weightStr) > 1)
817                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
818                        endif
819                        Wave weightW = $weightStr
820                        Wave abscissW = $zStr           // create the wave references
821                endif
822       
823                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
824                Return (0)
825        endif
826       
827End
828
829//
830//Smear_Model_20(SphereForm,s.coefW,s.yW,s.xW,s.resW)
831//
832//      Wave sigq               //std dev of resolution fn
833//      Wave qbar               //mean q-value
834//      Wave shad               //beamstop shadow factor
835//      Wave qvals      //q-values where R(q) is known
836//
837Function Smear_Model_20(fcn,w,x,answer,resW)                           
838        FUNCREF SANSModelAAO_proto fcn
839        Wave w                  //coefficients of function fcn(w,x)
840        Wave x  //x-value (q) for the calculation
841        Wave answer // ywave for calculation result
842        Wave resW               // Nx4 or NxN matrix of resolution
843        NVAR useTrap = root:Packages:NIST:USANSUseTrap
844
845        String weightStr,zStr
846        Variable nord=20
847       
848        if (dimsize(resW,1) > 4 && useTrap != 1)
849                if(dimsize(resW,1) != dimsize(answer,0) )
850                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
851                endif
852                //USANS Weighting matrix is present.
853                fcn(w,answer,x)
854       
855                MatrixOP/O  answer = resW x answer
856                //Duplicate/O answer,tmpMat
857                //MatrixOP/O answer = resW x tmpMat
858                Return(0)
859        else
860                weightStr = "gauss20wt"
861                zStr = "gauss20z"
862       
863        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
864                if (WaveExists($weightStr) == 0) // wave reference is not valid,
865                        Make/D/N=(nord) $weightStr,$zStr
866                        Wave weightW = $weightStr
867                        Wave abscissW = $zStr           // wave references to pass
868                        Make20GaussPoints(weightW,abscissW)     
869                else
870                        if(exists(weightStr) > 1)
871                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
872                        endif
873                        Wave weightW = $weightStr
874                        Wave abscissW = $zStr           // create the wave references
875                endif
876       
877                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
878                Return (0)
879        endif
880       
881End
882///////////////////////////////////////////////////////////////
883Function Smear_Model_76(fcn,w,x,answer,resW)                           
884        FUNCREF SANSModelAAO_proto fcn
885        Wave w                  //coefficients of function fcn(w,x)
886        Wave x  //x-value (q) for the calculation
887        Wave answer // ywave for calculation result
888        Wave resW               // Nx4 or NxN matrix of resolution
889        NVAR useTrap = root:Packages:NIST:USANSUseTrap
890
891        String weightStr,zStr
892        Variable nord=76
893       
894        if (dimsize(resW,1) > 4  && useTrap != 1)
895                if(dimsize(resW,1) != dimsize(answer,0) )
896                        Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0))
897                endif
898                //USANS Weighting matrix is present.
899                fcn(w,answer,x)
900       
901                MatrixOP/O  answer = resW x answer
902                //Duplicate/O answer,tmpMat
903                //MatrixOP/O answer = resW x tmpMat
904                Return(0)
905        else
906                weightStr = "gauss20wt"
907                zStr = "gauss20z"
908       
909        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
910                if (WaveExists($weightStr) == 0) // wave reference is not valid,
911                        Make/D/N=(nord) $weightStr,$zStr
912                        Wave weightW = $weightStr
913                        Wave abscissW = $zStr           // wave references to pass
914                        Make20GaussPoints(weightW,abscissW)     
915                else
916                        if(exists(weightStr) > 1)
917                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
918                        endif
919                        Wave weightW = $weightStr
920                        Wave abscissW = $zStr           // create the wave references
921                endif
922       
923                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
924                Return (0)
925        endif
926       
927End
928
929
930///////////////////////////////////////////////////////////////
931
932//typically, the first point (or any point) of sigQ is passed
933// if negative, it's USANS data...
934Function isSANSResolution(val)
935        Variable val
936        if(val>=0)
937                return(1)               //true, SANS data
938        else
939                return(0)               //false, USANS
940        endif
941End
942
943Function GenericQuadrature_proto(w,x,dum)
944        Wave w
945        Variable x,dum
946       
947        Print "in GenericQuadrature_proto function"
948        return(1)
949end
950 
951// prototype function for smearing routines, Smear_Model_N
952// and trapzd_USANS() and qtrap_USANS()
953// it intentionally does nothing
954Function SANSModel_proto(w,x)
955        Wave w
956        Variable x
957       
958        Print "in SANSModel_proto function"
959        return(1)
960end
961
962// prototype function for smearing routines, Smear_Model_N
963// and trapzd_USANS() and qtrap_USANS()
964// it intentionally does nothing
965Function SANSModelAAO_proto(w,yw,xw)
966        Wave w,yw,xw
967       
968        Print "in SANSModelAAO_proto function"
969        return(1)
970end
971
972// prototype function for fit wrapper
973// it intentionally does nothing
974Function SANSModelSTRUCT_proto(s)
975        Struct ResSmearAAOStruct &s     
976
977        Print "in SANSModelSTRUCT_proto function"
978        return(1)
979end
980
981//Numerical Recipes routine to calculate the nn(th) stage
982//refinement of a trapezoid integration
983//
984//must be called sequentially from nn=1...n from qtrap()
985// to cumulatively refine the integration value
986//
987// in the conversion:
988// -- s was replaced with sVal and declared global (rather than static)
989//  so that the nn-1 value would be available during the nn(th) call
990//
991// -- the specific coefficient wave for func() is passed in as a
992//  global string (then converted to a wave reference), since
993//  func() will eventually call sphereForm()
994//
995Function trapzd_USANS(fcn,aa,bb,nn)
996        FUNCREF SANSModelAAO_proto fcn
997        Variable aa,bb,nn
998       
999        Variable xx,tnm,summ,del
1000        Variable it,jj,arg1,arg2
1001        NVAR sVal=sVal          //calling function must initialize this global
1002        NVAR qval = gEvalQval
1003        SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
1004        Wave cw=$cwStr                 
1005        Variable temp=0
1006        Make/D/O/N=2 tmp_xw,tmp_yw
1007        if(nn==1)
1008                tmp_xw[0] = sqrt(qval^2 + aa^2)
1009                tmp_xw[1] = sqrt(qval^2 + bb^2)
1010                fcn(cw,tmp_yw,tmp_xw)
1011                temp = 0.5*(bb-aa)*(tmp_yw[0] + tmp_yw[1])
1012                sval = temp
1013                return(sVal)
1014        else
1015                it=1
1016                it= 2^(nn-2)  //done in NR with a bit shift <<=
1017                tnm = it
1018                del = (bb - aa)/tnm             //this is the spacing of points to add
1019                xx = aa+0.5*del
1020                summ=0
1021                for(jj=1;jj<=it;jj+=1)
1022                        tmp_xw = sqrt(qval^2 + xx^2)
1023                        fcn(cw,tmp_yw,tmp_xw)           //not the most efficient... but replaced by the matrix method
1024                        summ += tmp_yw[0]
1025                        xx += del
1026                endfor
1027                sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1028                return (sval)
1029        endif
1030        //KillWaves/Z tmp_xw,tmp_yw
1031End
1032
1033////Numerical Recipes routine to calculate the nn(th) stage
1034////refinement of a trapezoid integration
1035////
1036////must be called sequentially from nn=1...n from qtrap()
1037//// to cumulatively refine the integration value
1038////
1039//// in the conversion:
1040//// -- s was replaced with sVal and declared global (rather than static)
1041////  so that the nn-1 value would be available during the nn(th) call
1042////
1043//// -- the specific coefficient wave for func() is passed in as a
1044////  global string (then converted to a wave reference), since
1045////  func() will eventually call sphereForm()
1046////
1047//Function trapzd_USANS_point(fcn,aa,bb,nn)
1048//      FUNCREF SANSModel_proto fcn
1049//      Variable aa,bb,nn
1050//     
1051//      Variable xx,tnm,summ,del
1052//      Variable it,jj,arg1,arg2
1053//      NVAR sVal=sVal          //calling function must initialize this global
1054//      NVAR qval = gEvalQval
1055//      SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
1056//      Wave cw=$cwStr                 
1057//      Variable temp=0
1058//      if(nn==1)
1059//              arg1 = qval^2 + aa^2
1060//              arg2 = qval^2 + bb^2
1061//              temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2)))
1062//              sval = temp
1063//              return(sVal)
1064//      else
1065//              it=1
1066//              it= 2^(nn-2)  //done in NR with a bit shift <<=
1067//              tnm = it
1068//              del = (bb - aa)/tnm             //this is the spacing of points to add
1069//              xx = aa+0.5*del
1070//              summ=0
1071//              for(jj=1;jj<=it;jj+=1)
1072//                      arg1 = qval^2 + xx^2
1073//                      summ += fcn(cw,sqrt(arg1))
1074//                      xx += del
1075//              endfor
1076//              sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1077//              return (sval)
1078//      endif
1079//     
1080//End
1081
1082// Numerical Recipes routine to calculate the integral of a
1083// specified function, trapezoid rule is used to a user-specified
1084// level of refinement using sequential calls to trapzd()
1085//
1086// in NR, eps and maxIt were global, pass them in here...
1087// eps typically 1e-5
1088// maxit typically 20
1089Function qtrap_USANS(fcn,aa,bb,eps,maxIt)
1090        FUNCREF SANSModelAAO_proto fcn
1091        Variable aa,bb,eps,maxit
1092       
1093        Variable/G sVal=0               //create and initialize what trapzd will return
1094        Variable jj,ss,olds
1095       
1096        olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1097        for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1098                ss = trapzd_USANS(fcn,aa,bb,jj)
1099                if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1100                        return ss
1101                endif
1102                if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1103                        return ss
1104                endif
1105                olds = ss
1106        endfor
1107
1108        Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1109        return(ss)              //should never get here if function is well-behaved
1110       
1111End     
1112
1113//// Numerical Recipes routine to calculate the integral of a
1114//// specified function, trapezoid rule is used to a user-specified
1115//// level of refinement using sequential calls to trapzd()
1116////
1117//// in NR, eps and maxIt were global, pass them in here...
1118//// eps typically 1e-5
1119//// maxit typically 20
1120//Function qtrap_USANS_point(fcn,aa,bb,eps,maxIt)
1121//      FUNCREF SANSModel_proto fcn
1122//      Variable aa,bb,eps,maxit
1123//     
1124//      Variable/G sVal=0               //create and initialize what trapzd will return
1125//      Variable jj,ss,olds
1126//     
1127//      olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1128//      for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1129//              ss = trapzd_USANS_point(fcn,aa,bb,jj)
1130//              if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1131//                      return ss
1132//              endif
1133//              if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1134//                      return ss
1135//              endif
1136//              olds = ss
1137//      endfor
1138//
1139//      Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1140//      return(ss)              //should never get here if function is well-behaved
1141//     
1142//End   
1143
1144Proc RRW()
1145        ResetResolutionWaves()
1146End
1147
1148//utility procedures that are currently untied to any actions, although useful...
1149Proc ResetResolutionWaves(str)
1150        String Str
1151        Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","")
1152
1153
1154        Abort "This function is not data floder aware and does nothing..."
1155               
1156        String/G root:gQvals = str[0,strlen(str)-3]+"_q"
1157        String/G root:gSig_Q = str[0,strlen(str)-3]+"sq"
1158        String/G root:gQ_bar = str[0,strlen(str)-3]+"qb"
1159        String/G root:gShadow = str[0,strlen(str)-3]+"fs"
1160       
1161        //touch everything to make sure that the dependencies are
1162        //properly updated - especially the $gQvals reference in the
1163        // dependency assignment
1164        fKillDependencies("Smear*")
1165       
1166        //replace the q-values and intensity (so they're the right length)
1167        fResetSmearedModels("Smear*",root:gQvals)
1168       
1169        fRestoreDependencies("Smear*")
1170End
1171
1172// pass "*" as the matchString to do ALL dependent waves
1173// or "abc*" to get just the matching waves
1174//
1175Function fKillDependencies(matchStr)
1176        String matchStr
1177
1178        String str=WaveList(matchStr, ";", "" ),item,formula
1179        Variable ii
1180       
1181        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1182                item = StringFromList(ii, str ,";")
1183                formula = GetFormula($item)
1184                if(cmpstr("",formula)!=0)
1185                        Printf "wave %s had the formula %s removed\r",item,formula
1186                        Note $item, "FORMULA:"+formula
1187                        SetFormula $item, ""                    //clears the formula
1188                endif
1189        endfor
1190        return(0)
1191end
1192
1193// pass "*" as the matchString to do ALL dependent waves
1194// or "abc*" to get just the matching waves
1195//
1196Function fRestoreDependencies(matchStr)
1197        String matchStr
1198
1199        String str=WaveList(matchStr, ";", "" ),item,formula
1200        Variable ii
1201       
1202        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1203                item = StringFromList(ii, str ,";")
1204                formula = StringByKey("FORMULA", note($item),":",";")
1205                if(cmpstr("",formula)!=0)
1206                        Printf "wave %s had the formula %s restored\r",item,formula
1207                        Note/K $item
1208                        SetFormula $item, formula                       //restores the formula
1209                endif
1210        endfor
1211        return(0)
1212end
1213
1214Function fResetSmearedModels(matchStr,qStr)
1215        String matchStr,qStr
1216
1217        Duplicate/O $qStr root:smeared_qvals   
1218       
1219        String str=WaveList(matchStr, ";", "" ),item,formula
1220        Variable ii
1221       
1222        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1223                item = StringFromList(ii, str ,";")
1224                formula = StringByKey("FORMULA", note($item),":",";")
1225                if(cmpstr("",formula)!=0)
1226                        Printf "wave %s has been duplicated to gQvals\r",item
1227                        Duplicate/O $qStr $item
1228                        Note $item, "FORMULA:"+formula          //be sure to keep the formula note
1229                        Print "   and the formula is",formula
1230                endif
1231        endfor
1232        return(0)
1233end
Note: See TracBrowser for help on using the repository browser.