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

Last change on this file since 477 was 477, checked in by srkline, 14 years ago

Added save/restore buttons to SimpleGlobal? fit to be able to restore the state of the checkboxes if possible. Checkboxes are reset when new data sets are selected.

Typos in SANSRedn help file have been corrected.

Error bar styles have been changed to /T=0, a vertical line with no horizontal bar. My preference - I got tired of seeing more of the horizontal bar than the actual data.

  • Property eol-style set to native
File size: 52.6 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// 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,dumWave                //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        //actually open the file
594        Open refNum as fullpath
595       
596        fprintf refnum,"%s"+termStr,samStr
597        fprintf refnum,"%s"+termStr,str1
598        fprintf refnum,"%s"+termStr,str2
599        fprintf refnum,"%s"+termStr,dateStr
600       
601        wfprintf refnum, formatStr, tq,ti,te,res1,res2,res3
602       
603        Close refnum
604       
605        Killwaves/Z ti,tq,te,res1,res2,res3
606       
607        Return(0)
608End
609
610/// procedures to do the extrapolation to high Q
611// very similar to the procedures in the Invariant package
612//
613//create the wave to extrapolate
614// w is the input q-values
615Function DSM_SetExtrWaves(w)
616        Wave w
617
618        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
619
620        Variable num_extr=25
621       
622        SetDataFolder $(USANSFolder+":DSM")
623       
624        Make/O/D/N=(num_extr) extr_hqq,extr_hqi
625        extr_hqi=1              //default values
626
627        //set the q-range
628        Variable qmax,num
629//      qmax=0.03
630
631        num=numpnts(w)
632        qmax=6*w[num-1]
633       
634        extr_hqq = w[num-1] + x * (qmax-w[num-1])/num_extr
635       
636        SetDataFolder root:
637        return(0)
638End
639
640//creates I_ext,Q_ext,S_ext with extended q-values
641//
642// if num_extr == 0 , the input waves are returned as _ext
643//
644// !! uses simple linear extrapolation at low Q - just need something
645// reasonable in the waves to keep from getting smoothing artifacts
646// - uses power law at high q
647Function ExtendToSmooth(qw,iw,sw,nbeg,nend,num_extr)
648        Wave qw,iw,sw
649        Variable nbeg,nend,num_extr
650
651        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
652       
653        Variable/G V_FitMaxIters=300
654        Variable/G V_fitOptions=4               //suppress the iteration window
655        Variable num_new,qmax,num,qmin
656        num=numpnts(qw)
657       
658        num_new = num + 2*num_extr
659       
660//      Print "num,num_new",num,num_new
661       
662        SetDataFolder $(USANSFolder+":DSM")
663        Make/O/D/N=(num_new) Q_ext,I_ext,S_ext
664       
665        if(num_extr == 0)               //the extended waves are the input, get out
666                Q_ext = qw
667                I_ext = iw
668                S_ext = sw
669                setDatafolder root:
670                return (0)
671        endif
672       
673        //make the extensions
674        Make/O/D/N=(num_extr) hqq,hqi,lqq,lqi
675        hqi=1           //default values
676        lqi=0
677        //set the q-range based on the high/low values of the data set
678//      qmin=1e-6
679        qmin= qw[0]*0.5
680        qmax = qw[num-1]*2
681       
682        hqq = qw[num-1] + x * (qmax-qw[num-1])/num_extr
683        lqq = qmin + x * (qw[0]-qmin)/num_extr
684       
685        //do the fits
686        Duplicate/O iw dummy                    //use this as the destination to suppress fit_ from being appended
687       
688        //Use simple linear fits        line: y = K0+K1*x
689        CurveFit/Q line iw[0,(nbeg-1)] /X=qw /W=sw /D=dummy
690        Wave W_coef=W_coef
691        lqi=W_coef[0]+W_coef[1]*lqq
692//     
693// Guinier or Power-law fits
694//             
695//      Make/O/D G_coef={100,-100}              //input
696//      FuncFit DSM_Guinier_Fit G_coef iw[0,(nbeg-1)] /X=qw /W=sw /D
697//      lqi= DSM_Guinier_Fit(G_coef,lqq)
698       
699//      Printf "I(q=0) = %g (1/cm)\r",G_coef[0]
700//      Printf "Rg = %g (A)\r",sqrt(-3*G_coef[1])
701               
702        Make/O/D P_coef={0,1,-4}                        //input  --- (set background to zero and hold fixed)
703        CurveFit/Q/N/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /W=sw /D=dummy
704        hqi=P_coef[0]+P_coef[1]*hqq^P_coef[2]
705       
706//      Printf "Power law exponent = %g\r",P_coef[2]
707//      Printf "Pre-exponential = %g\r",P_coef[1]
708       
709        // concatenate the extensions
710        Q_ext[0,(num_extr-1)] = lqq[p]
711        Q_ext[num_extr,(num_extr+num-1)] = qw[p-num_extr]
712        Q_ext[(num_extr+num),(2*num_extr+num-1)] = hqq[p-num_extr-num]
713        I_ext[0,(num_extr-1)] = lqi[p]
714        I_ext[num_extr,(num_extr+num-1)] = iw[p-num_extr]
715        I_ext[(num_extr+num),(2*num_extr+num-1)] = hqi[p-num_extr-num]
716        S_ext[0,(num_extr-1)] = sw[0]
717        S_ext[num_extr,(num_extr+num-1)] = sw[p-num_extr]
718        S_ext[(num_extr+num),(2*num_extr+num-1)] = sw[num-1]
719       
720        killwaves/z dummy
721        SetDataFolder root:
722        return(0)
723End
724
725// pass in the (smeared) experimental data and find the power-law extrapolation
726// 10 points is usually good enough, unless the data is crummy
727// returns the prediction for the exponent
728//
729Function DSM_DoExtrapolate(qw,iw,sw,nend)
730        Wave qw,iw,sw
731        Variable nend
732
733        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
734       
735        Setdatafolder $(USANSFolder+":DSM")
736       
737//      Wave extr_lqi=extr_lqi
738//      Wave extr_lqq=extr_lqq
739        Wave extr_hqi=extr_hqi
740        Wave extr_hqq=extr_hqq
741        Variable/G V_FitMaxIters=300
742        Variable/G V_fitOptions=4               //suppress the iteration window
743        Variable num=numpnts(iw),retVal
744       
745        Make/O/D P_coef={0,1,-4}                        //input
746//      Make/O/T Constr={"K2<0","K2 > -20"}
747        //(set background to zero and hold fixed)
748        CurveFit/H="100" 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        // for the case of data with a background
752//      Make/O/D P_coef={0,1,-4}                        //input
753//      //(set background to iw[num-1], let it be a free parameter
754//      P_coef[0] = iw[num-1]
755//      CurveFit Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /W=sw /D
756//      extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
757       
758       
759        Printf "Smeared Power law exponent = %g\r",P_coef[2]
760        Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
761       
762        retVal = P_coef[2]-1
763        SetDataFolder root:
764        return(retVal)
765End
766
767Function DSM_Guinier_Fit(w,x) //: FitFunc
768        Wave w
769        Variable x
770       
771        //fit data to I(q) = A*exp(B*q^2)
772        // (B will be negative)
773        //two parameters
774        Variable a,b,ans
775        a=w[0]
776        b=w[1]
777        ans = a*exp(b*x*x)
778        return(ans)
779End
780
781Proc Desmear_Graph()
782        String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
783
784        PauseUpdate; Silent 1           // building window...
785        Display /W=(5,44,408,558) /K=1
786        ModifyGraph cbRGB=(51664,44236,58982)
787        DoWindow/C Desmear_Graph
788        ControlBar 160
789        // break into tabs
790        TabControl DSM_Tab,pos={5,3},size={392,128},proc=DSM_TabProc
791        TabControl DSM_Tab,labelBack=(49151,49152,65535),tabLabel(0)="Load"
792        TabControl DSM_Tab,tabLabel(1)="Mask",tabLabel(2)="Extrapolate"
793        TabControl DSM_Tab,tabLabel(3)="Smooth",tabLabel(4)="Desmear",value= 0
794       
795        //always visible - revert and save
796        //maybe the wrong place here?
797        Button DSMControlA,pos={225,135},size={80,20},proc=DSM_RevertButtonProc,title="Revert"
798        Button DSMControlA,help={"Revert the smeared data to its original state and start over"}
799        Button DSMControlB,pos={325,135},size={50,20},proc=DSM_SaveButtonProc,title="Save"
800        Button DSMControlB,help={"Save the desmeared data set"}
801        Button DSMControlC,pos={25,135},size={50,20},proc=DSM_HelpButtonProc,title="Help"
802        Button DSMControlC,help={"Show the help file for desmearing"}
803       
804        // add the controls to each tab ---- all names start with "DSMControl_"
805
806        //tab(0) Load - initially visible
807        Button DSMControl_0a,pos={23,39},size={80,20},proc=DSM_LoadButtonProc,title="Load Data"
808        Button DSMControl_0a,help={"Load slit-smeared USANS data = \".cor\" files"}
809        CheckBox DSMControl_0b,pos={26,74},size={80,14},proc=DSM_LoadCheckProc,title="Log Axes?"
810        CheckBox DSMControl_0b,help={"Toggle Log/Lin Q display"},value= 1
811        TitleBox DSMControl_0c,pos={120,37},size={104,19},font="Courier",fSize=10
812        TitleBox DSMControl_0c,variable= $(USANSFolder+":DSM:gStr1")
813        //second message string not used currently
814//      TitleBox DSMControl_0d,pos={120,57},size={104,19},font="Courier",fSize=10
815//      TitleBox DSMControl_0d,variable= root:DSM:gStr2
816       
817        //tab(1) Mask
818        Button DSMControl_1a,pos={20,35},size={90,20},proc=DSM_MyMaskProc,title="Mask Point"            //bMask
819        Button DSMControl_1a,help={"Toggles the masking of the selected data point"}
820        Button DSMControl_1a,disable=1
821        Button DSMControl_1b,pos={20,65},size={140,20},proc=DSM_MaskGTCursor,title="Mask Q >= Cursor"           //bMask
822        Button DSMControl_1b,help={"Toggles the masking of all q-values GREATER than the current cursor location"}
823        Button DSMControl_1b,disable=1
824        Button DSMControl_1c,pos={20,95},size={140,20},proc=DSM_MaskLTCursor,title="Mask Q <= Cursor"           //bMask
825        Button DSMControl_1c,help={"Toggles the masking of all q-values LESS than the current cursor location"}
826        Button DSMControl_1c,disable=1
827        Button DSMControl_1d,pos={180,35},size={90,20},proc=DSM_ClearMaskProc,title="Clear Mask"                //bMask
828        Button DSMControl_1d,help={"Clears all mask points"}
829        Button DSMControl_1d,disable=1
830//      Button DSMControl_1b,pos={144,66},size={110,20},proc=DSM_MaskDoneButton,title="Done Masking"
831//      Button DSMControl_1b,disable=1
832       
833        //tab(2) Extrapolate
834        Button DSMControl_2a,pos={31,42},size={90,20},proc=DSM_ExtrapolateButtonProc,title="Extrapolate"
835        Button DSMControl_2a,help={"Extrapolate the high-q region with a power-law"}
836        Button DSMControl_2a,disable=1
837        SetVariable DSMControl_2b,pos={31,70},size={100,15},title="# of points"
838        SetVariable DSMControl_2b,help={"Set the number of points for the power-law extrapolation"}
839        SetVariable DSMControl_2b,limits={5,100,1},value=  $(USANSFolder+":DSM:gNptsExtrap")
840        SetVariable DSMControl_2b,disable=1
841        CheckBox DSMControl_2c,pos={157,45},size={105,14},proc=DSM_ExtrapolationCheckProc,title="Show Extrapolation"
842        CheckBox DSMControl_2c,help={"Show or hide the high q extrapolation"},value= 1
843        CheckBox DSMControl_2c,disable=1
844        SetVariable DSMControl_2d,pos={31,96},size={150,15},title="Power Law Exponent"
845        SetVariable DSMControl_2d,help={"Power Law exponent from the fit = the DESMEARED slope - override as needed"}
846        SetVariable DSMControl_2d format="%5.2f"
847        SetVariable DSMControl_2d,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gPowerM")
848        SetVariable DSMControl_2d,disable=1
849       
850        //tab(3) Smooth
851        Button DSMControl_3a,pos={34,97},size={70,20},proc=DSM_SmoothButtonProc,title="Smooth"
852        Button DSMControl_3a,disable=1
853                //BoxCheck
854        CheckBox DSMControl_3b,pos={34,39},size={35,14},title="Box",value= 1
855        CheckBox DSMControl_3b,help={"Use a single pass of 3-point box smoothing"}
856        CheckBox DSMControl_3b,disable=1
857                //SSCheck
858        CheckBox DSMControl_3c,pos={34,60},size={45,14},title="Spline",value= 0
859        CheckBox DSMControl_3c,help={"Use a single pass of a smoothing spline"}
860        CheckBox DSMControl_3c,disable=1
861                //extendCheck
862        CheckBox DSMControl_3d,pos={268,60},size={71,14},title="Extend Data"
863        CheckBox DSMControl_3d,help={"extends the data at both low q and high q to avoid end effects in smoothing"}
864        CheckBox DSMControl_3d,value= 0
865        CheckBox DSMControl_3d,disable=1
866        Button DSMControl_3e,pos={125,97},size={90,20},proc=DSM_SmoothUndoButtonProc,title="Start Over"
867        Button DSMControl_3e,help={"Start the smoothing over again without needing to re-mask the data set"}
868        Button DSMControl_3e,disable=1
869        SetVariable DSMControl_3f,pos={94,60},size={150,15},title="Smoothing factor"
870        SetVariable DSMControl_3f,help={"Smoothing factor for the smoothing spline"}
871        SetVariable DSMControl_3f format="%5.4f"
872        SetVariable DSMControl_3f,limits={0.01,2,0.01},value= $(USANSFolder+":DSM:gSmoothFac")
873        SetVariable DSMControl_3f,disable=1
874        CheckBox DSMControl_3g,pos={268,39},size={90,14},title="Log-scale smoothing?"
875        CheckBox DSMControl_3g,help={"Use log-scaled intensity during smoothing (reverts to linear if negative intensity points found)"}
876        CheckBox DSMControl_3g,value=0
877        CheckBox DSMControl_3g,disable=1
878        ValDisplay DSMControl_3h pos={235,97},title="Chi^2",size={80,20},value=root:Packages:NIST:USANS:DSM:gChi2Smooth
879        ValDisplay DSMControl_3h,help={"This is the Chi^2 value for the smoothed data vs experimental data"}
880        ValDisplay DSMControl_3h,disable=1
881       
882        //tab(4) Desmear
883        Button DSMControl_4a,pos={35,93},size={80,20},proc=DSM_DesmearButtonProc,title="Desmear"
884        Button DSMControl_4a,help={"Do the desmearing - the result is in I_dsm"}
885        Button DSMControl_4a,disable=1
886        SetVariable DSMControl_4b,pos={35,63},size={120,15},title="Chi^2 target"
887        SetVariable DSMControl_4b,help={"Set the targetchi^2 for convergence (recommend chi^2=1)"}
888        SetVariable DSMControl_4b,limits={0,inf,0.1},value= $(USANSFolder+":DSM:gChi2Target")
889        SetVariable DSMControl_4b,disable=1
890        SetVariable DSMControl_4c,pos={35,35},size={80,15},title="dQv"
891        SetVariable DSMControl_4c,help={"Slit height as read in from the data file. 0.117 is the NIST value, override if necessary"}
892        SetVariable DSMControl_4c,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gDqv")
893        SetVariable DSMControl_4c,disable=1
894        TitleBox DSMControl_4d,pos={160,37},size={104,19},font="Courier",fSize=10
895        TitleBox DSMControl_4d,variable= $(USANSFolder+":DSM:gIterStr")
896        TitleBox DSMControl_4d,disable=1
897       
898       
899        SetDataFolder root:
900EndMacro
901
902// function to control the drawing of buttons in the TabControl on the main panel
903// Naming scheme for the buttons MUST be strictly adhered to... else buttons will
904// appear in odd places...
905// all buttons are named DSMControl_NA where N is the tab number and A is the letter denoting
906// the button's position on that particular tab.
907// in this way, buttons will always be drawn correctly :-)
908//
909Function DSM_TabProc(ctrlName,tab) //: TabControl
910        String ctrlName
911        Variable tab
912
913        String ctrlList = ControlNameList("",";"),item="",nameStr=""
914        Variable num = ItemsinList(ctrlList,";"),ii,onTab
915        for(ii=0;ii<num;ii+=1)
916                //items all start w/"DSMControl_"               //11 characters
917                item=StringFromList(ii, ctrlList ,";")
918                nameStr=item[0,10]
919                if(cmpstr(nameStr,"DSMControl_")==0)
920                        onTab = str2num(item[11])                       //12th is a number
921                        ControlInfo $item
922                        switch(abs(V_flag))     
923                                case 1:
924                                        Button $item,disable=(tab!=onTab)
925                                        break
926                                case 2:
927                                        CheckBox $item,disable=(tab!=onTab)
928                                        break
929                                case 5:
930                                        SetVariable $item,disable=(tab!=onTab)
931                                        break
932                                case 10:       
933                                        TitleBox $item,disable=(tab!=onTab)
934                                        break
935                                case 4:
936                                        ValDisplay $item,disable=(tab!=onTab)
937                                        break
938                                // add more items to the switch if different control types are used
939                        endswitch
940                       
941                endif
942        endfor
943       
944        if(tab==1)
945                DSM_MyMaskProc("")              //start maksing if you click on the tab
946        else
947                DSM_MaskDoneButton("")          //masking is done if you click off the tab
948        endif
949       
950        return 0
951End
952
953Proc AppendSmeared()
954        String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
955
956        SetDataFolder $(USANSFolder+":DSM")
957//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0,2) == -1)            //Igor 5
958        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0) == -1)
959                AppendToGraph/W=Desmear_Graph  I_exp_orig vs Q_exp_orig
960                ModifyGraph mode=4,marker=19            //3 is markers, 4 is markers and lines
961                ModifyGraph rgb(I_exp_orig)=(0,0,0)
962                ModifyGraph msize=2,grid=1,log=1,mirror=2,standoff=0,tickunit=1
963                ErrorBars/T=0 I_exp_orig Y,wave=(S_exp_orig,S_exp_orig)
964                Legend/N=text0/J "\\F'Courier'\\s(I_exp_orig) I_exp_orig"
965                Label left "Intensity"
966                Label bottom "Q (1/A)"
967        endif
968        //always update the textbox - kill the old one first
969        TextBox/K/N=text1
970//      TextBox/C/N=text1/F=0/A=MT/E=2/X=5.50/Y=0.00 root:DSM:gCurFile                  //Igor 5
971        TextBox/C/N=text1/F=0/A=MT/E=1/X=5.50/Y=0.00 $(USANSFolder+":DSM:gCurFile")
972End
973
974Proc AppendMask()
975        String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
976
977//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0,2) == -1)                      //Igor 5
978        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0) == -1)
979                SetDataFolder $(USANSFolder+":DSM:")
980                AppendToGraph/W=Desmear_Graph MaskData vs Q_exp_orig
981                ModifyGraph mode(MaskData)=3,marker(MaskData)=8,msize(MaskData)=2.5,opaque(MaskData)=1
982                ModifyGraph rgb(MaskData)=(65535,16385,16385)
983                setdatafolder root:
984        endif
985end
986
987Proc AppendSmoothed()
988        String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
989
990//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0,2) == -1)                        //Igor 5
991        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0) == -1)
992                SetDataFolder $(USANSFolder+":DSM:")
993                AppendToGraph/W=Desmear_Graph I_smth vs Q_smth
994                ModifyGraph/W=Desmear_Graph rgb(I_smth)=(3,52428,1),lsize(I_smth)=2
995                setdatafolder root:
996        endif
997end
998
999Function RemoveSmoothed()
1000        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1001
1002        SetDataFolder $(USANSFolder+":DSM:")
1003        RemoveFromGraph/W=Desmear_Graph/Z I_smth
1004        setdatafolder root:
1005end
1006
1007Function RemoveMask()
1008        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1009
1010        SetDataFolder $(USANSFolder+":DSM:")
1011        RemoveFromGraph/W=Desmear_Graph/Z MaskData
1012        setdatafolder root:
1013end
1014
1015Proc AppendDesmeared()
1016        String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1017
1018//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0,2) == -1)         //Igor 5
1019        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0) == -1)
1020                SetDataFolder $(USANSFolder+":DSM:")
1021                AppendToGraph/W=Desmear_Graph I_dsm vs Q_dsm
1022                ModifyGraph mode(I_dsm)=3,marker(I_dsm)=19
1023                ModifyGraph rgb(I_dsm)=(1,16019,65535),msize(I_dsm)=2
1024                ErrorBars/T=0 I_dsm Y,wave=(S_dsm,S_dsm)
1025                setdatafolder root:
1026        endif
1027end
1028
1029Function RemoveDesmeared()
1030        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1031
1032        SetDataFolder $(USANSFolder+":DSM:")
1033        RemoveFromGraph/W=Desmear_Graph/Z I_dsm
1034        setdatafolder root:
1035end
1036
1037Function AppendExtrapolation()
1038        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1039
1040//      if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0,2) == -1)              //Igor 5
1041        if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0) == -1)
1042                SetDataFolder $(USANSFolder+":DSM:")
1043                AppendToGraph/W=Desmear_Graph extr_hqi vs extr_hqq
1044                ModifyGraph/W=Desmear_Graph lSize(extr_hqi)=2
1045                setdatafolder root:
1046        endif
1047end
1048
1049Function RemoveExtrapolation()
1050        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1051
1052        SetDataFolder $(USANSFolder+":DSM:")
1053        RemoveFromGraph/W=Desmear_Graph/Z extr_hqi
1054        setdatafolder root:
1055end
1056
1057// step (1) - read in the data, and plot it
1058// clear out all of the "old" waves, remove them from the graph first
1059// reads in a fresh copy of the data
1060//
1061// produces Q_exp, I_exp, S_exp waves (and originals "_orig")
1062// add a dummy wave note that can be changed on later steps
1063//
1064Function DSM_LoadButtonProc(ctrlName) : ButtonControl
1065        String ctrlName
1066
1067        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1068
1069        String qStr,iStr,sStr,sqStr
1070        Variable nq,dqv,numBad,val
1071       
1072        // remove any of the old traces on the graph and delete the waves
1073        CleanUpJunk()
1074       
1075        SetDataFolder root:
1076        Execute "A_LoadOneDDataWithName(\"\",0)"
1077        //define the waves that the smoothing will be looking for...
1078        SVAR fname = $("root:Packages:NIST:gLastFileName")              //this changes as any data is loaded
1079        SVAR curFile = $(USANSFolder+":DSM:gCurFile")                                   //keep this for title, save
1080        curFile = fname
1081       
1082        qStr = CleanupName((fName + "_q"),0)            //the q-wave
1083        iStr = CleanupName((fName + "_i"),0)            //the i-wave
1084        sStr = CleanupName((fName + "_s"),0)            //the s-wave
1085        sqStr = CleanupName((fName + "sq"),0)           //the sq-wave, which should have -dQv as the elements
1086
1087        String DFStr= CleanupName(fname,0)
1088       
1089        Duplicate/O $("root:"+DFStr+":"+qStr) $(USANSFolder+":DSM:Q_exp ")             
1090        Duplicate/O $("root:"+DFStr+":"+iStr) $(USANSFolder+":DSM:I_exp")               
1091        Duplicate/O $("root:"+DFStr+":"+sStr) $(USANSFolder+":DSM:S_exp ")     
1092        wave Q_exp = $(USANSFolder+":DSM:Q_exp")
1093        Wave I_exp = $(USANSFolder+":DSM:I_exp")
1094        Wave S_exp = $(USANSFolder+":DSM:S_exp")
1095       
1096        // remove any negative q-values (and q=0 values!)(and report this)
1097        // ? and trim the low q to be >= 3.0e-5 (1/A), below this USANS is not reliable.
1098        NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,0)
1099        SVAR str1 = $(USANSFolder+":DSM:gStr1")
1100        sprintf str1,"%d negative q-values were removed",numBad
1101       
1102// don't trim off any positive q-values
1103//      val = 3e-5              //lowest "good" q-value from USANS
1104//      NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,val-1e-8)
1105//      SVAR str2 = root:DSM:gStr2
1106//      sprintf str2,"%d q-values below q = %g were removed",numBad,val
1107       
1108        Duplicate/O $(USANSFolder+":DSM:Q_exp") $(USANSFolder+":DSM:Q_exp_orig")
1109        Duplicate/O $(USANSFolder+":DSM:I_exp") $(USANSFolder+":DSM:I_exp_orig")
1110        Duplicate/O $(USANSFolder+":DSM:S_exp") $(USANSFolder+":DSM:S_exp_orig")
1111        wave I_exp_orig = $(USANSFolder+":DSM:I_exp_orig")
1112       
1113        nq = numpnts($(USANSFolder+":DSM:Q_exp"))
1114       
1115        dQv = NumVarOrDefault("root:"+DFStr+":USANS_dQv", 0.117 )
1116//      if(WaveExists(sigQ))                    //try to read dQv
1117////            dqv = -sigQ[0][0]
1118////            DoAlert 0,"Found dQv value of " + num2str(dqv)
1119//      else
1120//              dqv = 0.117
1121//      //      dqv = 0.037             //old value
1122//              DoAlert 0,"Could not find dQv in the data file - using " + num2str(dqv)
1123//      endif
1124        NVAR gDqv = $(USANSFolder+":DSM:gDqv")                          //needs to be global for Weights_L()
1125        NVAR gNq = $(USANSFolder+":DSM:gNq")
1126        //reset the globals
1127        gDqv = dqv
1128        gNq = nq
1129       
1130        // append the (blank) wave note to the intensity wave
1131        Note I_exp,"BOX=0;SPLINE=0;"
1132        Note I_exp_orig,"BOX=0;SPLINE=0;"
1133       
1134        //draw the graph
1135        Execute "AppendSmeared()"
1136       
1137        SetDataFolder root:
1138End
1139
1140// remove any q-values <= val
1141Function RemoveBadQPoints(qw,iw,sw,val)
1142        Wave qw,iw,sw
1143        Variable val
1144       
1145        Variable ii,num,numBad,qval
1146        num = numpnts(qw)
1147       
1148        ii=0
1149        numBad=0
1150        do
1151                qval = qw[ii]
1152                if(qval <= val)
1153                        numBad += 1
1154                else            //keep the points
1155                        qw[ii-numBad] = qval
1156                        iw[ii-numBad] = iw[ii]
1157                        sw[ii-numBad] = sw[ii]
1158                endif
1159                ii += 1
1160        while(ii<num)
1161        //trim the end of the waves
1162        DeletePoints num-numBad, numBad, qw,iw,sw
1163        return(numBad)
1164end
1165
1166// if mw = Nan, keep the point, if a numerical value, delete it
1167Function RemoveMaskedPoints(mw,qw,iw,sw)
1168        Wave mw,qw,iw,sw
1169       
1170        Variable ii,num,numBad,mask
1171        num = numpnts(qw)
1172       
1173        ii=0
1174        numBad=0
1175        do
1176                mask = mw[ii]
1177                if(numtype(mask) != 2)          //if not NaN
1178                        numBad += 1
1179                else            //keep the points that are NaN
1180                        qw[ii-numBad] = qw[ii]
1181                        iw[ii-numBad] = iw[ii]
1182                        sw[ii-numBad] = sw[ii]
1183                endif
1184                ii += 1
1185        while(ii<num)
1186        //trim the end of the waves
1187        DeletePoints num-numBad, numBad, qw,iw,sw
1188        return(numBad)
1189end
1190
1191// produces the _msk waves that have the new number of data points
1192//
1193Function DSM_MaskDoneButton(ctrlName) : ButtonControl
1194        String ctrlName
1195
1196        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1197
1198//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1199        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1200        if(!aExists)
1201                return(1)               //possibly reverted data, no cursor, no Mask wave
1202        endif
1203       
1204        Duplicate/O $(USANSFolder+":DSM:Q_exp_orig"),$(USANSFolder+":DSM:Q_msk")
1205        Duplicate/O $(USANSFolder+":DSM:I_exp_orig"),$(USANSFolder+":DSM:I_msk")
1206        Duplicate/O $(USANSFolder+":DSM:S_exp_orig"),$(USANSFolder+":DSM:S_msk")
1207        Wave Q_msk=$(USANSFolder+":DSM:Q_msk")
1208        Wave I_msk=$(USANSFolder+":DSM:I_msk")
1209        Wave S_msk=$(USANSFolder+":DSM:S_msk")
1210       
1211        //finish up - trim the data sets and reassign the working set
1212        Wave MaskData=$(USANSFolder+":DSM:MaskData")
1213       
1214        RemoveMaskedPoints(MaskData,Q_msk,I_msk,S_msk)
1215
1216        //reset the number of points
1217        NVAR gNq = $(USANSFolder+":DSM:gNq")
1218        gNq = numpnts(Q_msk)
1219       
1220        Cursor/K A
1221        HideInfo
1222       
1223        return(0)
1224End
1225
1226
1227// not quite the same as revert
1228Function DSM_ClearMaskProc(ctrlName) : ButtonControl
1229        String ctrlName
1230       
1231        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1232       
1233        Wave MaskData=$(USANSFolder+":DSM:MaskData")
1234        MaskData = NaN
1235       
1236        return(0)
1237end
1238
1239// when the mask button is pressed, A must be on the graph
1240// Displays MaskData wave on the graph
1241//
1242Function DSM_MyMaskProc(ctrlName) : ButtonControl
1243        String ctrlName
1244       
1245        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1246       
1247        Wave data=$(USANSFolder+":DSM:I_exp_orig")
1248       
1249//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1250        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1251       
1252// need to get rid of old smoothed data if data is re-masked
1253        Execute "RemoveSmoothed()"
1254        SetDataFolder $(USANSFolder+":DSM")
1255        Killwaves/Z I_smth,Q_smth,S_smth
1256               
1257        if(aExists)             //mask the selected point
1258                // toggle NaN (keep) or Data value (= masked)
1259                Wave MaskData
1260                MaskData[pcsr(A)] = (numType(MaskData[pcsr(A)])==0) ? NaN : data[pcsr(A)]               //if NaN, doesn't plot
1261        else
1262                Cursor /A=1/H=1/L=1/P/W=Desmear_Graph A I_exp_orig leftx(I_exp_orig)
1263                ShowInfo
1264                //if the mask wave does not exist, make one
1265                if(exists("MaskData") != 1)
1266                        Duplicate/O Q_exp_orig MaskData
1267                        MaskData = NaN          //use all data
1268                endif
1269                Execute "AppendMask()" 
1270        endif
1271
1272        SetDataFolder root:
1273
1274        return(0)
1275End
1276
1277// when the mask button is pressed, A must be on the graph
1278// Displays MaskData wave on the graph
1279//
1280Function DSM_MaskLTCursor(ctrlName) : ButtonControl
1281        String ctrlName
1282
1283        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1284               
1285//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1286        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1287       
1288        if(!aExists)
1289                return(1)
1290        endif
1291// need to get rid of old smoothed data if data is re-masked
1292        Execute "RemoveSmoothed()"
1293        SetDataFolder $(USANSFolder+":DSM")
1294        Killwaves/Z I_smth,Q_smth,S_smth
1295
1296        Wave data=I_exp_orig
1297
1298        Variable pt,ii
1299        pt = pcsr(A)
1300        for(ii=pt;ii>=0;ii-=1)
1301                // toggle NaN (keep) or Data value (= masked)
1302                Wave MaskData
1303                MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii]              //if NaN, doesn't plot
1304        endfor
1305
1306        SetDataFolder root:
1307        return(0)
1308End
1309
1310// when the mask button is pressed, A must be on the graph
1311// Displays MaskData wave on the graph
1312//
1313Function DSM_MaskGTCursor(ctrlName) : ButtonControl
1314        String ctrlName
1315       
1316        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1317
1318       
1319//      Variable aExists= strlen(CsrInfo(A)) > 0                        //Igor 5
1320        Variable aExists= strlen(CsrWave(A)) > 0                        //Igor 4
1321       
1322        if(!aExists)
1323                return(1)
1324        endif
1325// need to get rid of old smoothed data if data is re-masked
1326        Execute "RemoveSmoothed()"
1327        SetDataFolder $(USANSFolder+":DSM")
1328        Killwaves/Z I_smth,Q_smth,S_smth
1329
1330        Wave data=I_exp_orig
1331
1332        Variable pt,ii,endPt
1333        endPt=numpnts(MaskData)
1334        pt = pcsr(A)
1335        for(ii=pt;ii<endPt;ii+=1)
1336                // toggle NaN (keep) or Data value (= masked)
1337                Wave MaskData
1338                MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii]              //if NaN, doesn't plot
1339        endfor
1340
1341        SetDataFolder root:
1342
1343        return(0)
1344End
1345
1346Function CleanUpJunk()
1347        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1348
1349        // clean up the old junk on the graph, /Z for no error
1350        Execute "RemoveExtrapolation()"
1351        Execute "RemoveDesmeared()"
1352        Execute "RemoveSmoothed()"
1353        Execute "RemoveMask()"
1354       
1355        //remove the cursor
1356        Cursor/K A
1357       
1358        //always initialize these
1359        String/G $(USANSFolder+":DSM:gStr1") = ""
1360        String/G $(USANSFolder+":DSM:gStr2") = ""
1361        String/G $(USANSFolder+":DSM:gIterStr") = ""
1362       
1363        // clean up the old waves from smoothing and desmearing steps
1364        SetDataFolder $(USANSFolder+":DSM")
1365        Killwaves/Z I_smth,I_dsm,I_dsm_sm,Q_smth,Q_dsm,S_smth,S_dsm,Yi_SS,Yq_SS
1366        Killwaves/Z Weights,FW,R_wave,S_wave,H_wave,Di,Ci
1367        Killwaves/Z Is_old,I_old,err
1368        Killwaves/Z Q_ext,I_ext,S_ext,hqq,hqi,lqq,lqi
1369        Killwaves/Z MaskData,Q_msk,I_msk,S_msk
1370        Killwaves/Z Q_work,I_work,S_work                                //working waves for desmearing step
1371        SetDataFolder root:
1372End
1373
1374// does not alter the data sets - reports a power law
1375// exponent and makes it global so it will automatically
1376// be used during the desmearing
1377//
1378// generates extr_hqi vs extr_hqq that are Appended to the graph
1379Function DSM_ExtrapolateButtonProc(ctrlName) : ButtonControl
1380        String ctrlName
1381
1382        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1383
1384        NVAR nend = $(USANSFolder+":DSM:gNptsExtrap")
1385        NVAR m_pred = $(USANSFolder+":DSM:gPowerM")
1386       
1387        SetDataFolder $(USANSFolder+":DSM")
1388//use masked data if it exists
1389        if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1390                wave Qw = $(USANSFolder+":DSM:Q_msk")
1391                Wave Iw = $(USANSFolder+":DSM:I_msk")
1392                Wave Sw = $(USANSFolder+":DSM:S_msk")
1393        else
1394                //start from the "_exp" waves
1395                if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1396                        wave Qw = $(USANSFolder+":DSM:Q_exp")
1397                        Wave Iw = $(USANSFolder+":DSM:I_exp")
1398                        Wave Sw = $(USANSFolder+":DSM:S_exp")
1399                endif
1400        endif
1401        SetDataFolder root:
1402       
1403        DSM_SetExtrWaves(Qw)
1404        m_pred = DSM_DoExtrapolate(Qw,Iw,Sw,nend)
1405        AppendExtrapolation()
1406       
1407        return(0)
1408End
1409
1410//smooths the data in steps as requested...
1411//
1412// typically do a simple Box smooth first,
1413// then do a smoothing spline, keeping the same number of data points
1414//
1415// chi-squared is reported - so you can see how "bad" the smoothing is
1416// smoothing of the data is done on a log(I) scale if requested
1417// setting doLog variable to 0 will return to linear smoothing
1418// (I see little difference)
1419//
1420// start from the "_smth" waves if they exist
1421// otheriwse start from the working waves
1422//
1423// create Q_smth, I_smth, S_smth
1424// keep track of the smoothing types and passes in the wave note
1425//
1426Function DSM_SmoothButtonProc(ctrlName) : ButtonControl
1427        String ctrlName
1428
1429        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1430
1431        SetDataFolder $(USANSFolder+":DSM")     
1432
1433        Variable ii,new_n,pass,nq_ext,offset,doLog=1
1434        String noteStr=""
1435       
1436        // want log scaling of intensity during smoothing?
1437        ControlInfo DSMControl_3g
1438        doLog = V_value
1439                       
1440////    if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1 && ! freshMask)
1441        if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1)
1442                //start from the smoothed waves --- just need the wave references
1443        else
1444                //start from the "msk", creating smth waves
1445                if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1446                        wave Q_msk,I_msk,S_msk
1447                        Duplicate/O I_msk,I_smth,Q_smth,S_smth
1448                        I_smth = I_msk
1449                        Q_smth = Q_msk
1450                        S_smth = S_msk
1451                else
1452                        //start from the "_exp" waves
1453                        if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1454                                wave Q_exp,I_exp,S_exp
1455                                Duplicate/O I_exp,I_smth,Q_smth,S_smth
1456                                I_smth = I_exp
1457                                Q_smth = Q_exp
1458                                S_smth = S_exp
1459                        endif
1460                endif
1461        endif
1462       
1463        wave Q_smth,I_smth,S_smth
1464
1465        // extend the data to avoid end effects
1466        //creates I_ext,Q_ext,S_ext with extended q-values
1467        // fit 15 pts at each end, typically add 10 pts to each end
1468        // does nothing if num_extr = 0
1469        ControlInfo/W=Desmear_Graph DSMControl_3d               //extendCheck
1470        if(V_value == 1)
1471                ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,10)
1472        else
1473                ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,0)            //don't extend, just use the original data
1474        endif
1475       
1476        //whether extending or not, the working data is "_ext", set by ExtendToSmooth()
1477        SetDataFolder $(USANSFolder+":DSM")     
1478        wave Q_ext,I_ext ,S_ext
1479       
1480        noteStr=note(I_smth)
1481        Note I_ext , noteStr
1482       
1483        WaveStats/Q I_ext
1484        if(V_min<=0)
1485                Print "negative itensity found, using linear scale for smoothing"
1486                doLog = 0
1487        endif
1488       
1489        if(doLog)
1490                //convert to log scale
1491                Duplicate/O I_ext I_log,I_log_err
1492                Wave I_log ,I_log_err
1493                I_log = log(I_ext)
1494                WaveStats/Q I_log
1495                offset = 2*(-V_min)             
1496                I_log += offset
1497                I_log_err = S_ext/(2.30*I_ext)
1498                I_ext = I_log
1499        endif   
1500        //
1501       
1502        ControlInfo/W=Desmear_Graph DSMControl_3b               //BoxCheck
1503        if(V_value == 1)
1504                //do a simple Box smooth first - number of points does not change
1505                // fills ends with neighboring value (E=3) (n=3 points in smoothing window)     
1506                Smooth/E=3/B 3, I_ext
1507               
1508                noteStr=note(I_ext)
1509                pass = NumberByKey("BOX", noteStr, "=", ";")
1510                noteStr = ReplaceNumberByKey("BOX", noteStr, pass+1, "=", ";")
1511                Note/K I_ext
1512                Note I_ext , noteStr
1513//              Print "box = ",noteStr
1514        endif
1515
1516        NVAR sParam = gSmoothFac                //already in the right DF
1517       
1518        ControlInfo/W=Desmear_Graph DSMControl_3c               //SSCheck
1519        if(V_value == 1)
1520//              Interpolate2/T=3/N=(nq)/I=2/F=1/SWAV=S_ext/Y=Yi_SS/X=Yq_SS Q_ext, I_ext         //Igor 5
1521//              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
1522        //Igor 4
1523                String str=""
1524                nq_ext = numpnts(Q_ext)
1525                str = "Interpolate/T=3/N=("+num2str(nq_ext)+")/I=1/F=("+num2str(sParam)+")/Y=Yi_SS/X=Yq_SS I_ext /X=Q_ext"             
1526                Execute str
1527//              Print Str
1528//              Interpolate2/T=3/N=(nq_ext)/I=1/F=(sParam)/Y=Yi_SS/X=Yq_SS Q_ext, I_ext
1529        // end Igor 4   
1530                wave yi_ss = yi_ss              // already in the right DF
1531                wave yq_ss = yq_ss
1532                // reassign the "working" waves with the result of interpolate, which has the same I/Q values
1533                I_ext = yi_ss
1534                Q_ext = yq_ss
1535               
1536                noteStr=note(I_ext)
1537                pass = NumberByKey("SPLINE", noteStr, "=", ";")
1538                noteStr = ReplaceNumberByKey("SPLINE", noteStr, pass+1, "=", ";")
1539                Note/K I_ext
1540                Note I_ext , noteStr
1541//              Print "spline = ",noteStr
1542        endif
1543       
1544        //undo the scaling
1545        If(doLog)
1546                I_ext -= offset
1547                I_ext = 10^I_ext
1548        endif
1549       
1550        // at this point, I_ext has too many points - and we need to just return the
1551        // center chunk that is the original q-values
1552        // as assign this to the working set for the desmear step
1553        // and to the _smth set in case another smoothng pass is desired
1554        // Q_smth has not changed, S_smth has not changed
1555        I_smth = interp(Q_smth[p], Q_ext, I_ext )
1556       
1557        Note/K I_smth
1558        Note I_smth , noteStr
1559//      Print "end of smoothed, note = ",note(I_smth)
1560       
1561        Variable nq = numpnts($(USANSFolder+":DSM:Q_smth"))
1562//      Print "nq after smoothing = ",nq
1563
1564        //reset the global
1565        NVAR gNq = gNq
1566        gNq = nq
1567        //report the chi^2 difference between the smoothed curve and the experimental data
1568        NVAR chi2 = gChi2Smooth
1569        chi2 = SmoothedChi2(I_smth)
1570       
1571        Execute "AppendSmoothed()"
1572       
1573        setDataFolder root:
1574        return(0)
1575End
1576
1577// undo the smoothing and start over - useful if you've smoothed
1578// too aggressively and need to back off
1579Function DSM_SmoothUndoButtonProc(ctrlName) : ButtonControl
1580        String ctrlName
1581
1582        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1583       
1584        Execute "RemoveSmoothed()"
1585        SetDataFolder $(USANSFolder+":DSM")
1586        Killwaves/Z I_smth,Q_smth,S_smth,Q_ext,I_ext,S_ext,Yi_SS,Yq_SS
1587        SetDataFolder root:
1588        return (0)
1589end
1590
1591//calculate the chi^2 value for the smoothed data
1592Function SmoothedChi2(I_smth)
1593        Wave I_smth
1594       
1595        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1596       
1597        //start from the "msk", if they exist
1598        if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1599                Wave iw = $(USANSFolder+":DSM:I_msk")
1600                Wave sw = $(USANSFolder+":DSM:S_msk")
1601        else
1602                //start from the "_exp" waves
1603                if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1604                        Wave iw = $(USANSFolder+":DSM:I_exp")
1605                        Wave sw = $(USANSFolder+":DSM:S_exp")
1606                endif
1607        endif
1608       
1609        Variable ii,chi2,num
1610        chi2=0
1611        num = numpnts(iw)
1612        for(ii=0;ii<num;ii+=1)
1613                chi2 += (iw[ii] - I_smth[ii])^2 / (sw[ii])^2
1614        endfor
1615        Chi2 /= (num-1)
1616        return (chi2)
1617end
1618
1619// I_dsm is the desmeared data
1620//
1621// step (7) - desmearing is done, write out the result
1622//
1623Function DSM_SaveButtonProc(ctrlName) : ButtonControl
1624        String ctrlName
1625
1626        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1627
1628        String saveStr
1629        SVAR curFile = $(USANSFolder+":DSM:gCurFile")
1630        saveStr = CleanupName((curFile),0)              //the output filename
1631        //
1632        WriteUSANSDesmeared(saveStr,0,0,1)                      //use the full set (lo=hi=0) and present a dialog
1633       
1634        SetDataFolder root:
1635        return(0)
1636End
1637
1638Function DSM_HelpButtonProc(ctrlName) : ButtonControl
1639        String ctrlName
1640
1641        DisplayHelpTopic/Z/K=1 "Desmearing USANS Data"
1642        if(V_flag !=0)
1643                DoAlert 0,"The Desmearing USANS Data Help file could not be found"
1644        endif
1645        return(0)
1646End
1647
1648
1649//toggles the log/lin display of the loaded data set
1650Function DSM_LoadCheckProc(ctrlName,checked) : CheckBoxControl
1651        String ctrlName
1652        Variable checked
1653
1654        ModifyGraph log=(checked)
1655        return(0)
1656End
1657
1658
1659Function DSM_ExtrapolationCheckProc(ctrlName,checked) : CheckBoxControl
1660        String ctrlName
1661        Variable checked
1662
1663        if(checked)
1664                AppendExtrapolation()
1665        else
1666                RemoveExtrapolation()
1667        endif
1668        return(0)
1669End
1670
1671
1672// takes as input the "working waves"
1673//
1674// creates intermediate waves to work on
1675//
1676// output of Q_dsm,I_dsm,S_dsm (and I_dsm_sm, the result of smearing I_dsm)
1677//
1678Function DSM_DesmearButtonProc(ctrlName) : ButtonControl
1679        String ctrlName
1680       
1681        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1682       
1683        SetDataFolder $(USANSFolder+":DSM")
1684        if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1)
1685                wave Q_smth , I_smth ,S_smth
1686                Duplicate/O I_smth,I_work,Q_work,S_work
1687                I_work = I_smth
1688                Q_work = Q_smth
1689                S_work = S_smth
1690        else
1691                //start from the "msk", creating work waves
1692                if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1693                        wave Q_msk,I_msk,S_msk
1694                        Duplicate/O I_msk,I_work,Q_work,S_work
1695                        I_work = I_msk
1696                        Q_work = Q_msk
1697                        S_work = S_msk
1698                else
1699                        //start from the "_exp" waves
1700                        if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1701                                wave Q_exp,I_exp,S_exp
1702                                Duplicate/O I_exp,I_work,Q_work,S_work
1703                                I_work = I_exp
1704                                Q_work = Q_exp
1705                                S_work = S_exp
1706                        endif
1707                endif
1708        endif
1709        //SetDataFolder root:
1710       
1711        NVAR nq = gNq
1712        NVAR m = gPowerM
1713        NVAR chi2_target = gChi2Target
1714        NVAR maxFastIter = gMaxFastIter
1715        NVAR maxSlowIter = gMaxSlowIter
1716       
1717        //      SET WEIGHTING OF EXPERIMENTAL DATA.
1718        Duplicate/O Q_work weights
1719        Wave weights = weights
1720        weights = 1/S_work^2
1721
1722//      calculate weighting array for smearing of data
1723        Make/O/D/N=(nq,nq) FW
1724        Wave FW
1725        Weights_L(m,FW,Q_work)
1726
1727//      ^^^^   Iterative desmearing   ^^^^*
1728        Variable chi2_old,chi2_new,done,iter
1729//      FOR 0TH ITERATION, EXPERIMENTAL DATA IS USED FOR Y_OLD, create ys_old for result
1730//      y_old = I_old, y_new = I_dsm, I_dsm_sm = ys_new,
1731// duplicate preserves the wave note!
1732        Duplicate/O I_work I_old,Is_old,I_dsm,I_dsm_sm
1733        Duplicate/O Q_work Q_dsm,S_dsm          //sets Q_dsm correctly
1734        wave S_dsm,I_old,Is_old,I_dsm,I_dsm_sm
1735        I_old = I_work
1736        Is_old = 0
1737        I_dsm = 0
1738        I_dsm_sm = 0
1739       
1740//      Smear 0TH iter guess Y --> YS
1741        DSM_Smear(I_old,Is_old,NQ,FW)
1742        chi2_OLD = chi2(Is_old,I_work,weights,NQ)
1743       
1744        printf "starting chi^2 = %g\r",chi2_old
1745        Print "Starting fast convergence..... "
1746
1747        done = 0                //false
1748        iter = 0
1749        Do                                      // while not done, use Fast convergence
1750                I_dsm = I_work * I_old / Is_old //        Calculate corrected guess (no need for do-loop)
1751
1752//        Smear iter guess I_dsm --> I_dsm_sm and see how well I_dsm_sm matches experimental data
1753                DSM_Smear(I_dsm,I_dsm_sm,NQ,FW)
1754                chi2_new = chi2(I_dsm_sm,I_work,weights,NQ)
1755
1756//        Stop iteration if fit from new iteration has worse fit...
1757                If (chi2_new > chi2_old)
1758                        Done = 1
1759                Endif
1760               
1761//        Stop iteration if fit is better than target value...
1762                If (chi2_new < chi2_target)
1763                        Done = 1
1764                Else
1765                        Iter += 1
1766                        Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1767                        I_old = I_dsm
1768                        Is_old = I_dsm_sm
1769                        if(iter>maxFastIter)
1770                                break
1771                        endif
1772                        CHI2_OLD = CHI2_NEW
1773                EndIf
1774        while( !done )
1775       
1776        // append I_dsm,I_exp to the graph
1777        Execute "AppendDesmeared()"
1778        DoUpdate
1779       
1780        SetDataFolder $(USANSFolder+":DSM")
1781
1782       
1783// step (6) - refine the desmearing using slow convergence
1784        Print "Starting slow convergence..... "
1785        done = 0                //  ! reset flag for slow convergence
1786        Do
1787                I_dsm = I_old + (I_work - Is_old)       //        Calculate corrected guess
1788
1789//        Smear iter guess Y --> YS
1790                DSM_Smear(I_dsm,I_dsm_sm,NQ,FW)
1791                chi2_new = chi2(I_dsm_sm,I_work,weights,NQ)
1792
1793//        Stop iteration if fit from new iteration has worse fit...
1794                If (chi2_new > chi2_old)
1795                        Done = 1
1796                Endif
1797               
1798//        Stop iteration if fit is better than target value...
1799                If (chi2_new < chi2_target)
1800                        Done = 1
1801                Else
1802                        Iter += 1
1803                       
1804                        if(mod(iter, 50 ) ==0 )
1805                                Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1806                                DoUpdate
1807                        endif
1808                        I_old = I_dsm
1809                        Is_old = I_dsm_sm
1810
1811                        CHI2_OLD = CHI2_NEW
1812                EndIf
1813                if(iter>maxSlowIter)
1814                        break
1815                endif
1816        While ( !done )
1817
1818        Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1819
1820        // adjust the error
1821        SetDataFolder $(USANSFolder+":DSM")
1822        Duplicate/O S_work err
1823        S_dsm = abs(err/I_work*I_dsm)           //proportional error
1824//John simply keeps the same error as was read in from the smeared data set - is this right?
1825//      S_dsm = S_Work
1826       
1827        NVAR gChi =  gChi2Final         //chi^2 final
1828        NVAR gIter = gIterations                //total number of iterations
1829        gChi = chi2_new
1830        gIter = Iter
1831       
1832        SVAR str = gIterStr
1833        sprintf str,"%d iterations required",iter
1834       
1835        SetDataFolder root:
1836        return(0)
1837End
1838
1839Function DSM_RevertButtonProc(ctrlName) : ButtonControl
1840        String ctrlName
1841       
1842        SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder       
1843       
1844       
1845        CleanUpJunk()
1846
1847        SetDataFolder $(USANSFolder+":DSM")
1848
1849        //reset the working waves to the original
1850        wave Q_exp_orig,I_exp_orig,S_exp_orig
1851
1852        Duplicate/O Q_exp_orig Q_exp
1853        Duplicate/O I_exp_orig I_exp
1854        Duplicate/O S_exp_orig S_exp
1855        // kill the data folder
1856        // re-initialize?
1857
1858        SetDataFolder root:
1859       
1860        return(0)
1861End
Note: See TracBrowser for help on using the repository browser.