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

Last change on this file since 429 was 426, checked in by ajj, 14 years ago

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