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

Last change on this file since 693 was 665, checked in by srkline, 13 years ago

Made preferences a common panel (moved to PlotUtilsMacro?.ipf and globals to root:Packages:NIST:) and added menu items for all packages. Many files had to be modified so that the preferences could be properly accessed

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