source: sans/Dev/trunk/NCNR_User_Procedures/USANS/LakeDesmearing_JB.ipf @ 404

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

Refactoring USANS functions to use root:Packages:NIST:USANS as their base to comply with canSAS proposal.

  • Property eol-style set to native
File size: 52.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma Version=2.20
3#pragma IgorVersion=6.0
4
5//////////////////////////////////////////////
6// Igor conversion:     03 FEB 06 SRK
7//                      Revision:       03 MAR 06 SRK
8//
9//      Program uses Lake's iterative technique, J. A. Lake, Acta Cryst. 23 (1967) 191-4.
10//      to DESMEAR Infinite slit smeared USANS DATA.
11//
12//      TO SEE DESCRIPTION, CHECK COMPUTER SOFTWARE LOGBOOK, P13,41-46, 50
13//      J. BARKER, August, 2001
14//      Switches from fast to slow convergence near target chi JGB 2/2003
15//
16// steps are:
17// load
18// mask
19// extrapolate (find the power law of the desmeared set = smeared -1)
20// smooth (optional)
21// desmear based on dQv and chi^2 target
22// save result
23//
24/////////
25// Waves produced at each step: (Q/I/S with the following extensions)
26//
27// Load:                Uses -> nothing (Kills old waves)
28//                              Produces -> "_exp" and "_exp_orig"
29//
30// Mask:                Uses -> "_exp"
31//                              Produces -> "_msk"
32//
33// Extrapolate: Uses -> nothing
34//                                      Produces -> nothing
35//
36// Smooth:      Uses -> "_smth" OR "_msk" OR "_exp" (in that order)
37//                              Produces -> "_smth"
38//
39// Desmear:     Uses ->  "_smth" OR "_msk" OR "_exp" (in that order)
40//                              Produces -> "_dsm"
41//
42////////
43
44///***** TO FIX *******
45// - ? power law + background extrapolation (only useful for X-ray data...)
46//
47// see commented code lines for Igor 4 or Igor 5
48// - there are only a few options and calls that are not Igor 4 compatible.
49// Igor 4 routines are currently used.
50
51
52////////////////////////////////////////////////////////////////////////
53
54// main entry routine
55Proc Desmear()
56
57        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
58
59        //check for the correct data folder, initialize if necessary
60        //
61        if(DataFolderExists($(USANSFolder+":DSM")) == 0)
62                Execute "Init_Desmearing()"
63        endif
64       
65        //always initialize these
66        String/G $(USANSFolder+":DSM:gStr1") = ""
67        String/G $(USANSFolder+":DSM:gStr2") = ""
68        String/G $(USANSFolder+":DSM:gIterStr") = ""
69       
70        DoWindow/F Desmear_Graph
71        if(V_flag==0)
72                Execute "Desmear_Graph()"
73        endif
74End
75
76Proc Init_Desmearing()
77        //set up the folder(s) needed
78        NewDataFolder/O $(USANSFolder+":DSM")
79        NewDataFolder/O $(USANSFolder+":myGlobals")             //in case it wasn't created elsewhere
80       
81        String/G $(USANSFolder+":DSM:gCurFile")=""
82       
83        Variable/G $(USANSFolder+":DSM:gMaxFastIter") = 100                     //max number of iter in Fast convergence
84        Variable/G $(USANSFolder+":DSM:gMaxSlowIter") = 10000
85       
86        Variable/G $(USANSFolder+":DSM:gNptsExtrap") = 10               //points for high q extrapolation
87        Variable/G $(USANSFolder+":DSM:gChi2Target") = 1                //chi^2 target
88        Variable/G $(USANSFolder+":DSM:gPowerM") = -4
89        Variable/G $(USANSFolder+":DSM:gDqv") = 0.117                   //2005 measured slit height - see John
90        Variable/G $(USANSFolder+":DSM:gNq") = 1
91        Variable/G $(USANSFolder+":DSM:gS") = 0         // global varaible for Midpnt()
92        Variable/G $(USANSFolder+":DSM:gSmoothFac")=0.03
93       
94        Variable/G $(USANSFolder+":DSM:gChi2Final") = 0         //chi^2 final
95        Variable/G $(USANSFolder+":DSM:gIterations") = 0                //total number of iterations
96       
97        String/G $(USANSFolder+":DSM:gStr1") = ""                               //information strings
98        String/G $(USANSFolder+":DSM:gStr2") = ""
99        String/G $(USANSFolder+":DSM:gIterStr") = ""
100        Variable/G $(USANSFolder+":DSM:gChi2Smooth") = 0
101       
102        Variable/G $(USANSFolder+":DSM:gFreshMask")=1
103End
104
105//////////// Lake desmearing routines
106
107//      Smears guess Y --> YS using weighting array
108Function DSM_Smear(Y_wave,Ys,NQ,FW)
109        Wave Y_wave,Ys
110        Variable nq
111        Wave FW
112       
113        Variable ii,jj,summ
114        for(ii=0;ii<nq;ii+=1)
115                summ=0
116                for(jj=0;jj<nq;jj+=1)
117                        summ = summ + Y_wave[jj]*FW[ii][jj]
118                endfor
119                Ys[ii] = summ
120        endfor
121       
122        Return (0)
123End
124
125//      CALCULATES CHI^2
126FUNCTION CHI2(ys_new,ys_exp,w,NQ)
127        Wave ys_new,ys_exp,w
128        Variable nq
129
130        Variable CHI2,summ,ii
131
132        SUMm = 0.0
133        for(ii=0;ii<nq;ii+=1)
134//              if(numtype(YS_EXP[ii]) == 0 && numtype(YS_NEW[ii]) == 0)
135                SUMm=SUMm+W[ii]*(YS_EXP[ii]-YS_NEW[ii])^2
136//              endif
137        endfor
138       
139        CHI2=SUMm/(NQ-1)
140       
141        RETURN CHI2
142END
143
144//      Routine calculates the weights needed to convert a table
145//      representation of a scattering function I(q) into the infinite
146//      slit smeared table represented as Is(q)
147//      Is(qi) = Sum(j=i-n) Wij*I(qj)
148//      Program assumes data is extrapolated from qn to dqv using powerlaw
149//      I(q) = Aq^m
150//      Required input:
151//      N       : Number of data points {in Global common block }
152//      q(N)    : Array of q values     {in Global common block }
153//      dqv     : Limit of slit length  {in Global common block }
154//      m       : powerlaw for extrapolation.
155//
156//      Routine output:
157//      W(N,N)  : Weights array
158//
159//      The weights obtained from this routine can be used to calculate
160//      Is(q) for any I(q) function.
161//
162//      Calculation based upon linear interpolation between points.
163//      See p.44, Computer Software Log book.
164//      11/00 JGB
165Function Weights_L(m,FW,Q_exp)          //changed input w to FW (FW is two dimensional)
166        Variable m
167        Wave FW,Q_exp
168
169        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
170       
171        NVAR dqv = $(USANSFolder+":DSM:gDqv")
172        NVAR NN = $(USANSFolder+":DSM:gNq")
173
174//      Calculate Remainder fractions and put into separate array.
175        Variable lower,ii,ss,jj
176        Duplicate/O Q_exp $(USANSFolder+":DSM:R_wave")
177        wave r_wave = $(USANSFolder+":DSM:R_wave")
178       
179//      Make/O/D/N=75 root:DSM:SS_save  //debug
180//      wave SS_save = root:DSM:SS_save
181       
182        Print "calculating remainders by integration..."
183        for(ii=0;ii<NN-1;ii+=1)
184                lower = sqrt(Q_exp[NN-1]^2-Q_exp[ii]^2)
185                ss = Qromo(lower,dqv,Q_exp[ii],m)               //new parameter list for Qromo call
186//              SS_save[ii] = ss
187                R_wave[ii] = (Q_exp[NN-1]^(-m)/dqv)*SS
188//              Printf "I = %d R_wave[ii] =%g\r",ii,R_wave[ii]
189        endfor
190       
191        lower = 0.0
192        ss = Qromo(lower,dqv,Q_exp[NN-1],m)             //new parameter list for Qromo call
193        R_wave[NN-1] = (Q_exp[NN-1]^(-m)/dqv)*SS
194//      Printf "I = %d R_wave[ii] =%g\r",NN-1,R_wave[NN-1]
195//      SS_save[ii] = ss
196
197//      Make/O/D/N=(75,75) root:DSM:IG_save             //debug
198//      wave IG_save=root:DSM:IG_save
199       
200//      Zero weight matrix... then fill it
201        FW = 0
202        Print "calculating full weights...."
203//      Calculate weights
204        for(ii=0;ii<NN;ii+=1)
205                for(jj=ii+1;jj<NN-1;jj+=1)
206                        FW[ii][jj] = DU(ii,jj)*(1.0+Q_exp[jj]/(Q_exp[jj+1]-Q_exp[jj]))
207                        FW[ii][jj] -= (1.0/(Q_exp[jj+1]-Q_exp[jj]))*IG(ii,jj)
208                        FW[ii][jj] -= DU(ii,jj-1)*Q_exp[jj-1]/(Q_exp[jj]-Q_exp[jj-1])
209                        FW[ii][jj] += (1.0/(Q_exp[jj]-Q_exp[jj-1]))*IG(ii,jj-1)
210                        FW[ii][jj] *= (1.0/dqv)                 
211//              Printf "FW[%d][%d] = %g\r",ii,jj,FW[ii][jj]
212                endfor
213        endfor
214//
215//      special case: I=J,I=N
216        for(ii=0;ii<NN-1;ii+=1)
217                FW[ii][ii] = DU(ii,ii)*(1.0+Q_exp[ii]/(Q_exp[ii+1]-Q_exp[ii]))
218                FW[ii][ii] -= (1.0/(Q_exp[ii+1]-Q_exp[ii]))*IG(ii,ii)
219                FW[ii][ii] *= (1.0/dqv)
220//              Printf "FW[%d][%d] = %g\r",ii,jj,FW[ii][jj]
221//     following line corrected for N -> NN-1 since Igor indexes from 0! (Q_exp[NN] DNE!)
222                FW[ii][NN-1] = -DU(ii,NN-2)*Q_exp[NN-2]/(Q_exp[NN-1]-Q_exp[NN-2])
223                FW[ii][NN-1] += (1.0/(Q_exp[NN-1]-Q_exp[NN-2]))*IG(ii,NN-2)
224                FW[ii][NN-1] *= (1.0/dqv)
225                FW[ii][NN-1] += R_wave[ii]
226//              Printf "FW[%d][%d] = %g\r",ii,NN-1,FW[ii][NN-1]
227        endfor
228//
229//      special case: I=J=N
230        FW[NN-1][NN-1] = R_wave[NN-1]
231//
232        SetDataFolder root:
233        Return (0)
234End
235
236////
237Function Integrand_dsm(u,qi,m)
238        Variable u,qi,m
239       
240        Variable integrand
241        Integrand = (u^2+qi^2)^(m/2)
242        return integrand
243end
244
245///
246Function DU(ii,jj)
247        Variable ii,jj
248
249        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
250       
251        Wave Q_exp = $(USANSFolder+":DSM:Q_exp")
252        Variable DU
253
254//      Wave DU_save=root:DSM:DU_save
255        If (ii == jj)
256          DU = sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2)
257        Else
258          DU = sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2) - sqrt(Q_exp[jj]^2-Q_exp[ii]^2)
259        EndIf
260       
261//      DU_save[ii][jj] = DU
262        Return DU
263End
264
265Function IG(ii,jj)
266        Variable ii,jj
267       
268        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
269       
270        Wave Q_exp=$(USANSFolder+":DSM:Q_exp")
271        Variable IG,UL,UU
272
273//      WAVE IG_save = root:DSM:IG_save
274        If (ii == jj)
275          UL=0.0
276        Else
277          UL=sqrt(Q_exp[jj]^2-Q_exp[ii]^2)
278        EndIf
279        UU=sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2)
280       
281        //in FORTRAN log = natural log....     
282        IG = UU*Q_exp[jj+1]+Q_exp[ii]^2*ln(UU+Q_exp[jj+1])
283        IG -= UL*Q_exp[jj]
284        IG -= Q_exp[ii]^2*ln(UL+Q_exp[jj])
285        IG *= 0.5
286//      IG_save[ii][jj] = IG
287       
288        Return IG
289End
290
291//
292Function FF(ii,jj)
293        Variable ii,jj
294
295        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
296
297        Variable FF
298        NVAR dqv = $(USANSFolder+":DSM:gDqv")
299
300        FF = (1.0/dqv)*(0.5+HH(ii,jj))
301        Return FF
302End
303
304//
305Function GG(ii,jj)
306        Variable ii,jj
307
308        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
309
310        Variable GG
311        NVAR dqv = $(USANSFolder+":DSM:gDqv")
312       
313        GG = (1.0/dqv)*(0.5-HH(ii,jj))
314        Return  GG
315End
316//
317Function HH(ii,jj)
318        Variable ii,jj
319
320        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
321
322        Wave Q_exp=$(USANSFolder+":DSM:Q_exp")
323        Variable HH
324       
325        HH = 0.5*(Q_exp[jj+1]+Q_exp[jj])/(Q_exp[jj+1]-Q_exp[jj])
326        HH -= (1.0/(Q_exp[jj+1]-Q_exp[jj]))*(CC(ii,jj+1)-CC(ii,jj))
327        return HH
328End
329//
330//
331Function CC(ii,jj)
332        Variable ii,jj
333
334        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
335       
336        wave Q_exp = $(USANSFolder+":DSM:Q_exp")
337        Variable CC
338       
339        If (ii == jj)
340          CC = 0.0
341        Else
342          CC = (Q_exp[jj]-Q_exp[ii])^0.5
343        EndIf
344        Return CC
345End
346
347// QROMO is a gerneric NR routine that takes function arguments
348// Call Qromo(Integrand,lower,dqv,Q_exp(N),m,SS,Midpnt)
349// -- here, it is always called with Integrand and Midpnt as the functions
350// -- so rewrite in a simpler way....
351//
352// SS is the returned value?
353// H_wave, S_wave (original names H,S)
354//
355// modified using c-version in NR (pg 143)
356Function QROMO(A,B,qi,m)
357        Variable A,B,qi,m
358
359        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
360       
361        Variable EPS,JMAX,JMAXP,KM,K
362        EPS=1.E-6
363        JMAX=14
364        KM=4
365        K=KM+1
366       
367        Make/O/D/N=(JMAX) $(USANSFolder+":DSM:S_wave")
368        Make/O/D/N=(JMAX+1) $(USANSFolder+":DSM:H_wave")
369        wave S_wave=$(USANSFolder+":DSM:S_wave")
370        wave H_wave=$(USANSFolder+":DSM:H_wave")
371        S_wave=0
372        H_wave=0
373       
374        H_Wave[0] = 1
375       
376        variable jj,SS,DSS
377       
378        for(jj=0;jj<jmax;jj+=1)
379                S_wave[jj] = Midpnt(A,B,jj,qi,m)                //remove FUNC, always call Integrand from within Midpnt
380                IF (jj>=KM)             //after 1st 5 points calculated
381                    POLINT(H_wave,S_wave,jj-KM,KM,0.0,SS,DSS)           //ss, dss returned
382                    IF (ABS(DSS) < EPS*ABS(SS))
383                        RETURN ss
384                    endif
385                ENDIF
386//        S_wave[jj+1]=S_wave[jj]
387          H_wave[jj+1]=H_wave[jj]/9.
388        endfor
389       
390        DoAlert 0,"Too many steps in QROMO"
391        return 1                        //error if you get here
392END
393
394//
395// see NR pg 109
396// - Given input arrays xa[1..n] and ya[1..n], and a given value x
397// return a value y and an error estimate dy. if P(x) is the polynomial of
398// degree N-1 such that P(xai) = yai, then the returned value y=P(x)
399//
400// arrays XA[] and YA[] are passed in with an offset of the index
401// of where to start the interpolation
402Function POLINT(XA,YA,offset,N,X,Y,DY)
403        Wave XA,YA
404        Variable offset,N,X,&Y,&DY
405       
406        Variable ii,mm,nmax,ns,dif,den,ho,hp,wi,dift
407
408        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
409       
410        NMAX=10
411       
412        Make/O/D/N=(NMAX) $(USANSFolder+":DSM:Ci"),$(USANSFolder+":DSM:Di")
413        wave Ci = $(USANSFolder+":DSM:Ci")
414        wave Di = $(USANSFolder+":DSM:Di")
415
416        NS=1
417        DIF=ABS(X-XA[0])
418        for(ii=0;ii<n;ii+=1)
419                DIFT=ABS(X-XA[ii+offset])
420                IF (DIFT < DIF)
421                        NS=ii
422                        DIF=DIFT
423                ENDIF
424                Ci[ii]=YA[ii+offset]
425                Di[ii]=YA[ii+offset]
426        endfor
427       
428        Y=YA[NS+offset]
429        NS=NS-1
430        for(mm=1;mm<n-1;mm+=1)
431                for(ii=0;ii<(n-mm);ii+=1)
432                        HO=XA[ii+offset]-X
433                        HP=XA[ii+offset+mm]-X
434                        Wi=Ci[ii+1]-Di[ii]
435                        DEN=HO-HP
436                        IF(DEN == 0.)
437                                print "den == 0 in POLINT - ERROR!!!"
438                        endif
439                        DEN=Wi/DEN
440                        Di[ii]=HP*DEN
441                        Ci[ii]=HO*DEN
442                endfor  //ii
443                IF (2*NS < (N-mm) )
444                        DY=Ci[NS+1]
445                ELSE
446                        DY=Di[NS]
447                        NS=NS-1
448                ENDIF
449                Y=Y+DY
450        endfor  //mm
451        RETURN (0)              //y and dy are returned as pass-by-reference
452END
453
454//
455// FUNC is always Integrand()
456// again, see the c-version, NR pg 142
457Function MIDPNT(A,B,N,qi,m)
458        Variable A,B,N,qi,m
459       
460        Variable it,tnm,del,ddel,x,summ,jj
461
462        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
463       
464        NVAR S_ret = $(USANSFolder+":DSM:gS")
465       
466        IF (N == 0)
467                S_ret=(B-A)*Integrand_dsm(0.5*(A+B),qi,m)
468//              Print "N==0, S_ret = ",s_ret
469                return(S_ret)
470                        //IT=1
471        ELSE
472                //AJJ This code is confusing!
473//              it = 1
474//              for(jj=1;jj<n-1;jj+=1)
475//                      it *= 3
476//              endfor
477                //AJJ Equivalent and simpler   
478                it = 3^(N-1)   
479                //
480                TNM=IT
481                DEL=(B-A)/(3.*TNM)
482                DDEL=DEL+DEL
483                X=A+0.5*DEL
484                SUMm=0.
485                for(jj=1;jj<=it;jj+=1)
486                        SUMm=SUMm+Integrand_dsm(X,qi,m)
487                        X=X+DDEL
488                        SUMm=SUMm+Integrand_dsm(X,qi,m)
489                        X=X+DEL
490                endfor
491                S_ret=(S_ret+(B-A)*SUMm/TNM)/3.
492                IT=3*IT
493        ENDIF
494        return S_ret
495END
496
497//////////// end of Lake desmearing routines
498
499// - Almost everything below here is utility routines for data handling, panel
500// controls, and all the fluff that makes this useable. Note that this is
501// a lot more code than the actual guts of the Lake method!
502// (the guts of the iteration are in the DemsearButtonProc)
503
504
505// add three "fake" resolution columns to the output file so that it
506// looks like a regular SANS data set, but make the sigQ column 100x smaller than Q
507// so that there is effectively no resolution smearing. Set Qbar = Q, and fs = 1
508//
509// SRK 06 FEB 06
510//
511Function WriteUSANSDesmeared(fullpath,lo,hi,dialog)
512        String fullpath
513        Variable lo,hi,dialog           //=1 will present dialog for name
514       
515        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
516       
517        String termStr="\r\n"
518        String destStr = USANSFolder+":DSM:"
519        String formatStr = "%15.6g %15.6g %15.6g %15.6g %15.6g %15.6g"+termStr
520       
521        Variable refNum,integer,realval
522       
523        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
524        WAVE Q_dsm =$(destStr + "Q_dsm")
525        WAVE I_dsm=$(destStr + "I_dsm")
526        WAVE S_dsm=$(destStr + "S_dsm")
527       
528        //check each wave
529        If(!(WaveExists(Q_dsm)))
530                Abort "Q_dsm DNExist in WriteUSANSDesmeared()"
531        Endif
532        If(!(WaveExists(I_dsm)))
533                Abort "I_dsm DNExist in WriteUSANSDesmeared()"
534        Endif
535        If(!(WaveExists(S_dsm)))
536                Abort "S_dsm DNExist in WriteUSANSDesmeared()"
537        Endif
538       
539        // 06 FEB 06 SRK
540        // make dummy waves to hold the "fake" resolution, and write it as the last 3 columns
541        //
542        Duplicate/O Q_dsm,res1,res2,res3
543        res3 = 1                // "fake" beamstop shadowing
544        res1 /= 100             //make the sigmaQ so small that there is no smearing
545       
546        if(dialog)
547                Open/D refnum as fullpath+".dsm"                //won't actually open the file
548                If(cmpstr(S_filename,"")==0)
549                        //user cancel, don't write out a file
550                        Close/A
551                        Abort "no data file was written"
552                Endif
553                fullpath = S_filename
554        Endif
555       
556        //write out partial set?
557        Duplicate/O Q_dsm,tq,ti,te
558        ti=I_dsm
559        te=S_dsm
560        if( (lo!=hi) && (lo<hi))
561                redimension/N=(hi-lo+1) tq,ti,te,dumWave                //lo to hi, inclusive
562                tq=Q_dsm[p+lo]
563                ti=I_dsm[p+lo]
564                te=S_dsm[p+lo]
565        endif
566       
567        //tailor the output given the type of data written out...
568        String samStr="",dateStr="",str1,str2
569       
570        NVAR m = $(USANSFolder+":DSM:gPowerM")                          // power law exponent
571        NVAR chiFinal = $(USANSFolder+":DSM:gChi2Final")                //chi^2 final
572        NVAR iter = $(USANSFolder+":DSM:gIterations:")          //total number of iterations
573       
574        //get the number of spline passes from the wave note
575        String noteStr
576        Variable boxPass,SplinePass
577        noteStr=note(I_dsm)
578        BoxPass = NumberByKey("BOX", noteStr, "=", ";")
579        splinePass = NumberByKey("SPLINE", noteStr, "=", ";")
580       
581        samStr = fullpath
582        dateStr="CREATED: "+date()+" at  "+time()
583        sprintf str1,"Chi^2 = %g   PowerLaw m = %4.2f   Iterations = %d",chiFinal,m,iter
584        sprintf str2,"%d box smooth passes and %d smoothing spline passes",boxPass,splinePass
585       
586        //actually open the file
587        Open refNum as fullpath
588       
589        fprintf refnum,"%s"+termStr,samStr
590        fprintf refnum,"%s"+termStr,str1
591        fprintf refnum,"%s"+termStr,str2
592        fprintf refnum,"%s"+termStr,dateStr
593       
594        wfprintf refnum, formatStr, tq,ti,te,res1,res2,res3
595       
596        Close refnum
597       
598        Killwaves/Z ti,tq,te,res1,res2,res3
599       
600        Return(0)
601End
602
603/// procedures to do the extrapolation to high Q
604// very similar to the procedures in the Invariant package
605//
606//create the wave to extrapolate
607// w is the input q-values
608Function DSM_SetExtrWaves(w)
609        Wave w
610
611        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
612
613        Variable num_extr=25
614       
615        SetDataFolder $(USANSFolder+":DSM")
616       
617        Make/O/D/N=(num_extr) extr_hqq,extr_hqi
618        extr_hqi=1              //default values
619
620        //set the q-range
621        Variable qmax,num
622//      qmax=0.03
623
624        num=numpnts(w)
625        qmax=6*w[num-1]
626       
627        extr_hqq = w[num-1] + x * (qmax-w[num-1])/num_extr
628       
629        SetDataFolder root:
630        return(0)
631End
632
633//creates I_ext,Q_ext,S_ext with extended q-values
634//
635// if num_extr == 0 , the input waves are returned as _ext
636//
637// !! uses simple linear extrapolation at low Q - just need something
638// reasonable in the waves to keep from getting smoothing artifacts
639// - uses power law at high q
640Function ExtendToSmooth(qw,iw,sw,nbeg,nend,num_extr)
641        Wave qw,iw,sw
642        Variable nbeg,nend,num_extr
643
644        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
645       
646        Variable/G V_FitMaxIters=300
647        Variable/G V_fitOptions=4               //suppress the iteration window
648        Variable num_new,qmax,num,qmin
649        num=numpnts(qw)
650       
651        num_new = num + 2*num_extr
652       
653//      Print "num,num_new",num,num_new
654       
655        SetDataFolder $(USANSFolder+":DSM")
656        Make/O/D/N=(num_new) Q_ext,I_ext,S_ext
657       
658        if(num_extr == 0)               //the extended waves are the input, get out
659                Q_ext = qw
660                I_ext = iw
661                S_ext = sw
662                setDatafolder root:
663                return (0)
664        endif
665       
666        //make the extensions
667        Make/O/D/N=(num_extr) hqq,hqi,lqq,lqi
668        hqi=1           //default values
669        lqi=0
670        //set the q-range based on the high/low values of the data set
671//      qmin=1e-6
672        qmin= qw[0]*0.5
673        qmax = qw[num-1]*2
674       
675        hqq = qw[num-1] + x * (qmax-qw[num-1])/num_extr
676        lqq = qmin + x * (qw[0]-qmin)/num_extr
677       
678        //do the fits
679        Duplicate/O iw dummy                    //use this as the destination to suppress fit_ from being appended
680       
681        //Use simple linear fits        line: y = K0+K1*x
682        CurveFit/Q line iw[0,(nbeg-1)] /X=qw /W=sw /D=dummy
683        Wave W_coef=W_coef
684        lqi=W_coef[0]+W_coef[1]*lqq
685//     
686// Guinier or Power-law fits
687//             
688//      Make/O/D G_coef={100,-100}              //input
689//      FuncFit DSM_Guinier_Fit G_coef iw[0,(nbeg-1)] /X=qw /W=sw /D
690//      lqi= DSM_Guinier_Fit(G_coef,lqq)
691       
692//      Printf "I(q=0) = %g (1/cm)\r",G_coef[0]
693//      Printf "Rg = %g (A)\r",sqrt(-3*G_coef[1])
694               
695        Make/O/D P_coef={0,1,-4}                        //input  --- (set background to zero and hold fixed)
696        CurveFit/Q/N/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /W=sw /D=dummy
697        hqi=P_coef[0]+P_coef[1]*hqq^P_coef[2]
698       
699//      Printf "Power law exponent = %g\r",P_coef[2]
700//      Printf "Pre-exponential = %g\r",P_coef[1]
701       
702        // concatenate the extensions
703        Q_ext[0,(num_extr-1)] = lqq[p]
704        Q_ext[num_extr,(num_extr+num-1)] = qw[p-num_extr]
705        Q_ext[(num_extr+num),(2*num_extr+num-1)] = hqq[p-num_extr-num]
706        I_ext[0,(num_extr-1)] = lqi[p]
707        I_ext[num_extr,(num_extr+num-1)] = iw[p-num_extr]
708        I_ext[(num_extr+num),(2*num_extr+num-1)] = hqi[p-num_extr-num]
709        S_ext[0,(num_extr-1)] = sw[0]
710        S_ext[num_extr,(num_extr+num-1)] = sw[p-num_extr]
711        S_ext[(num_extr+num),(2*num_extr+num-1)] = sw[num-1]
712       
713        killwaves/z dummy
714        SetDataFolder root:
715        return(0)
716End
717
718// pass in the (smeared) experimental data and find the power-law extrapolation
719// 10 points is usually good enough, unless the data is crummy
720// returns the prediction for the exponent
721//
722Function DSM_DoExtrapolate(qw,iw,sw,nend)
723        Wave qw,iw,sw
724        Variable nend
725
726        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
727       
728        Setdatafolder $(USANSFolder+":DSM")
729       
730//      Wave extr_lqi=extr_lqi
731//      Wave extr_lqq=extr_lqq
732        Wave extr_hqi=extr_hqi
733        Wave extr_hqq=extr_hqq
734        Variable/G V_FitMaxIters=300
735        Variable/G V_fitOptions=4               //suppress the iteration window
736        Variable num=numpnts(iw),retVal
737       
738        Make/O/D P_coef={0,1,-4}                        //input
739//      Make/O/T Constr={"K2<0","K2 > -20"}
740        //(set background to zero and hold fixed)
741        CurveFit/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /W=sw /D
742        extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
743       
744        // for the case of data with a background
745//      Make/O/D P_coef={0,1,-4}                        //input
746//      //(set background to iw[num-1], let it be a free parameter
747//      P_coef[0] = iw[num-1]
748//      CurveFit Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /W=sw /D
749//      extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
750       
751       
752        Printf "Smeared Power law exponent = %g\r",P_coef[2]
753        Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
754       
755        retVal = P_coef[2]-1
756        SetDataFolder root:
757        return(retVal)
758End
759
760Function DSM_Guinier_Fit(w,x) //: FitFunc
761        Wave w
762        Variable x
763       
764        //fit data to I(q) = A*exp(B*q^2)
765        // (B will be negative)
766        //two parameters
767        Variable a,b,ans
768        a=w[0]
769        b=w[1]
770        ans = a*exp(b*x*x)
771        return(ans)
772End
773
774Proc Desmear_Graph()
775        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
776
777        PauseUpdate; Silent 1           // building window...
778        Display /W=(5,44,408,558) /K=1
779        ModifyGraph cbRGB=(51664,44236,58982)
780        DoWindow/C Desmear_Graph
781        ControlBar 160
782        // break into tabs
783        TabControl DSM_Tab,pos={5,3},size={392,128},proc=DSM_TabProc
784        TabControl DSM_Tab,labelBack=(49151,49152,65535),tabLabel(0)="Load"
785        TabControl DSM_Tab,tabLabel(1)="Mask",tabLabel(2)="Extrapolate"
786        TabControl DSM_Tab,tabLabel(3)="Smooth",tabLabel(4)="Desmear",value= 0
787       
788        //always visible - revert and save
789        //maybe the wrong place here?
790        Button DSMControlA,pos={225,135},size={80,20},proc=DSM_RevertButtonProc,title="Revert"
791        Button DSMControlA,help={"Revert the smeared data to its original state and start over"}
792        Button DSMControlB,pos={325,135},size={50,20},proc=DSM_SaveButtonProc,title="Save"
793        Button DSMControlB,help={"Save the desmeared data set"}
794        Button DSMControlC,pos={25,135},size={50,20},proc=DSM_HelpButtonProc,title="Help"
795        Button DSMControlC,help={"Show the help file for desmearing"}
796       
797        // add the controls to each tab ---- all names start with "DSMControl_"
798
799        //tab(0) Load - initially visible
800        Button DSMControl_0a,pos={23,39},size={80,20},proc=DSM_LoadButtonProc,title="Load Data"
801        Button DSMControl_0a,help={"Load slit-smeared USANS data = \".cor\" files"}
802        CheckBox DSMControl_0b,pos={26,74},size={80,14},proc=DSM_LoadCheckProc,title="Log Axes?"
803        CheckBox DSMControl_0b,help={"Toggle Log/Lin Q display"},value= 1
804        TitleBox DSMControl_0c,pos={120,37},size={104,19},font="Courier",fSize=10
805        TitleBox DSMControl_0c,variable= $(USANSFolder+":DSM:gStr1")
806        //second message string not used currently
807//      TitleBox DSMControl_0d,pos={120,57},size={104,19},font="Courier",fSize=10
808//      TitleBox DSMControl_0d,variable= root:DSM:gStr2
809       
810        //tab(1) Mask
811        Button DSMControl_1a,pos={20,35},size={90,20},proc=DSM_MyMaskProc,title="Mask Point"            //bMask
812        Button DSMControl_1a,help={"Toggles the masking of the selected data point"}
813        Button DSMControl_1a,disable=1
814        Button DSMControl_1b,pos={20,65},size={140,20},proc=DSM_MaskGTCursor,title="Mask Q >= Cursor"           //bMask
815        Button DSMControl_1b,help={"Toggles the masking of all q-values GREATER than the current cursor location"}
816        Button DSMControl_1b,disable=1
817        Button DSMControl_1c,pos={20,95},size={140,20},proc=DSM_MaskLTCursor,title="Mask Q <= Cursor"           //bMask
818        Button DSMControl_1c,help={"Toggles the masking of all q-values LESS than the current cursor location"}
819        Button DSMControl_1c,disable=1
820        Button DSMControl_1d,pos={180,35},size={90,20},proc=DSM_ClearMaskProc,title="Clear Mask"                //bMask
821        Button DSMControl_1d,help={"Clears all mask points"}
822        Button DSMControl_1d,disable=1
823//      Button DSMControl_1b,pos={144,66},size={110,20},proc=DSM_MaskDoneButton,title="Done Masking"
824//      Button DSMControl_1b,disable=1
825       
826        //tab(2) Extrapolate
827        Button DSMControl_2a,pos={31,42},size={90,20},proc=DSM_ExtrapolateButtonProc,title="Extrapolate"
828        Button DSMControl_2a,help={"Extrapolate the high-q region with a power-law"}
829        Button DSMControl_2a,disable=1
830        SetVariable DSMControl_2b,pos={31,70},size={100,15},title="# of points"
831        SetVariable DSMControl_2b,help={"Set the number of points for the power-law extrapolation"}
832        SetVariable DSMControl_2b,limits={5,100,1},value=  $(USANSFolder+":DSM:gNptsExtrap")
833        SetVariable DSMControl_2b,disable=1
834        CheckBox DSMControl_2c,pos={157,45},size={105,14},proc=DSM_ExtrapolationCheckProc,title="Show Extrapolation"
835        CheckBox DSMControl_2c,help={"Show or hide the high q extrapolation"},value= 1
836        CheckBox DSMControl_2c,disable=1
837        SetVariable DSMControl_2d,pos={31,96},size={150,15},title="Power Law Exponent"
838        SetVariable DSMControl_2d,help={"Power Law exponent from the fit = the DESMEARED slope - override as needed"}
839        SetVariable DSMControl_2d format="%5.2f"
840        SetVariable DSMControl_2d,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gPowerM")
841        SetVariable DSMControl_2d,disable=1
842       
843        //tab(3) Smooth
844        Button DSMControl_3a,pos={34,97},size={70,20},proc=DSM_SmoothButtonProc,title="Smooth"
845        Button DSMControl_3a,disable=1
846                //BoxCheck
847        CheckBox DSMControl_3b,pos={34,39},size={35,14},title="Box",value= 1
848        CheckBox DSMControl_3b,help={"Use a single pass of 3-point box smoothing"}
849        CheckBox DSMControl_3b,disable=1
850                //SSCheck
851        CheckBox DSMControl_3c,pos={34,60},size={45,14},title="Spline",value= 0
852        CheckBox DSMControl_3c,help={"Use a single pass of a smoothing spline"}
853        CheckBox DSMControl_3c,disable=1
854                //extendCheck
855        CheckBox DSMControl_3d,pos={268,60},size={71,14},title="Extend Data"
856        CheckBox DSMControl_3d,help={"extends the data at both low q and high q to avoid end effects in smoothing"}
857        CheckBox DSMControl_3d,value= 0
858        CheckBox DSMControl_3d,disable=1
859        Button DSMControl_3e,pos={125,97},size={90,20},proc=DSM_SmoothUndoButtonProc,title="Start Over"
860        Button DSMControl_3e,help={"Start the smoothing over again without needing to re-mask the data set"}
861        Button DSMControl_3e,disable=1
862        SetVariable DSMControl_3f,pos={94,60},size={150,15},title="Smoothing factor"
863        SetVariable DSMControl_3f,help={"Smoothing factor for the smoothing spline"}
864        SetVariable DSMControl_3f format="%5.4f"
865        SetVariable DSMControl_3f,limits={0.01,2,0.01},value= $(USANSFolder+"DSM:gSmoothFac")
866        SetVariable DSMControl_3f,disable=1
867        CheckBox DSMControl_3g,pos={268,39},size={90,14},title="Log-scale smoothing?"
868        CheckBox DSMControl_3g,help={"Use log-scaled intensity during smoothing (reverts to linear if negative intensity points found)"}
869        CheckBox DSMControl_3g,value=0
870        CheckBox DSMControl_3g,disable=1
871        ValDisplay DSMControl_3h pos={235,97},title="Chi^2",size={80,20},value=root:DSM:gChi2Smooth
872        ValDisplay DSMControl_3h,help={"This is the Chi^2 value for the smoothed data vs experimental data"}
873        ValDisplay DSMControl_3h,disable=1
874       
875        //tab(4) Desmear
876        Button DSMControl_4a,pos={35,93},size={80,20},proc=DSM_DesmearButtonProc,title="Desmear"
877        Button DSMControl_4a,help={"Do the desmearing - the result is in I_dsm"}
878        Button DSMControl_4a,disable=1
879        SetVariable DSMControl_4b,pos={35,63},size={120,15},title="Chi^2 target"
880        SetVariable DSMControl_4b,help={"Set the targetchi^2 for convergence (recommend chi^2=1)"}
881        SetVariable DSMControl_4b,limits={0,inf,0.1},value= $(USANSFolder+":DSM:gChi2Target")
882        SetVariable DSMControl_4b,disable=1
883        SetVariable DSMControl_4c,pos={35,35},size={80,15},title="dQv"
884        SetVariable DSMControl_4c,help={"Slit height as read in from the data file. 0.117 is the NIST value, override if necessary"}
885        SetVariable DSMControl_4c,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gDqv")
886        SetVariable DSMControl_4c,disable=1
887        TitleBox DSMControl_4d,pos={160,37},size={104,19},font="Courier",fSize=10
888        TitleBox DSMControl_4d,variable= $(USANSFolder+":DSM:gIterStr")
889        TitleBox DSMControl_4d,disable=1
890       
891       
892        SetDataFolder root:
893EndMacro
894
895// function to control the drawing of buttons in the TabControl on the main panel
896// Naming scheme for the buttons MUST be strictly adhered to... else buttons will
897// appear in odd places...
898// all buttons are named DSMControl_NA where N is the tab number and A is the letter denoting
899// the button's position on that particular tab.
900// in this way, buttons will always be drawn correctly :-)
901//
902Function DSM_TabProc(ctrlName,tab) //: TabControl
903        String ctrlName
904        Variable tab
905
906        String ctrlList = ControlNameList("",";"),item="",nameStr=""
907        Variable num = ItemsinList(ctrlList,";"),ii,onTab
908        for(ii=0;ii<num;ii+=1)
909                //items all start w/"DSMControl_"               //11 characters
910                item=StringFromList(ii, ctrlList ,";")
911                nameStr=item[0,10]
912                if(cmpstr(nameStr,"DSMControl_")==0)
913                        onTab = str2num(item[11])                       //12th is a number
914                        ControlInfo $item
915                        switch(abs(V_flag))     
916                                case 1:
917                                        Button $item,disable=(tab!=onTab)
918                                        break
919                                case 2:
920                                        CheckBox $item,disable=(tab!=onTab)
921                                        break
922                                case 5:
923                                        SetVariable $item,disable=(tab!=onTab)
924                                        break
925                                case 10:       
926                                        TitleBox $item,disable=(tab!=onTab)
927                                        break
928                                case 4:
929                                        ValDisplay $item,disable=(tab!=onTab)
930                                        break
931                                // add more items to the switch if different control types are used
932                        endswitch
933                       
934                endif
935        endfor
936       
937        if(tab==1)
938                DSM_MyMaskProc("")              //start maksing if you click on the tab
939        else
940                DSM_MaskDoneButton("")          //masking is done if you click off the tab
941        endif
942       
943        return 0
944End
945
946Proc AppendSmeared()
947        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
948
949        SetDataFolder $(USANSFolder+":DSM")
950//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0,2) == -1)            //Igor 5
951        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0) == -1)
952                AppendToGraph/W=Desmear_Graph  I_exp_orig vs Q_exp_orig
953                ModifyGraph mode=4,marker=19            //3 is markers, 4 is markers and lines
954                ModifyGraph rgb(I_exp_orig)=(0,0,0)
955                ModifyGraph msize=2,grid=1,log=1,mirror=2,standoff=0,tickunit=1
956                ErrorBars I_exp_orig Y,wave=(S_exp_orig,S_exp_orig)
957                Legend/N=text0/J "\\F'Courier'\\s(I_exp_orig) I_exp_orig"
958                Label left "Intensity"
959                Label bottom "Q (1/A)"
960        endif
961        //always update the textbox - kill the old one first
962        TextBox/K/N=text1
963//      TextBox/C/N=text1/F=0/A=MT/E=2/X=5.50/Y=0.00 root:DSM:gCurFile                  //Igor 5
964        TextBox/C/N=text1/F=0/A=MT/E=1/X=5.50/Y=0.00 $(USANSFolder+"DSM:gCurFile")
965End
966
967Proc AppendMask()
968        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
969
970//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0,2) == -1)                      //Igor 5
971        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0) == -1)
972                SetDataFolder $(USANSFolder+":DSM:")
973                AppendToGraph/W=Desmear_Graph MaskData vs Q_exp_orig
974                ModifyGraph mode(MaskData)=3,marker(MaskData)=8,msize(MaskData)=2.5,opaque(MaskData)=1
975                ModifyGraph rgb(MaskData)=(65535,16385,16385)
976                setdatafolder root:
977        endif
978end
979
980Proc AppendSmoothed()
981        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
982
983//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0,2) == -1)                        //Igor 5
984        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0) == -1)
985                SetDataFolder $(USANSFolder+":DSM:")
986                AppendToGraph/W=Desmear_Graph I_smth vs Q_smth
987                ModifyGraph/W=Desmear_Graph rgb(I_smth)=(3,52428,1),lsize(I_smth)=2
988                setdatafolder root:
989        endif
990end
991
992Function RemoveSmoothed()
993        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
994
995        SetDataFolder $(USANSFolder+":DSM:")
996        RemoveFromGraph/W=Desmear_Graph/Z I_smth
997        setdatafolder root:
998end
999
1000Function RemoveMask()
1001        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1002
1003        SetDataFolder $(USANSFolder+":DSM:")
1004        RemoveFromGraph/W=Desmear_Graph/Z MaskData
1005        setdatafolder root:
1006end
1007
1008Proc AppendDesmeared()
1009        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1010
1011//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0,2) == -1)         //Igor 5
1012        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0) == -1)
1013                SetDataFolder $(USANSFolder+":DSM:")
1014                AppendToGraph/W=Desmear_Graph I_dsm vs Q_dsm
1015                ModifyGraph mode(I_dsm)=3,marker(I_dsm)=19
1016                ModifyGraph rgb(I_dsm)=(1,16019,65535),msize(I_dsm)=2
1017                ErrorBars I_dsm Y,wave=(S_dsm,S_dsm)
1018                setdatafolder root:
1019        endif
1020end
1021
1022Function RemoveDesmeared()
1023        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1024
1025        SetDataFolder $(USANSFolder+":DSM:")
1026        RemoveFromGraph/W=Desmear_Graph/Z I_dsm
1027        setdatafolder root:
1028end
1029
1030Function AppendExtrapolation()
1031        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1032
1033//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0,2) == -1)              //Igor 5
1034        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0) == -1)
1035                SetDataFolder $(USANSFolder+":DSM:")
1036                AppendToGraph/W=Desmear_Graph extr_hqi vs extr_hqq
1037                ModifyGraph/W=Desmear_Graph lSize(extr_hqi)=2
1038                setdatafolder root:
1039        endif
1040end
1041
1042Function RemoveExtrapolation()
1043        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1044
1045        SetDataFolder $(USANSFolder+":DSM:")
1046        RemoveFromGraph/W=Desmear_Graph/Z extr_hqi
1047        setdatafolder root:
1048end
1049
1050// step (1) - read in the data, and plot it
1051// clear out all of the "old" waves, remove them from the graph first
1052// reads in a fresh copy of the data
1053//
1054// produces Q_exp, I_exp, S_exp waves (and originals "_orig")
1055// add a dummy wave note that can be changed on later steps
1056//
1057Function DSM_LoadButtonProc(ctrlName) : ButtonControl
1058        String ctrlName
1059
1060        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1061
1062        String qStr,iStr,sStr,sqStr
1063        Variable nq,dqv,numBad,val
1064       
1065        // remove any of the old traces on the graph and delete the waves
1066        CleanUpJunk()
1067       
1068        SetDataFolder root:
1069        Execute "U_LoadOneDDataWithName(\"\")"
1070        //define the waves that the smoothing will be looking for...
1071        SVAR fname = $(USANSFolder+":myGlobals:gLastFileName")          //this changes as any data is loaded
1072        SVAR curFile = $(USANSFolder+":DSM:gCurFile")                                   //keep this for title, save
1073        curFile = fname
1074       
1075        qStr = CleanupName((fName + "_q"),0)            //the q-wave
1076        iStr = CleanupName((fName + "_i"),0)            //the i-wave
1077        sStr = CleanupName((fName + "_s"),0)            //the s-wave
1078        sqStr = CleanupName((fName + "sq"),0)           //the sq-wave, which should have -dQv as the elements
1079
1080        Duplicate/O $qStr $(USANSFolder+":DSM:Q_exp     ")             
1081        Duplicate/O $iStr $(USANSFolder+":DSM:I_exp")           
1082        Duplicate/O $sStr $(USANSFolder+":DSM:S_exp     ")     
1083        wave Q_exp = $(USANSFolder+":DSM:Q_exp")
1084        Wave I_exp = $(USANSFolder+":DSM:I_exp")
1085        Wave S_exp = $(USANSFolder+":DSM:S_exp")
1086        Wave/Z sigQ = $sqStr
1087       
1088        // remove any negative q-values (and q=0 values!)(and report this)
1089        // ? and trim the low q to be >= 3.0e-5 (1/A), below this USANS is not reliable.
1090        NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,0)
1091        SVAR str1 = $(USANSFolder+":DSM:gStr1")
1092        sprintf str1,"%d negative q-values were removed",numBad
1093       
1094// don't trim off any positive q-values
1095//      val = 3e-5              //lowest "good" q-value from USANS
1096//      NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,val-1e-8)
1097//      SVAR str2 = root:DSM:gStr2
1098//      sprintf str2,"%d q-values below q = %g were removed",numBad,val
1099       
1100        Duplicate/O $(USANSFolder+":DSM:Q_exp") $(USANSFolder+":DSM:Q_exp_orig")
1101        Duplicate/O $(USANSFolder+":DSM:I_exp") $(USANSFolder+":DSM:I_exp_orig")
1102        Duplicate/O $(USANSFolder+":DSM:S_exp") $(USANSFolder+":DSM:S_exp_orig")
1103        wave I_exp_orig = $(USANSFolder+":DSM:I_exp_orig")
1104       
1105        nq = numpnts($(USANSFolder+":DSM:Q_exp"))
1106       
1107        if(WaveExists(sigQ))                    //try to read dQv
1108                dqv = -sigQ[0]
1109//              DoAlert 0,"Found dQv value of " + num2str(dqv)
1110        else
1111                dqv = 0.117
1112        //      dqv = 0.037             //old value
1113                DoAlert 0,"Could not find dQv in the data file - using " + num2str(dqv)
1114        endif
1115        NVAR gDqv = $(USANSFolder+":DSM:gDqv")                          //needs to be global for Weights_L()
1116        NVAR gNq = $(USANSFolder+":DSM:gNq")
1117        //reset the globals
1118        gDqv = dqv
1119        gNq = nq
1120       
1121        // append the (blank) wave note to the intensity wave
1122        Note I_exp,"BOX=0;SPLINE=0;"
1123        Note I_exp_orig,"BOX=0;SPLINE=0;"
1124       
1125        //draw the graph
1126        Execute "AppendSmeared()"
1127       
1128        SetDataFolder root:
1129End
1130
1131// remove any q-values <= val
1132Function RemoveBadQPoints(qw,iw,sw,val)
1133        Wave qw,iw,sw
1134        Variable val
1135       
1136        Variable ii,num,numBad,qval
1137        num = numpnts(qw)
1138       
1139        ii=0
1140        numBad=0
1141        do
1142                qval = qw[ii]
1143                if(qval <= val)
1144                        numBad += 1
1145                else            //keep the points
1146                        qw[ii-numBad] = qval
1147                        iw[ii-numBad] = iw[ii]
1148                        sw[ii-numBad] = sw[ii]
1149                endif
1150                ii += 1
1151        while(ii<num)
1152        //trim the end of the waves
1153        DeletePoints num-numBad, numBad, qw,iw,sw
1154        return(numBad)
1155end
1156
1157// if mw = Nan, keep the point, if a numerical value, delete it
1158Function RemoveMaskedPoints(mw,qw,iw,sw)
1159        Wave mw,qw,iw,sw
1160       
1161        Variable ii,num,numBad,mask
1162        num = numpnts(qw)
1163       
1164        ii=0
1165        numBad=0
1166        do
1167                mask = mw[ii]
1168                if(numtype(mask) != 2)          //if not NaN
1169                        numBad += 1
1170                else            //keep the points that are NaN
1171                        qw[ii-numBad] = qw[ii]
1172                        iw[ii-numBad] = iw[ii]
1173                        sw[ii-numBad] = sw[ii]
1174                endif
1175                ii += 1
1176        while(ii<num)
1177        //trim the end of the waves
1178        DeletePoints num-numBad, numBad, qw,iw,sw
1179        return(numBad)
1180end
1181
1182// produces the _msk waves that have the new number of data points
1183//
1184Function DSM_MaskDoneButton(ctrlName) : ButtonControl
1185        String ctrlName
1186
1187        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1188
1189//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1190        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1191        if(!aExists)
1192                return(1)               //possibly reverted data, no cursor, no Mask wave
1193        endif
1194       
1195        Duplicate/O $(USANSFolder+":DSM:Q_exp_orig"),$(USANSFolder+":DSM:Q_msk")
1196        Duplicate/O $(USANSFolder+":DSM:I_exp_orig"),$(USANSFolder+":DSM:I_msk")
1197        Duplicate/O $(USANSFolder+":DSM:S_exp_orig"),$(USANSFolder+":DSM:S_msk")
1198        Wave Q_msk=$(USANSFolder+":DSM:Q_msk")
1199        Wave I_msk=$(USANSFolder+":DSM:I_msk")
1200        Wave S_msk=$(USANSFolder+":DSM:S_msk")
1201       
1202        //finish up - trim the data sets and reassign the working set
1203        Wave MaskData=$(USANSFolder+":DSM:MaskData")
1204       
1205        RemoveMaskedPoints(MaskData,Q_msk,I_msk,S_msk)
1206
1207        //reset the number of points
1208        NVAR gNq = $(USANSFolder+":DSM:gNq")
1209        gNq = numpnts(Q_msk)
1210       
1211        Cursor/K A
1212        HideInfo
1213       
1214        return(0)
1215End
1216
1217
1218// not quite the same as revert
1219Function DSM_ClearMaskProc(ctrlName) : ButtonControl
1220        String ctrlName
1221       
1222        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1223       
1224        Wave MaskData=$(USANSFolder+":DSM:MaskData")
1225        MaskData = NaN
1226       
1227        return(0)
1228end
1229
1230// when the mask button is pressed, A must be on the graph
1231// Displays MaskData wave on the graph
1232//
1233Function DSM_MyMaskProc(ctrlName) : ButtonControl
1234        String ctrlName
1235       
1236        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1237       
1238        Wave data=$(USANSFolder+":DSM:I_exp_orig")
1239       
1240//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1241        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1242       
1243// need to get rid of old smoothed data if data is re-masked
1244        Execute "RemoveSmoothed()"
1245        SetDataFolder $(USANSFolder+":DSM")
1246        Killwaves/Z I_smth,Q_smth,S_smth
1247               
1248        if(aExists)             //mask the selected point
1249                // toggle NaN (keep) or Data value (= masked)
1250                Wave MaskData
1251                MaskData[pcsr(A)] = (numType(MaskData[pcsr(A)])==0) ? NaN : data[pcsr(A)]               //if NaN, doesn't plot
1252        else
1253                Cursor /A=1/H=1/L=1/P/W=Desmear_Graph A I_exp_orig leftx(I_exp_orig)
1254                ShowInfo
1255                //if the mask wave does not exist, make one
1256                if(exists("MaskData") != 1)
1257                        Duplicate/O Q_exp_orig MaskData
1258                        MaskData = NaN          //use all data
1259                endif
1260                Execute "AppendMask()" 
1261        endif
1262
1263        SetDataFolder root:
1264
1265        return(0)
1266End
1267
1268// when the mask button is pressed, A must be on the graph
1269// Displays MaskData wave on the graph
1270//
1271Function DSM_MaskLTCursor(ctrlName) : ButtonControl
1272        String ctrlName
1273
1274        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1275               
1276//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1277        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1278       
1279        if(!aExists)
1280                return(1)
1281        endif
1282// need to get rid of old smoothed data if data is re-masked
1283        Execute "RemoveSmoothed()"
1284        SetDataFolder $(USANSFolder+":DSM")
1285        Killwaves/Z I_smth,Q_smth,S_smth
1286
1287        Wave data=I_exp_orig
1288
1289        Variable pt,ii
1290        pt = pcsr(A)
1291        for(ii=pt;ii>=0;ii-=1)
1292                // toggle NaN (keep) or Data value (= masked)
1293                Wave MaskData
1294                MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii]              //if NaN, doesn't plot
1295        endfor
1296
1297        SetDataFolder root:
1298        return(0)
1299End
1300
1301// when the mask button is pressed, A must be on the graph
1302// Displays MaskData wave on the graph
1303//
1304Function DSM_MaskGTCursor(ctrlName) : ButtonControl
1305        String ctrlName
1306       
1307        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1308
1309       
1310//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1311        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1312       
1313        if(!aExists)
1314                return(1)
1315        endif
1316// need to get rid of old smoothed data if data is re-masked
1317        Execute "RemoveSmoothed()"
1318        SetDataFolder root:DSM
1319        Killwaves/Z I_smth,Q_smth,S_smth
1320
1321        Wave data=I_exp_orig
1322
1323        Variable pt,ii,endPt
1324        endPt=numpnts(MaskData)
1325        pt = pcsr(A)
1326        for(ii=pt;ii<endPt;ii+=1)
1327                // toggle NaN (keep) or Data value (= masked)
1328                Wave MaskData
1329                MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii]              //if NaN, doesn't plot
1330        endfor
1331
1332        SetDataFolder root:
1333
1334        return(0)
1335End
1336
1337Function CleanUpJunk()
1338        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1339
1340        // clean up the old junk on the graph, /Z for no error
1341        Execute "RemoveExtrapolation()"
1342        Execute "RemoveDesmeared()"
1343        Execute "RemoveSmoothed()"
1344        Execute "RemoveMask()"
1345       
1346        //remove the cursor
1347        Cursor/K A
1348       
1349        //always initialize these
1350        String/G $(USANSFolder+":DSM:gStr1") = ""
1351        String/G $(USANSFolder+":DSM:gStr2") = ""
1352        String/G $(USANSFolder+":DSM:gIterStr") = ""
1353       
1354        // clean up the old waves from smoothing and desmearing steps
1355        SetDataFolder $(USANSFolder+":DSM")
1356        Killwaves/Z I_smth,I_dsm,I_dsm_sm,Q_smth,Q_dsm,S_smth,S_dsm,Yi_SS,Yq_SS
1357        Killwaves/Z Weights,FW,R_wave,S_wave,H_wave,Di,Ci
1358        Killwaves/Z Is_old,I_old,err
1359        Killwaves/Z Q_ext,I_ext,S_ext,hqq,hqi,lqq,lqi
1360        Killwaves/Z MaskData,Q_msk,I_msk,S_msk
1361        Killwaves/Z Q_work,I_work,S_work                                //working waves for desmearing step
1362        SetDataFolder root:
1363End
1364
1365// does not alter the data sets - reports a power law
1366// exponent and makes it global so it will automatically
1367// be used during the desmearing
1368//
1369// generates extr_hqi vs extr_hqq that are Appended to the graph
1370Function DSM_ExtrapolateButtonProc(ctrlName) : ButtonControl
1371        String ctrlName
1372
1373        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1374
1375        NVAR nend = $(USANSFolder+":DSM:gNptsExtrap")
1376        NVAR m_pred = $(USANSFolder+":DSM:gPowerM")
1377       
1378        SetDataFolder $(USANSFolder+":DSM")
1379//use masked data if it exists
1380        if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1381                wave Qw = $(USANSFolder+":DSM:Q_msk")
1382                Wave Iw = $(USANSFolder+":DSM:I_msk")
1383                Wave Sw = $(USANSFolder+":DSM:S_msk")
1384        else
1385                //start from the "_exp" waves
1386                if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1387                        wave Qw = $(USANSFolder+":DSM:Q_exp")
1388                        Wave Iw = $(USANSFolder+":DSM:I_exp")
1389                        Wave Sw = $(USANSFolder+":DSM:S_exp")
1390                endif
1391        endif
1392        SetDataFolder root:
1393       
1394        DSM_SetExtrWaves(Qw)
1395        m_pred = DSM_DoExtrapolate(Qw,Iw,Sw,nend)
1396        AppendExtrapolation()
1397       
1398        return(0)
1399End
1400
1401//smooths the data in steps as requested...
1402//
1403// typically do a simple Box smooth first,
1404// then do a smoothing spline, keeping the same number of data points
1405//
1406// chi-squared is reported - so you can see how "bad" the smoothing is
1407// smoothing of the data is done on a log(I) scale if requested
1408// setting doLog variable to 0 will return to linear smoothing
1409// (I see little difference)
1410//
1411// start from the "_smth" waves if they exist
1412// otheriwse start from the working waves
1413//
1414// create Q_smth, I_smth, S_smth
1415// keep track of the smoothing types and passes in the wave note
1416//
1417Function DSM_SmoothButtonProc(ctrlName) : ButtonControl
1418        String ctrlName
1419
1420        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1421
1422        SetDataFolder $(USANSFolder+":DSM")     
1423
1424        Variable ii,new_n,pass,nq_ext,offset,doLog=1
1425        String noteStr=""
1426       
1427        // want log scaling of intensity during smoothing?
1428        ControlInfo DSMControl_3g
1429        doLog = V_value
1430                       
1431////    if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1 && ! freshMask)
1432        if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1)
1433                //start from the smoothed waves --- just need the wave references
1434        else
1435                //start from the "msk", creating smth waves
1436                if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1437                        wave Q_msk,I_msk,S_msk
1438                        Duplicate/O I_msk,I_smth,Q_smth,S_smth
1439                        I_smth = I_msk
1440                        Q_smth = Q_msk
1441                        S_smth = S_msk
1442                else
1443                        //start from the "_exp" waves
1444                        if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1445                                wave Q_exp,I_exp,S_exp
1446                                Duplicate/O I_exp,I_smth,Q_smth,S_smth
1447                                I_smth = I_exp
1448                                Q_smth = Q_exp
1449                                S_smth = S_exp
1450                        endif
1451                endif
1452        endif
1453       
1454        wave Q_smth,I_smth,S_smth
1455
1456        // extend the data to avoid end effects
1457        //creates I_ext,Q_ext,S_ext with extended q-values
1458        // fit 15 pts at each end, typically add 10 pts to each end
1459        // does nothing if num_extr = 0
1460        ControlInfo/W=Desmear_Graph DSMControl_3d               //extendCheck
1461        if(V_value == 1)
1462                ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,10)
1463        else
1464                ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,0)            //don't extend, just use the original data
1465        endif
1466       
1467        //whether extending or not, the working data is "_ext", set by ExtendToSmooth()
1468        SetDataFolder $(USANSFolder+":DSM")     
1469        wave Q_ext,I_ext ,S_ext
1470       
1471        noteStr=note(I_smth)
1472        Note I_ext , noteStr
1473       
1474        WaveStats/Q I_ext
1475        if(V_min<=0)
1476                Print "negative itensity found, using linear scale for smoothing"
1477                doLog = 0
1478        endif
1479       
1480        if(doLog)
1481                //convert to log scale
1482                Duplicate/O I_ext I_log,I_log_err
1483                Wave I_log ,I_log_err
1484                I_log = log(I_ext)
1485                WaveStats/Q I_log
1486                offset = 2*(-V_min)             
1487                I_log += offset
1488                I_log_err = S_ext/(2.30*I_ext)
1489                I_ext = I_log
1490        endif   
1491        //
1492       
1493        ControlInfo/W=Desmear_Graph DSMControl_3b               //BoxCheck
1494        if(V_value == 1)
1495                //do a simple Box smooth first - number of points does not change
1496                // fills ends with neighboring value (E=3) (n=3 points in smoothing window)     
1497                Smooth/E=3/B 3, I_ext
1498               
1499                noteStr=note(I_ext)
1500                pass = NumberByKey("BOX", noteStr, "=", ";")
1501                noteStr = ReplaceNumberByKey("BOX", noteStr, pass+1, "=", ";")
1502                Note/K I_ext
1503                Note I_ext , noteStr
1504//              Print "box = ",noteStr
1505        endif
1506
1507        NVAR sParam = $(USANSFolder+":DSM:gSmoothFac")
1508       
1509        ControlInfo/W=Desmear_Graph DSMControl_3c               //SSCheck
1510        if(V_value == 1)
1511//              Interpolate2/T=3/N=(nq)/I=2/F=1/SWAV=S_ext/Y=Yi_SS/X=Yq_SS Q_ext, I_ext         //Igor 5
1512//              Interpolate/T=3/N=(nq)/I=2/F=1/S=S_ext/Y=Yi_SS/X=Yq_SS I_ext /X=Q_ext                   //Igor 4
1513        //Igor 4
1514                String str=""
1515                nq_ext = numpnts(Q_ext)
1516                str = "Interpolate/T=3/N=("+num2str(nq_ext)+")/I=1/F=("+num2str(sParam)+")/Y=Yi_SS/X=Yq_SS I_ext /X=Q_ext"             
1517                Execute str
1518//              Print Str
1519//              Interpolate2/T=3/N=(nq_ext)/I=1/F=(sParam)/Y=Yi_SS/X=Yq_SS Q_ext, I_ext
1520        // end Igor 4   
1521                wave yi_ss = root:DSM:yi_ss
1522                wave yq_ss = root:DSM:yq_ss
1523                // reassign the "working" waves with the result of interpolate, which has the same I/Q values
1524                I_ext = yi_ss
1525                Q_ext = yq_ss
1526               
1527                noteStr=note(I_ext)
1528                pass = NumberByKey("SPLINE", noteStr, "=", ";")
1529                noteStr = ReplaceNumberByKey("SPLINE", noteStr, pass+1, "=", ";")
1530                Note/K I_ext
1531                Note I_ext , noteStr
1532//              Print "spline = ",noteStr
1533        endif
1534       
1535        //undo the scaling
1536        If(doLog)
1537                I_ext -= offset
1538                I_ext = 10^I_ext
1539        endif
1540       
1541        // at this point, I_ext has too many points - and we need to just return the
1542        // center chunk that is the original q-values
1543        // as assign this to the working set for the desmear step
1544        // and to the _smth set in case another smoothng pass is desired
1545        // Q_smth has not changed, S_smth has not changed
1546        I_smth = interp(Q_smth[p], Q_ext, I_ext )
1547       
1548        Note/K I_smth
1549        Note I_smth , noteStr
1550//      Print "end of smoothed, note = ",note(I_smth)
1551       
1552        Variable nq = numpnts(root:DSM:Q_smth)
1553//      Print "nq after smoothing = ",nq
1554
1555        //reset the global
1556        NVAR gNq = $(USANSFolder+":DSM:gNq")
1557        gNq = nq
1558        //report the chi^2 difference between the smoothed curve and the experimental data
1559        NVAR chi2 = root:DSM:gChi2Smooth
1560        chi2 = SmoothedChi2(I_smth)
1561       
1562        Execute "AppendSmoothed()"
1563       
1564        setDataFolder root:
1565        return(0)
1566End
1567
1568// undo the smoothing and start over - useful if you've smoothed
1569// too aggressively and need to back off
1570Function DSM_SmoothUndoButtonProc(ctrlName) : ButtonControl
1571        String ctrlName
1572
1573        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1574       
1575        Execute "RemoveSmoothed()"
1576        SetDataFolder $(USANSFolder+":DSM")
1577        Killwaves/Z I_smth,Q_smth,S_smth,Q_ext,I_ext,S_ext,Yi_SS,Yq_SS
1578        SetDataFolder root:
1579        return (0)
1580end
1581
1582//calculate the chi^2 value for the smoothed data
1583Function SmoothedChi2(I_smth)
1584        Wave I_smth
1585       
1586        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1587       
1588        //start from the "msk", if they exist
1589        if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1590                Wave iw = $(USANSFolder+":DSM:I_msk")
1591                Wave sw = $(USANSFolder+":DSM:S_msk")
1592        else
1593                //start from the "_exp" waves
1594                if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1595                        Wave iw = $(USANSFolder+":DSM:I_exp")
1596                        Wave sw = $(USANSFolder+":DSM:S_exp")
1597                endif
1598        endif
1599       
1600        Variable ii,chi2,num
1601        chi2=0
1602        num = numpnts(iw)
1603        for(ii=0;ii<num;ii+=1)
1604                chi2 += (iw[ii] - I_smth[ii])^2 / (sw[ii])^2
1605        endfor
1606        Chi2 /= (num-1)
1607        return (chi2)
1608end
1609
1610// I_dsm is the desmeared data
1611//
1612// step (7) - desmearing is done, write out the result
1613//
1614Function DSM_SaveButtonProc(ctrlName) : ButtonControl
1615        String ctrlName
1616
1617        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1618
1619        String saveStr
1620        SVAR curFile = $(USANSFolder+":DSM:gCurFile")
1621        saveStr = CleanupName((curFile),0)              //the output filename
1622        //
1623        WriteUSANSDesmeared(saveStr,0,0,1)                      //use the full set (lo=hi=0) and present a dialog
1624       
1625        SetDataFolder root:
1626        return(0)
1627End
1628
1629Function DSM_HelpButtonProc(ctrlName) : ButtonControl
1630        String ctrlName
1631
1632        DisplayHelpTopic/Z/K=1 "Desmearing USANS Data"
1633        if(V_flag !=0)
1634                DoAlert 0,"The Desmearing USANS Data Help file could not be found"
1635        endif
1636        return(0)
1637End
1638
1639
1640//toggles the log/lin display of the loaded data set
1641Function DSM_LoadCheckProc(ctrlName,checked) : CheckBoxControl
1642        String ctrlName
1643        Variable checked
1644
1645        ModifyGraph log=(checked)
1646        return(0)
1647End
1648
1649
1650Function DSM_ExtrapolationCheckProc(ctrlName,checked) : CheckBoxControl
1651        String ctrlName
1652        Variable checked
1653
1654        if(checked)
1655                AppendExtrapolation()
1656        else
1657                RemoveExtrapolation()
1658        endif
1659        return(0)
1660End
1661
1662
1663// takes as input the "working waves"
1664//
1665// creates intermediate waves to work on
1666//
1667// output of Q_dsm,I_dsm,S_dsm (and I_dsm_sm, the result of smearing I_dsm)
1668//
1669Function DSM_DesmearButtonProc(ctrlName) : ButtonControl
1670        String ctrlName
1671       
1672        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1673       
1674        SetDataFolder $(USANSFolder+":DSM")
1675        if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1)
1676                wave Q_smth , I_smth ,S_smth
1677                Duplicate/O I_smth,I_work,Q_work,S_work
1678                I_work = I_smth
1679                Q_work = Q_smth
1680                S_work = S_smth
1681        else
1682                //start from the "msk", creating work waves
1683                if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1684                        wave Q_msk,I_msk,S_msk
1685                        Duplicate/O I_msk,I_work,Q_work,S_work
1686                        I_work = I_msk
1687                        Q_work = Q_msk
1688                        S_work = S_msk
1689                else
1690                        //start from the "_exp" waves
1691                        if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1692                                wave Q_exp,I_exp,S_exp
1693                                Duplicate/O I_exp,I_work,Q_work,S_work
1694                                I_work = I_exp
1695                                Q_work = Q_exp
1696                                S_work = S_exp
1697                        endif
1698                endif
1699        endif
1700        //SetDataFolder root:
1701       
1702        NVAR nq = gNq
1703        NVAR m = gPowerM
1704        NVAR chi2_target = gChi2Target
1705        NVAR maxFastIter = gMaxFastIter
1706        NVAR maxSlowIter = gMaxSlowIter
1707       
1708        //      SET WEIGHTING OF EXPERIMENTAL DATA.
1709        Duplicate/O Q_work weights
1710        Wave weights = weights
1711        weights = 1/S_work^2
1712
1713//      calculate weighting array for smearing of data
1714        Make/O/D/N=(nq,nq) FW
1715        Wave FW
1716        Weights_L(m,FW,Q_work)
1717
1718//      ^^^^   Iterative desmearing   ^^^^*
1719        Variable chi2_old,chi2_new,done,iter
1720//      FOR 0TH ITERATION, EXPERIMENTAL DATA IS USED FOR Y_OLD, create ys_old for result
1721//      y_old = I_old, y_new = I_dsm, I_dsm_sm = ys_new,
1722// duplicate preserves the wave note!
1723        Duplicate/O I_work I_old,Is_old,root:DSM:I_dsm,root:DSM:I_dsm_sm
1724        Duplicate/O Q_work Q_dsm,S_dsm          //sets Q_dsm correctly
1725        wave S_dsm,I_old,Is_old,I_dsm,I_dsm_sm
1726        I_old = I_work
1727        Is_old = 0
1728        I_dsm = 0
1729        I_dsm_sm = 0
1730       
1731//      Smear 0TH iter guess Y --> YS
1732        DSM_Smear(I_old,Is_old,NQ,FW)
1733        chi2_OLD = chi2(Is_old,I_work,weights,NQ)
1734       
1735        printf "starting chi^2 = %g\r",chi2_old
1736        Print "Starting fast convergence..... "
1737
1738        done = 0                //false
1739        iter = 0
1740        Do                                      // while not done, use Fast convergence
1741                I_dsm = I_work * I_old / Is_old //        Calculate corrected guess (no need for do-loop)
1742
1743//        Smear iter guess I_dsm --> I_dsm_sm and see how well I_dsm_sm matches experimental data
1744                DSM_Smear(I_dsm,I_dsm_sm,NQ,FW)
1745                chi2_new = chi2(I_dsm_sm,I_work,weights,NQ)
1746
1747//        Stop iteration if fit from new iteration has worse fit...
1748                If (chi2_new > chi2_old)
1749                        Done = 1
1750                Endif
1751               
1752//        Stop iteration if fit is better than target value...
1753                If (chi2_new < chi2_target)
1754                        Done = 1
1755                Else
1756                        Iter += 1
1757                        Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1758                        I_old = I_dsm
1759                        Is_old = I_dsm_sm
1760                        if(iter>maxFastIter)
1761                                break
1762                        endif
1763                        CHI2_OLD = CHI2_NEW
1764                EndIf
1765        while( !done )
1766       
1767        // append I_dsm,I_exp to the graph
1768        Execute "AppendDesmeared()"
1769        DoUpdate
1770        SetDataFolder $(USANSFolder+":DSM")
1771
1772       
1773// step (6) - refine the desmearing using slow convergence
1774        Print "Starting slow convergence..... "
1775        done = 0                //  ! reset flag for slow convergence
1776        Do
1777                I_dsm = I_old + (I_work - Is_old)       //        Calculate corrected guess
1778
1779//        Smear iter guess Y --> YS
1780                DSM_Smear(I_dsm,I_dsm_sm,NQ,FW)
1781                chi2_new = chi2(I_dsm_sm,I_work,weights,NQ)
1782
1783//        Stop iteration if fit from new iteration has worse fit...
1784                If (chi2_new > chi2_old)
1785                        Done = 1
1786                Endif
1787               
1788//        Stop iteration if fit is better than target value...
1789                If (chi2_new < chi2_target)
1790                        Done = 1
1791                Else
1792                        Iter += 1
1793                       
1794                        if(mod(iter, 50 ) ==0 )
1795                                Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1796                                DoUpdate
1797                        endif
1798                        I_old = I_dsm
1799                        Is_old = I_dsm_sm
1800
1801                        CHI2_OLD = CHI2_NEW
1802                EndIf
1803                if(iter>maxSlowIter)
1804                        break
1805                endif
1806        While ( !done )
1807
1808        Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1809
1810        // adjust the error
1811        SetDataFolder $(USANSFolder+":DSM")
1812        Duplicate/O S_work err
1813        S_dsm = abs(err/I_work*I_dsm)           //proportional error
1814//John simply keeps the same error as was read in from the smeared data set - is this right?
1815//      S_dsm = S_Work
1816       
1817        NVAR gChi =  gChi2Final         //chi^2 final
1818        NVAR gIter = gIterations                //total number of iterations
1819        gChi = chi2_new
1820        gIter = Iter
1821       
1822        SVAR str = gIterStr
1823        sprintf str,"%d iterations required",iter
1824       
1825        SetDataFolder root:
1826        return(0)
1827End
1828
1829Function DSM_RevertButtonProc(ctrlName) : ButtonControl
1830        String ctrlName
1831       
1832        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1833       
1834       
1835        CleanUpJunk()
1836
1837        SetDataFolder $(USANSFolder+":DSM")
1838
1839        //reset the working waves to the original
1840        wave Q_exp_orig,I_exp_orig,S_exp_orig
1841
1842        Duplicate/O Q_exp_orig Q_exp
1843        Duplicate/O I_exp_orig I_exp
1844        Duplicate/O S_exp_orig S_exp
1845        // kill the data folder
1846        // re-initialize?
1847
1848        SetDataFolder root:
1849       
1850        return(0)
1851End
Note: See TracBrowser for help on using the repository browser.