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

Last change on this file since 1103 was 742, checked in by srkline, 12 years ago

fix for ticket #279: calculation of the slope when desmearing.

Now the behavior is for the desmearing panel to use the slope as found when loading the .cor data. For a fresh file, it proceeds as usual. For a previously-loaded file, be sure to NOT re-load the file (<no> from the dialog) so that the slope is not recalculated, and the data is simply plotted. Then the slope information is copied along with the data to the desmearing folder.

The changed behavior is when selecting the "Extrapolate" tab. Now the extrapolation is automatically calculated from the loading parameters (not re-fit) and displayed on the graph. If the extrapolation is good, nothing needs to be done. In fact, since the extrapolation (slope) should be acceptable from the load, the mask, extrapolate, and smooth tabs can all be skipped, proceeding directly to the desmear tab.

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