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

Last change on this file was 1188, checked in by srkline, 4 years ago

increasing setVariable control size on the masking panel so that the sector parameters are more visible on windows computers.

Added the capability to save smoothed (not desmeared) data files from the USANs desmearing panel. These smoothed data files are re-saved with the -dQv value that they were loaded with, since they are still smeared USANS data sets. They are saved with the ".smth" extension.

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