source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/GaussUtils.ipf @ 56

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

initial checkin of Analysis v.3.00 files

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