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

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

Updated version numbers in files (4.0) and corrected a typo in SA_includes_v301

File size: 35.0 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
740        String weightStr,zStr
741        Variable nord=5
742       
743        if (dimsize(resW,1) > 4)
744                //USANS Weighting matrix is present.
745                fcn(w,answer,x)
746       
747                MatrixOP/O  answer = resW x answer
748                //Duplicate/O answer,tmpMat
749                //MatrixOP/O answer = resW x tmpMat
750                Return(0)
751        else
752                weightStr = "gauss20wt"
753                zStr = "gauss20z"
754       
755        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
756                if (WaveExists($weightStr) == 0) // wave reference is not valid,
757                        Make/D/N=(nord) $weightStr,$zStr
758                        Wave weightW = $weightStr
759                        Wave abscissW = $zStr           // wave references to pass
760                        Make20GaussPoints(weightW,abscissW)     
761                else
762                        if(exists(weightStr) > 1)
763                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
764                        endif
765                        Wave weightW = $weightStr
766                        Wave abscissW = $zStr           // create the wave references
767                endif
768       
769                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
770                Return (0)
771        endif
772       
773End
774
775//resolution smearing, using only 10 Gauss points
776//
777//
778Function Smear_Model_10(fcn,w,x,answer,resW)                           
779        FUNCREF SANSModelAAO_proto fcn
780        Wave w                  //coefficients of function fcn(w,x)
781        Wave x  //x-value (q) for the calculation
782        Wave answer // ywave for calculation result
783        Wave resW               // Nx4 or NxN matrix of resolution
784
785        String weightStr,zStr
786        Variable nord=10
787       
788        if (dimsize(resW,1) > 4)
789                //USANS Weighting matrix is present.
790                fcn(w,answer,x)
791       
792                MatrixOP/O  answer = resW x answer
793                //Duplicate/O answer,tmpMat
794                //MatrixOP/O answer = resW x tmpMat
795                Return(0)
796        else
797                weightStr = "gauss20wt"
798                zStr = "gauss20z"
799       
800        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
801                if (WaveExists($weightStr) == 0) // wave reference is not valid,
802                        Make/D/N=(nord) $weightStr,$zStr
803                        Wave weightW = $weightStr
804                        Wave abscissW = $zStr           // wave references to pass
805                        Make20GaussPoints(weightW,abscissW)     
806                else
807                        if(exists(weightStr) > 1)
808                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
809                        endif
810                        Wave weightW = $weightStr
811                        Wave abscissW = $zStr           // create the wave references
812                endif
813       
814                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
815                Return (0)
816        endif
817       
818End
819
820//
821//Smear_Model_20(SphereForm,s.coefW,s.yW,s.xW,s.resW)
822//
823//      Wave sigq               //std dev of resolution fn
824//      Wave qbar               //mean q-value
825//      Wave shad               //beamstop shadow factor
826//      Wave qvals      //q-values where R(q) is known
827//
828Function Smear_Model_20(fcn,w,x,answer,resW)                           
829        FUNCREF SANSModelAAO_proto fcn
830        Wave w                  //coefficients of function fcn(w,x)
831        Wave x  //x-value (q) for the calculation
832        Wave answer // ywave for calculation result
833        Wave resW               // Nx4 or NxN matrix of resolution
834
835        String weightStr,zStr
836        Variable nord=20
837       
838        if (dimsize(resW,1) > 4)
839                //USANS Weighting matrix is present.
840                fcn(w,answer,x)
841       
842                MatrixOP/O  answer = resW x answer
843                //Duplicate/O answer,tmpMat
844                //MatrixOP/O answer = resW x tmpMat
845                Return(0)
846        else
847                weightStr = "gauss20wt"
848                zStr = "gauss20z"
849       
850        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
851                if (WaveExists($weightStr) == 0) // wave reference is not valid,
852                        Make/D/N=(nord) $weightStr,$zStr
853                        Wave weightW = $weightStr
854                        Wave abscissW = $zStr           // wave references to pass
855                        Make20GaussPoints(weightW,abscissW)     
856                else
857                        if(exists(weightStr) > 1)
858                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
859                        endif
860                        Wave weightW = $weightStr
861                        Wave abscissW = $zStr           // create the wave references
862                endif
863       
864                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
865                Return (0)
866        endif
867       
868End
869///////////////////////////////////////////////////////////////
870Function Smear_Model_76(fcn,w,x,answer,resW)                           
871        FUNCREF SANSModelAAO_proto fcn
872        Wave w                  //coefficients of function fcn(w,x)
873        Wave x  //x-value (q) for the calculation
874        Wave answer // ywave for calculation result
875        Wave resW               // Nx4 or NxN matrix of resolution
876
877        String weightStr,zStr
878        Variable nord=76
879       
880        if (dimsize(resW,1) > 4)
881                //USANS Weighting matrix is present.
882                fcn(w,answer,x)
883       
884                MatrixOP/O  answer = resW x answer
885                //Duplicate/O answer,tmpMat
886                //MatrixOP/O answer = resW x tmpMat
887                Return(0)
888        else
889                weightStr = "gauss20wt"
890                zStr = "gauss20z"
891       
892        //      if wt,z waves don't exist, create them (only check for weight, should really check for both)
893                if (WaveExists($weightStr) == 0) // wave reference is not valid,
894                        Make/D/N=(nord) $weightStr,$zStr
895                        Wave weightW = $weightStr
896                        Wave abscissW = $zStr           // wave references to pass
897                        Make20GaussPoints(weightW,abscissW)     
898                else
899                        if(exists(weightStr) > 1)
900                                 Abort "wave name is already in use"            //executed only if name is in use elsewhere
901                        endif
902                        Wave weightW = $weightStr
903                        Wave abscissW = $zStr           // create the wave references
904                endif
905       
906                answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord)
907                Return (0)
908        endif
909       
910End
911
912
913///////////////////////////////////////////////////////////////
914
915//typically, the first point (or any point) of sigQ is passed
916// if negative, it's USANS data...
917Function isSANSResolution(val)
918        Variable val
919        if(val>=0)
920                return(1)               //true, SANS data
921        else
922                return(0)               //false, USANS
923        endif
924End
925
926Function GenericQuadrature_proto(w,x,dum)
927        Wave w
928        Variable x,dum
929       
930        Print "in GenericQuadrature_proto function"
931        return(1)
932end
933 
934// prototype function for smearing routines, Smear_Model_N
935// and trapzd_USANS() and qtrap_USANS()
936// it intentionally does nothing
937Function SANSModel_proto(w,x)
938        Wave w
939        Variable x
940       
941        Print "in SANSModel_proto function"
942        return(1)
943end
944
945// prototype function for smearing routines, Smear_Model_N
946// and trapzd_USANS() and qtrap_USANS()
947// it intentionally does nothing
948Function SANSModelAAO_proto(w,yw,xw)
949        Wave w,yw,xw
950       
951        Print "in SANSModelAAO_proto function"
952        return(1)
953end
954
955// prototype function for fit wrapper
956// it intentionally does nothing
957Function SANSModelSTRUCT_proto(s)
958        Struct ResSmearAAOStruct &s     
959
960        Print "in SANSModelSTRUCT_proto function"
961        return(1)
962end
963
964//Numerical Recipes routine to calculate the nn(th) stage
965//refinement of a trapezoid integration
966//
967//must be called sequentially from nn=1...n from qtrap()
968// to cumulatively refine the integration value
969//
970// in the conversion:
971// -- s was replaced with sVal and declared global (rather than static)
972//  so that the nn-1 value would be available during the nn(th) call
973//
974// -- the specific coefficient wave for func() is passed in as a
975//  global string (then converted to a wave reference), since
976//  func() will eventually call sphereForm()
977//
978Function trapzd_USANS(fcn,aa,bb,nn)
979        FUNCREF SANSModelAAO_proto fcn
980        Variable aa,bb,nn
981       
982        Variable xx,tnm,summ,del
983        Variable it,jj,arg1,arg2
984        NVAR sVal=sVal          //calling function must initialize this global
985        NVAR qval = gEvalQval
986        SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
987        Wave cw=$cwStr                 
988        Variable temp=0
989        Make/D/O/N=2 tmp_xw,tmp_yw
990        if(nn==1)
991                tmp_xw[0] = sqrt(qval^2 + aa^2)
992                tmp_xw[1] = sqrt(qval^2 + bb^2)
993                fcn(cw,tmp_yw,tmp_xw)
994                temp = 0.5*(bb-aa)*(tmp_yw[0] + tmp_yw[1])
995                sval = temp
996                return(sVal)
997        else
998                it=1
999                it= 2^(nn-2)  //done in NR with a bit shift <<=
1000                tnm = it
1001                del = (bb - aa)/tnm             //this is the spacing of points to add
1002                xx = aa+0.5*del
1003                summ=0
1004                for(jj=1;jj<=it;jj+=1)
1005                        tmp_xw = sqrt(qval^2 + xx^2)
1006                        fcn(cw,tmp_yw,tmp_xw)           //not the most efficient... but replaced by the matrix method
1007                        summ += tmp_yw[0]
1008                        xx += del
1009                endfor
1010                sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1011                return (sval)
1012        endif
1013        //KillWaves/Z tmp_xw,tmp_yw
1014End
1015
1016////Numerical Recipes routine to calculate the nn(th) stage
1017////refinement of a trapezoid integration
1018////
1019////must be called sequentially from nn=1...n from qtrap()
1020//// to cumulatively refine the integration value
1021////
1022//// in the conversion:
1023//// -- s was replaced with sVal and declared global (rather than static)
1024////  so that the nn-1 value would be available during the nn(th) call
1025////
1026//// -- the specific coefficient wave for func() is passed in as a
1027////  global string (then converted to a wave reference), since
1028////  func() will eventually call sphereForm()
1029////
1030//Function trapzd_USANS_point(fcn,aa,bb,nn)
1031//      FUNCREF SANSModel_proto fcn
1032//      Variable aa,bb,nn
1033//     
1034//      Variable xx,tnm,summ,del
1035//      Variable it,jj,arg1,arg2
1036//      NVAR sVal=sVal          //calling function must initialize this global
1037//      NVAR qval = gEvalQval
1038//      SVAR cwStr = gTrap_CoefStr              //pass in the coefficient wave (string)
1039//      Wave cw=$cwStr                 
1040//      Variable temp=0
1041//      if(nn==1)
1042//              arg1 = qval^2 + aa^2
1043//              arg2 = qval^2 + bb^2
1044//              temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2)))
1045//              sval = temp
1046//              return(sVal)
1047//      else
1048//              it=1
1049//              it= 2^(nn-2)  //done in NR with a bit shift <<=
1050//              tnm = it
1051//              del = (bb - aa)/tnm             //this is the spacing of points to add
1052//              xx = aa+0.5*del
1053//              summ=0
1054//              for(jj=1;jj<=it;jj+=1)
1055//                      arg1 = qval^2 + xx^2
1056//                      summ += fcn(cw,sqrt(arg1))
1057//                      xx += del
1058//              endfor
1059//              sval = 0.5*(sval+(bb-aa)*summ/tnm)      //replaces sval with its refined value
1060//              return (sval)
1061//      endif
1062//     
1063//End
1064
1065// Numerical Recipes routine to calculate the integral of a
1066// specified function, trapezoid rule is used to a user-specified
1067// level of refinement using sequential calls to trapzd()
1068//
1069// in NR, eps and maxIt were global, pass them in here...
1070// eps typically 1e-5
1071// maxit typically 20
1072Function qtrap_USANS(fcn,aa,bb,eps,maxIt)
1073        FUNCREF SANSModelAAO_proto fcn
1074        Variable aa,bb,eps,maxit
1075       
1076        Variable/G sVal=0               //create and initialize what trapzd will return
1077        Variable jj,ss,olds
1078       
1079        olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1080        for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1081                ss = trapzd_USANS(fcn,aa,bb,jj)
1082                if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1083                        return ss
1084                endif
1085                if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1086                        return ss
1087                endif
1088                olds = ss
1089        endfor
1090
1091        Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1092        return(ss)              //should never get here if function is well-behaved
1093       
1094End     
1095
1096//// Numerical Recipes routine to calculate the integral of a
1097//// specified function, trapezoid rule is used to a user-specified
1098//// level of refinement using sequential calls to trapzd()
1099////
1100//// in NR, eps and maxIt were global, pass them in here...
1101//// eps typically 1e-5
1102//// maxit typically 20
1103//Function qtrap_USANS_point(fcn,aa,bb,eps,maxIt)
1104//      FUNCREF SANSModel_proto fcn
1105//      Variable aa,bb,eps,maxit
1106//     
1107//      Variable/G sVal=0               //create and initialize what trapzd will return
1108//      Variable jj,ss,olds
1109//     
1110//      olds = -1e30            //any number that is not the avg. of endpoints of the funciton
1111//      for(jj=1;jj<=maxit;jj+=1)       //call trapzd() repeatedly until convergence
1112//              ss = trapzd_USANS_point(fcn,aa,bb,jj)
1113//              if( abs(ss-olds) < eps*abs(olds) )              // good enough, stop now
1114//                      return ss
1115//              endif
1116//              if( (ss == 0.0) && (olds == 0.0) && (jj>6) )            //no progress?
1117//                      return ss
1118//              endif
1119//              olds = ss
1120//      endfor
1121//
1122//      Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap"
1123//      return(ss)              //should never get here if function is well-behaved
1124//     
1125//End   
1126
1127Proc RRW()
1128        ResetResolutionWaves()
1129End
1130
1131//utility procedures that are currently untied to any actions, although useful...
1132Proc ResetResolutionWaves(str)
1133        String Str
1134        Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","")
1135
1136
1137        Abort "This function is not data floder aware and does nothing..."
1138               
1139        String/G root:gQvals = str[0,strlen(str)-3]+"_q"
1140        String/G root:gSig_Q = str[0,strlen(str)-3]+"sq"
1141        String/G root:gQ_bar = str[0,strlen(str)-3]+"qb"
1142        String/G root:gShadow = str[0,strlen(str)-3]+"fs"
1143       
1144        //touch everything to make sure that the dependencies are
1145        //properly updated - especially the $gQvals reference in the
1146        // dependency assignment
1147        fKillDependencies("Smear*")
1148       
1149        //replace the q-values and intensity (so they're the right length)
1150        fResetSmearedModels("Smear*",root:gQvals)
1151       
1152        fRestoreDependencies("Smear*")
1153End
1154
1155// pass "*" as the matchString to do ALL dependent waves
1156// or "abc*" to get just the matching waves
1157//
1158Function fKillDependencies(matchStr)
1159        String matchStr
1160
1161        String str=WaveList(matchStr, ";", "" ),item,formula
1162        Variable ii
1163       
1164        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1165                item = StringFromList(ii, str ,";")
1166                formula = GetFormula($item)
1167                if(cmpstr("",formula)!=0)
1168                        Printf "wave %s had the formula %s removed\r",item,formula
1169                        Note $item, "FORMULA:"+formula
1170                        SetFormula $item, ""                    //clears the formula
1171                endif
1172        endfor
1173        return(0)
1174end
1175
1176// pass "*" as the matchString to do ALL dependent waves
1177// or "abc*" to get just the matching waves
1178//
1179Function fRestoreDependencies(matchStr)
1180        String matchStr
1181
1182        String str=WaveList(matchStr, ";", "" ),item,formula
1183        Variable ii
1184       
1185        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1186                item = StringFromList(ii, str ,";")
1187                formula = StringByKey("FORMULA", note($item),":",";")
1188                if(cmpstr("",formula)!=0)
1189                        Printf "wave %s had the formula %s restored\r",item,formula
1190                        Note/K $item
1191                        SetFormula $item, formula                       //restores the formula
1192                endif
1193        endfor
1194        return(0)
1195end
1196
1197Function fResetSmearedModels(matchStr,qStr)
1198        String matchStr,qStr
1199
1200        Duplicate/O $qStr root:smeared_qvals   
1201       
1202        String str=WaveList(matchStr, ";", "" ),item,formula
1203        Variable ii
1204       
1205        for(ii=0;ii<ItemsInList(str ,";");ii+=1)
1206                item = StringFromList(ii, str ,";")
1207                formula = StringByKey("FORMULA", note($item),":",";")
1208                if(cmpstr("",formula)!=0)
1209                        Printf "wave %s has been duplicated to gQvals\r",item
1210                        Duplicate/O $qStr $item
1211                        Note $item, "FORMULA:"+formula          //be sure to keep the formula note
1212                        Print "   and the formula is",formula
1213                endif
1214        endfor
1215        return(0)
1216end
Note: See TracBrowser for help on using the repository browser.