source: sans/Release/trunk/NCNR_User_Procedures/USANS Procedures v2.21/LakeDesmearing_JB.ipf @ 284

Last change on this file since 284 was 265, checked in by srkline, 15 years ago

Major structural changes, hopefully correct, to prepare for v2.21 release (not tagged yet)

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