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

Last change on this file since 611 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

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