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

Last change on this file since 709 was 708, checked in by srkline, 13 years ago

SA_Includes_v410 : now include Smear_2D

PeakGauss_2D, Sphere_2D : included threaded resolution smearing calculations for testing

DataSetHandling? : Included a quick and dirty batch converter for XML->6-col. See the top
of the file for the command to run

GaussUtils? : re-define the ResSemear_2D_AAOStruct. Relocated q-value and phi calculations from
RawWindowHook? to this file so they would be available for reduction and analysis

Smear_2D : now has a generic (non-threaded) smearing routine. Threading must be done in
individual functions since FUNCREF can't be passed to threads (plus a few other issues)

PlotUtils_2D : updated loader for new QxQy? columns. Fixes to Wrapper_2D to enable smeared fits

RawWindowHook? : removed Q-value calculation functions and moved these to GaussUtils?

WriteQIS : now writes out 8-columns for QxQy? data, defining the resolution
function in terms of directions parallel and perpendicular to Q. TEMPORARILY in the data
file an error in intensity is generated that is SQRT(I), being careful to
replace any NaN, inf, or zero with an average error value

MultiScatter_MonteCarlo_2D : 4-processor aware

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