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

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

Two changes:
(1) Changed the behavior of DIV file creation. Now a panel is presented that asks for all of the information - run numbers, transmission, XY box, etc. and then does all of the reduction/replacing/saving with one click. Cedric's ticket #195

(2) Some fixes in USANS notebook generation - filtering out junk.

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