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

Last change on this file since 127 was 127, checked in by srkline, 16 years ago

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

24 JUL 2007 SRK

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