source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/GaussUtils.ipf @ 153

Last change on this file since 153 was 145, checked in by ajj, 15 years ago

Changes:

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