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

Last change on this file since 599 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

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