source: sans/USANSReduction/trunk/Put in User Procedures Folder/USANS Procedures v2.20/LakeDesmearing_JB.ipf @ 53

Last change on this file since 53 was 53, checked in by srkline, 16 years ago

changed the desmeared error calculation to be abs(err) so that
none of the errors were negative

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