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

Last change on this file since 476 was 444, checked in by srkline, 14 years ago

Changes to the analysis package to add a few more model functions. Documentation and XOPs are to follow later.

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