source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/WorkFileUtils.ipf @ 890

Last change on this file since 890 was 890, checked in by srkline, 10 years ago

Added a procedure file to allow "scripting" of the MonteCarlo? simulation: MC_SimulationScripting.ipf. See the beginning of the file for instructions, and then the examples for numbered instructions of what needs to be edited to write your own script. Not casual-user-friendly, but makes it easire for more advanced users to simulate an entire experiment.

File size: 50.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//*******************
6// Vers 1.2 100901
7//
8//*******************
9// Utility procedures for handling of workfiles (each is housed in a separate datafolder)
10//
11// - adding data to a workfile
12// - copying workfiles to another folder
13//
14// - absolute scaling
15// - DIV detector sensitivity corrections
16//
17// - the WorkFile Math panel for simple image math
18// -
19// - adding work.drk data without normalizing to monitor counts
20//***************************
21
22//
23//
24// MARCH 2011 - changed the references that manipulate the data to explcitly work on the linear_data wave
25//                                      at the end of routines that maniplulate linear_data, it is copied back to "data" which is displayed
26//
27//
28
29//testing procedure, not called anymore
30Proc Add_to_Workfile(type, add)
31        String type,add
32        Prompt type,"WORK data type",popup,"SAM;EMP;BGD"
33        Prompt add,"Add to current WORK contents?",popup,"No;Yes"
34       
35        //macro will take whatever is in RAW folder and "ADD" it to the folder specified
36        //in the popup menu
37       
38        //"add" = yes/no, don't add to previous runs
39        //switch here - two separate functions to avoid (my) confusion
40        Variable err
41        if(cmpstr(add,"No")==0)
42                //don't add to prev work contents, copy RAW contents to work and convert
43                err = Raw_to_work(type)
44        else
45                //yes, add RAW to the current work folder contents
46                err = Add_raw_to_work(type)
47        endif
48       
49        String newTitle = "WORK_"+type
50        DoWindow/F SANS_Data
51        DoWindow/T SANS_Data, newTitle
52        KillStrings/Z newTitle
53       
54        //need to update the display with "data" from the correct dataFolder
55        fRawWindowHook()
56       
57End
58
59//will "ADD" the current contents of the RAW folder to the newType work folder
60//and will ADD the RAW contents to the existing content of the newType folder
61// - used when adding multiple runs together
62//(the function Raw_to_work(type) makes a fresh workfile)
63//
64//the current display type is updated to newType (global)
65Function Add_raw_to_work(newType)
66        String newType
67               
68        String destPath=""
69       
70        // if the desired workfile doesn't exist, let the user know, and just make a new one
71        if(WaveExists($("root:Packages:NIST:"+newType + ":data")) == 0)
72                Print "There is no old work file to add to - a new one will be created"
73                //call Raw_to_work(), then return from this function
74                Raw_to_Work(newType)
75                Return(0)               //does not generate an error - a single file was converted to work.newtype
76        Endif
77       
78        NVAR pixelsX = root:myGlobals:gNPixelsX
79        NVAR pixelsY = root:myGlobals:gNPixelsY
80       
81        //now make references to data in newType folder
82        DestPath="root:Packages:NIST:"+newType 
83        WAVE data=$(destPath +":linear_data")                   // these wave references point to the EXISTING work data
84        WAVE data_copy=$(destPath +":data")                     // these wave references point to the EXISTING work data
85        WAVE dest_data_err=$(destPath +":linear_data_error")                    // these wave references point to the EXISTING work data
86        WAVE/T textread=$(destPath + ":textread")
87        WAVE integersread=$(destPath + ":integersread")
88        WAVE realsread=$(destPath + ":realsread")
89       
90        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
91        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy,xshift,yshift
92
93// 08/01 detector constants are now returned from a function, based on the detector type and beamline
94//      dt_ornl = 3.4e-6                //deadtime of Ordella detectors as of 30-AUG-99
95//      dt_ill=3.0e-6                   //Cerca detector deadtime constant as of 30-AUG-99
96
97        defmon=1e8                      //default monitor counts
98       
99        //Yes, add to previous run(s) in work, that does exist
100        //use the actual monitor count run.savmon rather than the normalized monitor count
101        //in run.moncnt and unscale the work data
102       
103        total_mon = realsread[1]        //saved monitor count
104        uscale = total_mon/defmon               //unscaling factor
105        total_det = uscale*realsread[2]         //unscaled detector count
106        total_trn = uscale*realsread[39]        //unscaled trans det count
107        total_numruns = integersread[3] //number of runs in workfile
108        total_rtime = integersread[2]           //total counting time in workfile
109        //retrieve workfile beamcenter
110        wrk_beamx = realsread[16]
111        wrk_beamy = realsread[17]
112        //unscale the workfile data in "newType"
113        //
114        //check for log-scaling and adjust if necessary
115        ConvertFolderToLinearScale(newType)
116        //
117        //then unscale the data array
118        data *= uscale
119        dest_data_err *= uscale
120       
121        //DetCorr() has not been applied to the data in RAW , do it now in a local reference to the raw data
122        WAVE raw_data = $"root:Packages:NIST:RAW:linear_data"
123        WAVE raw_data_err = $"root:Packages:NIST:RAW:linear_data_error"
124        WAVE raw_reals =  $"root:Packages:NIST:RAW:realsread"
125        WAVE/T raw_text = $"root:Packages:NIST:RAW:textread"
126        WAVE raw_ints = $"root:Packages:NIST:RAW:integersread"
127       
128        //check for log-scaling of the raw data - make sure it's linear
129        ConvertFolderToLinearScale("RAW")
130       
131        // switches to control what is done, don't do the transmission correction for the BGD measurement
132        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCorr
133        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCorr
134        Variable doTrans = gDoTrans
135        if(cmpstr("BGD",newtype) == 0)
136                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
137        endif   
138       
139        DetCorr(raw_data,raw_data_err,raw_reals,doEfficiency,doTrans)   //applies correction to raw_data, and overwrites it
140       
141        //if RAW data is ILL type detector, correct raw_data for same counts being written to 4 pixels
142        if(cmpstr(raw_text[9], "ILL   ") == 0 )         //text field in header is 6 characters "ILL---"
143                raw_data /= 4
144                raw_data_err /= 4
145        endif
146       
147        //deadtime corrections to raw data
148        deadTime = DetectorDeadtime(raw_text[3],raw_text[9],dateAndTimeStr=raw_text[1])         //pick the correct detector deadtime, switch on date too
149        itim = raw_ints[2]
150        cntrate = sum(raw_data,-inf,inf)/itim           //080802 use data sum, rather than scaler value
151        dscale = 1/(1-deadTime*cntrate)
152
153#if (exists("ILL_D22")==6)
154        Variable tubeSum
155        // for D22 detector might need to use cntrate/128 as it is the tube response
156        for(ii=0;ii<pixelsX;ii+=1)
157                //sum the counts in each tube
158                tubeSum = 0
159                for(jj=0;jj<pixelsY;jj+=1)
160                        tubeSum += data[jj][ii]
161                endfor
162                // countrate in tube ii
163                cntrate = tubeSum/itim
164                // deadtime scaling in tube ii
165                dscale = 1/(1-deadTime*cntrate)
166                // multiply data[ii][] by the dead time
167                raw_data[][ii] *= dscale
168                raw_data_err[][ii] *= dscale
169        endfor
170#else
171        // dead time correction on all other RAW data, including NCNR
172        raw_data *= dscale
173        raw_data_err *= dscale
174#endif
175
176        //update totals by adding RAW values to the local ones (write to work header at end of function)
177        total_mon += raw_reals[0]
178#if (exists("ILL_D22")==6)
179        total_det += sum(raw_data,-inf,inf)                     //add the newly scaled detector array
180#else
181        total_det += dscale*raw_reals[2]
182#endif
183        total_trn += raw_reals[39]
184        total_rtime += raw_ints[2]
185        total_numruns +=1
186       
187        //do the beamcenter shifting if there is a mismatch
188        //and then add the two data sets together, changing "data" since it is the workfile data
189        xshift = raw_reals[16] - wrk_beamx
190        yshift = raw_reals[17] - wrk_beamy
191       
192        If((xshift != 0) || (yshift != 0))
193                DoAlert 1,"Do you want to ignore the beam center mismatch?"
194                if(V_flag==1)
195                        xshift=0
196                        yshift=0
197                endif
198        endif
199       
200        If((xshift == 0) && (yshift == 0))              //no shift, just add them
201                data += raw_data                //deadtime correction has already been done to the raw data
202                dest_data_err = sqrt(dest_data_err^2 + raw_data_err^2)                  // error of the sum
203        else
204                //shift the beamcenter, then add
205                Make/O/N=1 $(destPath + ":noadd")               //needed to get noadd condition back from ShiftSum()
206                WAVE noadd = $(destPath + ":noadd")
207                Variable sh_sum                 //returned value
208                Print "BEAM CENTER MISMATCH - - BEAM CENTER WILL BE SHIFTED TO THIS FILE'S VALUE"
209                //ii,jj are just indices here, not physical locations - so [0,127] is fine
210                ii=0
211                do
212                        jj=0
213                        do
214                                //get the contribution of shifted data
215                                sh_sum = ShiftSum(data,ii,jj,xshift,yshift,noadd)
216                                if(noadd[0])
217                                        //don't do anything to data[][]
218                                else
219                                        //add the raw_data + shifted sum (and do the deadtime correction on both)
220                                        data[ii][jj] += (raw_data[ii][jj]+sh_sum)               //do the deadtime correction on RAW here
221                                        dest_data_err[ii][jj] = sqrt(dest_data_err[ii][jj]^2 + raw_data_err[ii][jj]^2)                  // error of the sum
222                                Endif
223                                jj+=1
224                        while(jj<pixelsY)
225                        ii+=1
226                while(ii<pixelsX)
227        Endif
228       
229        //scale the data to the default montor counts
230        scale = defmon/total_mon
231        data *= scale
232        dest_data_err *= scale
233       
234        // keep "data" and linear_data in sync in the destination folder
235        data_copy = data
236       
237        //all is done, except for the bookkeeping of updating the header info in the work folder
238        textread[1] = date() + " " + time()             //date + time stamp
239        integersread[3] = total_numruns                                         //numruns = more than one
240        realsread[1] = total_mon                        //save the true monitor count
241        realsread[0] = defmon                                   //monitor ct = defmon
242        integersread[2] = total_rtime                   // total counting time
243        realsread[2] = scale*total_det                  //scaled detector counts
244        realsread[39] = scale*total_trn                 //scaled transmission counts
245       
246        //Add the added raw filename to the list of files in the workfile
247        String newfile = ";" + raw_text[0]
248        SVAR oldList = $(destPath + ":fileList")
249        String/G $(destPath + ":fileList") = oldList + newfile
250       
251        //reset the current displaytype to "newtype"
252        String/G root:myGlobals:gDataDisplayType=newType
253       
254        //return to root folder (redundant)
255        SetDataFolder root:
256       
257        Return(0)
258End
259
260//will copy the current contents of the RAW folder to the newType work folder
261//and do the geometric corrections and normalization to monitor counts
262//(the function Add_Raw_to_work(type) adds multiple runs together)
263//
264//the current display type is updated to newType (global)
265//
266Function Raw_to_work(newType)
267        String newType
268       
269        Variable deadTime,defmon,total_mon,total_det,total_trn,total_numruns,total_rtime
270        Variable ii,jj,itim,cntrate,dscale,scale,uscale,wrk_beamx,wrk_beamy
271        String destPath
272       
273// 08/01 detector constants are now returned from a function, based on the detector type and beamline
274//      dt_ornl = 3.4e-6                //deadtime of Ordella detectors as of 30-AUG-99
275//      dt_ill=3.0e-6                   //Cerca detector deadtime constant as of 30-AUG-99
276        defmon=1e8                      //default monitor counts
277       
278        //initialize values before normalization
279        total_mon=0
280        total_det=0
281        total_trn=0
282        total_numruns=0
283        total_rtime=0
284       
285        //Not adding multiple runs, so wipe out the old contents of the work folder and
286        // replace with the contents of raw
287
288        destPath = "root:Packages:NIST:" + newType
289       
290        //check for log-scaling of the RAW data and adjust if necessary
291        ConvertFolderToLinearScale("RAW")
292        //then continue
293       
294        NVAR pixelsX = root:myGlobals:gNPixelsX
295        NVAR pixelsY = root:myGlobals:gNPixelsY
296       
297        //copy from current dir (RAW) to work, defined by destpath
298        DestPath = "root:Packages:NIST:"+newType
299        Duplicate/O $"root:Packages:NIST:RAW:data",$(destPath + ":data")
300        Duplicate/O $"root:Packages:NIST:RAW:linear_data_error",$(destPath + ":linear_data_error")
301        Duplicate/O $"root:Packages:NIST:RAW:linear_data",$(destPath + ":linear_data")
302//      Duplicate/O $"root:Packages:NIST:RAW:vlegend",$(destPath + ":vlegend")
303        Duplicate/O $"root:Packages:NIST:RAW:textread",$(destPath + ":textread")
304        Duplicate/O $"root:Packages:NIST:RAW:integersread",$(destPath + ":integersread")
305        Duplicate/O $"root:Packages:NIST:RAW:realsread",$(destPath + ":realsread")
306        Variable/G $(destPath + ":gIsLogscale")=0                       //overwite flag in newType folder, data converted (above) to linear scale
307
308        WAVE data=$(destPath + ":linear_data")                                                  // these wave references point to the data in work
309        WAVE data_copy=$(destPath + ":data")                                                    // these wave references point to the data in work
310        WAVE data_err=$(destPath + ":linear_data_error")
311        WAVE/T textread=$(destPath + ":textread")                               //that are to be directly operated on
312        WAVE integersread=$(destPath + ":integersread")
313        WAVE realsread=$(destPath + ":realsread")
314        String/G $(destPath + ":fileList") = textread[0]                        //a list of names of the files in the work file (1)
315       
316        //apply nonlinear, Jacobian corrections ---
317        // switches to control what is done, don't do the transmission correction for the BGD measurement
318        NVAR doEfficiency = root:Packages:NIST:gDoDetectorEffCorr
319        NVAR gDoTrans = root:Packages:NIST:gDoTransmissionCorr
320        Variable doTrans = gDoTrans
321        if(cmpstr("BGD",newtype) == 0)
322                doTrans = 0             //skip the trans correction for the BGD file but don't change the value of the global
323        endif
324       
325        DetCorr(data,data_err,realsread,doEfficiency,doTrans)           //the parameters are waves, and will be changed by the function
326
327        //if ILL type detector, correct for same counts being written to 4 pixels
328        if(cmpstr(textread[9], "ILL   ") == 0 )         //text field in header is 6 characters "ILL---"
329                data /= 4
330                data_err /= 4           //rescale error
331        endif
332       
333        //deadtime corrections
334        itim = integersread[2]
335        cntrate = sum(data,-inf,inf)/itim               //use sum of detector counts rather than scaler value
336        deadtime = DetectorDeadtime(textread[3],textread[9],dateAndTimeStr=textRead[1]) //pick the correct deadtime
337        dscale = 1/(1-deadTime*cntrate)
338       
339#if (exists("ILL_D22")==6)
340        Variable tubeSum
341        // for D22 detector might need to use cntrate/128 as it is the tube response
342        for(ii=0;ii<pixelsX;ii+=1)
343                //sum the counts in each tube
344                tubeSum = 0
345                for(jj=0;jj<pixelsY;jj+=1)
346                        tubeSum += data[jj][ii]
347                endfor
348                // countrate in tube ii
349                cntrate = tubeSum/itim
350                // deadtime scaling in tube ii
351                dscale = 1/(1-deadTime*cntrate)
352                // multiply data[ii][] by the dead time
353                data[][ii] *= dscale
354                data_err[][ii] *= dscale
355        endfor
356#endif
357
358        // NO xcenter,ycenter shifting is done - this is the first (and only) file in the work folder
359       
360        //only ONE data file- no addition of multiple runs in this function, so data is
361        //just simply corrected for deadtime.
362#if (exists("ILL_D22")==0)              //ILL correction done tube-by-tube above
363        data *= dscale          //deadtime correction for everyone else, including NCNR
364        data_err *= dscale
365#endif
366       
367        //update totals to put in the work header (at the end of the function)
368        total_mon += realsread[0]
369#if (exists("ILL_D22")==6)
370        total_det += sum(data,-inf,inf)                 //add the newly scaled detector array
371#else
372        total_det += dscale*realsread[2]
373#endif
374        total_trn += realsread[39]
375        total_rtime += integersread[2]
376        total_numruns +=1
377       
378       
379        //scale the data to the default montor counts
380        scale = defmon/total_mon
381        data *= scale
382        data_err *= scale               //assumes total monitor count is so large there is essentially no error
383       
384        // to keep "data" and linear_data in sync
385        data_copy = data
386       
387        //all is done, except for the bookkeeping, updating the header information in the work folder
388        textread[1] = date() + " " + time()             //date + time stamp
389        integersread[3] = total_numruns                                         //numruns = 1
390        realsread[1] = total_mon                        //save the true monitor count
391        realsread[0] = defmon                                   //monitor ct = defmon
392        integersread[2] = total_rtime                   // total counting time
393        realsread[2] = scale*total_det                  //scaled detector counts
394        realsread[39] = scale*total_trn                 //scaled transmission counts
395       
396        //reset the current displaytype to "newtype"
397        String/G root:myGlobals:gDataDisplayType=newType
398       
399        //return to root folder (redundant)
400        SetDataFolder root:
401       
402        Return(0)
403End
404
405//used for adding DRK (beam shutter CLOSED) data to a workfile
406//force the monitor count to 1, since it's irrelevant
407// run data through normal "add" step, then unscale default monitor counts
408//to get the data back on a simple time basis
409//
410Function Raw_to_Work_NoNorm(type)
411        String type
412       
413        WAVE reals=$("root:Packages:NIST:RAW:realsread")
414        reals[1]=1              //true monitor counts, still in raw
415        Raw_to_work(type)
416        //data is now in "type" folder
417        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
418        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
419        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
420        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
421       
422        Variable norm_mon,tot_mon,scale
423       
424        norm_mon = new_reals[0]         //should be 1e8
425        tot_mon = new_reals[1]          //should be 1
426        scale= norm_mon/tot_mon
427       
428        data /= scale           //unscale the data
429        data_err /= scale
430       
431        // to keep "data" and linear_data in sync
432        data_copy = data
433       
434        return(0)
435End
436
437//used for adding DRK (beam shutter CLOSED) data to a workfile
438//force the monitor count to 1, since it's irrelevant
439// run data through normal "add" step, then unscale default monitor counts
440//to get the data back on a simple time basis
441//
442Function Add_Raw_to_Work_NoNorm(type)
443        String type
444       
445        WAVE reals=$("root:Packages:NIST:RAW:realsread")
446        reals[1]=1              //true monitor counts, still in raw
447        Add_Raw_to_work(type)
448        //data is now in "type" folder
449        WAVE data=$("root:Packages:NIST:"+type+":linear_data")
450        WAVE data_copy=$("root:Packages:NIST:"+type+":data")
451        WAVE data_err=$("root:Packages:NIST:"+type+":linear_data_error")
452        WAVE new_reals=$("root:Packages:NIST:"+type+":realsread")
453       
454        Variable norm_mon,tot_mon,scale
455       
456        norm_mon = new_reals[0]         //should be 1e8
457        tot_mon = new_reals[1]          //should be equal to the number of runs (1 count per run)
458        scale= norm_mon/tot_mon
459       
460        data /= scale           //unscale the data
461        data_err /= scale
462       
463        // to keep "data" and linear_data in sync
464        data_copy = data
465       
466        return(0)
467End
468
469//performs solid angle and non-linear detector corrections to raw data as it is "added" to a work folder
470//function is called by Raw_to_work() and Add_raw_to_work() functions
471//works on the actual data array, assumes that is is already on LINEAR scale
472//
473Function DetCorr(data,data_err,realsread,doEfficiency,doTrans)
474        Wave data,data_err,realsread
475        Variable doEfficiency,doTrans
476       
477        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,xx0,yy0
478        Variable ii,jj,dtdist,dtdis2
479        Variable xi,xd,yd,rad,ratio,domega,xy
480        Variable lambda,trans,trans_err,lat_err,tmp_err,lat_corr
481       
482//      Print "...doing jacobian and non-linear corrections"
483
484        NVAR pixelsX = root:myGlobals:gNPixelsX
485        NVAR pixelsY = root:myGlobals:gNPixelsY
486       
487        //set up values to send to auxiliary trig functions
488        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
489        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
490
491        x0 = realsread[16]
492        y0 = realsread[17]
493        sx = realsread[10]
494        sx3 = realsread[11]
495        sy = realsread[13]
496        sy3 = realsread[14]
497       
498        dtdist = 1000*realsread[18]     //sdd in mm
499        dtdis2 = dtdist^2
500       
501        lambda = realsRead[26]
502        trans = RealsRead[4]
503        trans_err = RealsRead[41]               //new, March 2011
504       
505        xx0 = dc_fx(x0,sx,sx3,xcenter)
506        yy0 = dc_fy(y0,sy,sy3,ycenter)
507       
508
509        //waves to contain repeated function calls
510        Make/O/N=(pixelsX) fyy,xx,yy            //Assumes square detector !!!
511        ii=0
512        do
513                xi = ii
514                fyy[ii] = dc_fy(ii+1,sy,sy3,ycenter)
515                xx[ii] = dc_fxn(ii+1,sx,sx3,xcenter)
516                yy[ii] = dc_fym(ii+1,sy,sy3,ycenter)
517                ii+=1
518        while(ii<pixelsX)
519       
520        Make/O/N=(pixelsX,pixelsY) SolidAngle           // testing only
521       
522        ii=0
523        do
524                xi = ii
525                xd = dc_fx(ii+1,sx,sx3,xcenter)-xx0
526                jj=0
527                do
528                        yd = fyy[jj]-yy0
529                        //rad is the distance of pixel ij from the sample
530                        //domega is the ratio of the solid angle of pixel ij versus center pixel
531                        // product xy = 1 for a detector with a linear spatial response (modern Ordela)
532                        // solid angle calculated, dW^3 >=1, so multiply data to raise measured values to correct values.
533                        rad = sqrt(dtdis2 + xd^2 + yd^2)
534                        domega = rad/dtdist
535                        ratio = domega^3
536                        xy = xx[ii]*yy[jj]
537                       
538                        data[ii][jj] *= xy*ratio
539                       
540                        solidAngle[ii][jj] = xy*ratio           //testing only 
541                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
542                       
543                       
544                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
545                        // correction inserted 11/2007 SRK
546                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
547                        // so divide here to get the correct answer (5/22/08 SRK)
548                        if(doEfficiency)
549#if (exists("ILL_D22")==6)
550                                data[ii][jj]  /= DetEffCorrILL(lambda,dtdist,xd)                //tube-by-tube corrections
551                                data_err[ii][jj] /= DetEffCorrILL(lambda,dtdist,xd)                     //assumes correction is exact
552//                solidAngle[ii][jj] = DetEffCorrILL(lambda,dtdist,xd)
553#else
554                                data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
555                                data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
556//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
557#endif
558                        endif
559                       
560                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
561                        // so divide here to get the correct answer
562                        if(doTrans)
563                       
564                                if(trans<0.1 && ii==0 && jj==0)
565                                        Print "***transmission is less than 0.1*** and is a significant correction"
566                                endif
567                               
568                                if(trans==0)
569                                        if(ii==0 && jj==0)
570                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
571                                        endif
572                                        trans = 1
573                                endif
574                               
575                                // pass in the transmission error, and the error in the correction is returned as the last parameter
576                                lat_corr = LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)             //moved from 1D avg SRK 11/2007
577                                data[ii][jj] /= lat_corr                        //divide by the correction factor
578                                //
579                                //
580                                //
581                                // relative errors add in quadrature
582                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
583                                tmp_err = sqrt(tmp_err)
584                               
585                                data_err[ii][jj] = tmp_err
586                               
587//                              solidAngle[ii][jj] = lat_err
588
589                               
590                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
591                        endif
592                       
593                        jj+=1
594                while(jj<pixelsX)
595                ii+=1
596        while(ii<pixelsX)
597       
598        //clean up waves
599//      Killwaves/Z fyy,xx,yy
600       
601        Return(0)
602End
603
604//trig function used by DetCorr()
605Function dc_fx(x,sx,sx3,xcenter)
606        Variable x,sx,sx3,xcenter
607       
608        Variable result
609       
610        result = sx3*tan((x-xcenter)*sx/sx3)
611        Return(result)
612End
613
614//trig function used by DetCorr()
615Function dc_fy(y,sy,sy3,ycenter)
616        Variable y,sy,sy3,ycenter
617       
618        Variable result
619       
620        result = sy3*tan((y-ycenter)*sy/sy3)
621        Return(result)
622End
623
624//trig function used by DetCorr()
625Function dc_fxn(x,sx,sx3,xcenter)
626        Variable x,sx,sx3,xcenter
627       
628        Variable result
629       
630        result = (cos((x-xcenter)*sx/sx3))^2
631        Return(result)
632End
633
634//trig function used by DetCorr()
635Function dc_fym(y,sy,sy3,ycenter)
636        Variable y,sy,sy3,ycenter
637       
638        Variable result
639       
640        result = (cos((y-ycenter)*sy/sy3))^2
641        Return(result)
642End
643
644//distances passed in are in mm
645// dtdist is SDD
646// xd and yd are distances from the beam center to the current pixel
647//
648Function DetEffCorr(lambda,dtdist,xd,yd)
649        Variable lambda,dtdist,xd,yd
650       
651        Variable theta,cosT,ff,stAl,stHe
652       
653        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )
654        cosT = cos(theta)
655       
656        stAl = 0.00967*lambda*0.8               //dimensionless, constants from JGB memo
657        stHe = 0.146*lambda*2.5
658       
659        ff = exp(-stAl/cosT)*(1-exp(-stHe/cosT)) / ( exp(-stAl)*(1-exp(-stHe)) )
660               
661        return(ff)
662End
663
664// DIVIDE the intensity by this correction to get the right answer
665Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err)
666        Variable trans,dtdist,xd,yd,trans_err,&err
667
668        //angle dependent transmission correction
669        Variable uval,arg,cos_th,correction,theta
670       
671        ////this section is the trans_correct() VAX routine
672//      if(trans<0.1)
673//              Print "***transmission is less than 0.1*** and is a significant correction"
674//      endif
675//      if(trans==0)
676//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
677//              trans = 1
678//      endif
679       
680        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel
681       
682        //optical thickness
683        uval = -ln(trans)               //use natural logarithm
684        cos_th = cos(theta)
685        arg = (1-cos_th)/cos_th
686       
687        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
688        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
689        // OR
690        if((uval<0.01) || (cos_th>0.99))       
691                //small arg, approx correction
692                correction= 1-0.5*uval*arg
693        else
694                //large arg, exact correction
695                correction = (1-exp(-uval*arg))/(uval*arg)
696        endif
697
698        Variable tmp
699       
700        if(trans == 1)
701                err = 0         //no correction, no error
702        else
703                //sigT, calculated from the Taylor expansion
704                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
705                tmp *= tmp
706                tmp *= trans_err^2
707                tmp = sqrt(tmp)         //sigT
708               
709                err = tmp
710        endif
711       
712//      Printf "trans error = %g\r",trans_err
713//      Printf "correction = %g +/- %g\r", correction, err
714       
715        //end of transmission/pathlength correction
716
717        return(correction)
718end
719
720//******************
721//direct port of the FORTRAN code for calculating the weighted
722//shifted element to add when beam centers in data headers do not match
723//(indices updated to [0,n-1] indexing rather than (1,n) of fortran
724//
725// as of IGOR 4.0, could be rewritten to pass-by-reference noadd, rather than wave, but the function
726// is so little used, it's not worth the time
727Function ShiftSum(DATA,ip,jp,XSHIFT,YSHIFT,noadd)
728        Wave data
729        Variable ip,jp,xshift,yshift
730        Wave noadd
731//
732//       COMPUTE WEIGHTED OFFSET ELEMENT SUM FOR USE IN SANS DATA
733//       ANALYSIS MODULES.
734//
735// "data" wave passed in is the current contents of the work file
736// sum_val is the return value of the function
737// "noadd" is passed back to the calling function as a one-point wave
738
739        Variable XDELTA,YDELTA,kk,II,JJ,ISHIFT,JSHIFT,sum_val
740        Make/O/N=4 iii,jjj,a
741
742//       -----------------------------------------------------------------
743
744        ISHIFT = trunc(XSHIFT)          // INTEGER PART, trunc gives int closest in dierction of zero
745        XDELTA = XSHIFT - ISHIFT        //FRACTIONAL PART.
746        JSHIFT = trunc(YSHIFT)
747        YDELTA = YSHIFT - JSHIFT
748        II = ip + ISHIFT
749        JJ = jp + JSHIFT
750
751//       SHIFT IS DEFINED AS A VECTOR ANCHORED AT THE STATIONARY CENTER
752//       AND POINTING TO THE MOVABLE CENTER.  THE MOVABLE FIELD IS THUS
753//       ACTUALLY MOVED BY -SHIFT.
754//
755        IF ((XDELTA>= 0) && (YDELTA >= 0))              // CASE I ---- "%&" is "and"
756                III[0] = II
757                JJJ[0] = JJ
758                III[1] = II + 1
759                JJJ[1] = JJ
760                III[2] = II + 1
761                JJJ[2] = JJ + 1
762                III[3] = II
763                JJJ[3] = JJ + 1
764                A[0] = (1. - XDELTA)*(1. - YDELTA)
765                A[1] = XDELTA*(1. - YDELTA)
766                A[2] = XDELTA*YDELTA
767                A[3] = (1. - XDELTA)*YDELTA
768        Endif
769        IF ((XDELTA >= 0) && (YDELTA < 0))              // CASE II.
770                III[0] = II
771                JJJ[0] = JJ
772                III[1] = II
773                JJJ[1] = JJ - 1
774                III[2] = II + 1
775                JJJ[2] = JJ - 1
776                III[3] = II + 1
777                JJJ[3] = JJ
778                A[0] = (1. - XDELTA)*(1. + YDELTA)
779                A[1] = (1. - XDELTA)*(-YDELTA)
780                A[2] = XDELTA*(-YDELTA)
781                A[3] = XDELTA*(1. + YDELTA)
782        Endif
783        IF ((XDELTA < 0) && (YDELTA >= 0))              // CASE III.
784                III[0] = II
785                JJJ[0] = JJ
786                III[1] = II
787                JJJ[1] = JJ + 1
788                III[2] = II - 1
789                JJJ[2] = JJ + 1
790                III[3] = II - 1
791                JJJ[3] = JJ
792                A[0] = (1. + XDELTA)*(1 - YDELTA)
793                A[1] = (1. + XDELTA)*YDELTA
794                A[2] = -XDELTA*YDELTA
795                A[3] = -XDELTA*(1. - YDELTA)
796        Endif
797        IF ((XDELTA < 0) && (YDELTA < 0))               //CASE IV.
798                III[0] = II
799                JJJ[0] = JJ
800                III[1] = II - 1
801                JJJ[1] = JJ
802                III[2] = II - 1
803                JJJ[2] = JJ - 1
804                III[3] = II
805                JJJ[3] = JJ - 1
806                A[0] = (1. + XDELTA)*(1. + YDELTA)
807                A[1] = -XDELTA*(1. + YDELTA)
808                A[2] = (-XDELTA)*(-YDELTA)
809                A[3] = (1. + XDELTA)*(-YDELTA)
810        Endif
811
812        NVAR pixelsX = root:myGlobals:gNPixelsX
813        NVAR pixelsY = root:myGlobals:gNPixelsY
814//check to see if iii[0],jjj[0] are valid detector elements, in [0,127]
815//if not set noadd[0] to 1, to let calling routine know NOT to add
816//        CALL TESTIJ(III(1),JJJ(1),OKIJ)
817        NOADD[0] = 0
818        if( (iii[0]<0) || (iii[0]>(pixelsX-1)) )
819                noadd[0] = 1
820        endif
821        if((jjj[0]<0) || (jjj[0]>(pixelsY-1)) )
822                noadd[0] = 1
823        endif
824       
825
826       
827        sum_val = 0.
828        kk = 0
829        Do
830                IF(JJJ[kk] == pixelsX)
831                        //do nothing
832                else
833                        sum_val += A[kk]*DATA[III[kk]][JJJ[kk]]
834                endif
835                kk+=1
836        while(kk<4)
837       
838        //clean up waves
839        KillWaves/z iii,jjj,a
840       
841        RETURN (sum_val)
842       
843End             //function ShiftSum
844
845//************************
846//unused testing procedure, may not be up-to-date with other procedures
847//check before re-implementing
848//
849Proc DIV_a_Workfile(type)
850        String type
851        Prompt type,"WORK data type",popup,"COR;SAM;EMP;BGD"
852       
853        //macro will take whatever is in SELECTED folder and DIVide it by the current
854        //contents of the DIV folder - the function will check for existence
855        //before proceeding
856       
857        Variable err
858        err = Divide_work(type)         //returns err = 1 if data doesn't exist in specified folders
859       
860        if(err)
861                Abort "error in Divide_work"
862        endif
863       
864        //contents are always dumped to CAL
865        type = "CAL"
866       
867        String newTitle = "WORK_"+type
868        DoWindow/F SANS_Data
869        DoWindow/T SANS_Data, newTitle
870        KillStrings/Z newTitle
871       
872        //need to update the display with "data" from the correct dataFolder
873        //reset the current displaytype to "type"
874        String/G root:myGlobals:gDataDisplayType=Type
875       
876        fRawWindowHook()
877       
878End
879
880//function will divide the contents of "type" folder with the contents of
881//the DIV folder
882// all data is converted to linear scale for the calculation
883//
884Function Divide_work(type)
885        String type
886       
887        //check for existence of data in type and DIV
888        // if the desired workfile doesn't exist, let the user know, and abort
889        String destPath=""
890
891        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
892                Print "There is no work file in "+type+"--Aborting"
893                Return(1)               //error condition
894        Endif
895        //check for DIV
896        // if the DIV workfile doesn't exist, let the user know,and abort
897
898        if(WaveExists($"root:Packages:NIST:DIV:data") == 0)
899                Print "There is no work file in DIV --Aborting"
900                Return(1)               //error condition
901        Endif
902        //files exist, proceed
903       
904        //check for log-scaling of the "DIV" data and adjust if necessary
905        ConvertFolderToLinearScale("DIV")
906       
907        //copy type information to CAL, wiping out the old contents of the CAL folder first
908       
909        //destPath = "root:Packages:NIST:CAL"
910        //SetDataFolder destPath
911        //KillWaves/A/Z                 //get rid of the old data in CAL folder
912
913        //check for log-scaling of the "type" data and adjust if necessary
914        ConvertFolderToLinearScale(type)
915        //then continue
916
917        //copy from current dir (type)=destPath to CAL, overwriting CAL contents
918        destPath = "root:Packages:NIST:" + type
919        Duplicate/O $(destPath + ":data"),$"root:Packages:NIST:CAL:data"
920        Duplicate/O $(destPath + ":linear_data"),$"root:Packages:NIST:CAL:linear_data"
921        Duplicate/O $(destPath + ":linear_data_error"),$"root:Packages:NIST:CAL:linear_data_error"
922//      Duplicate/O $(destPath + ":vlegend"),$"root:Packages:NIST:CAL:vlegend"
923        Duplicate/O $(destPath + ":textread"),$"root:Packages:NIST:CAL:textread"
924        Duplicate/O $(destPath + ":integersread"),$"root:Packages:NIST:CAL:integersread"
925        Duplicate/O $(destPath + ":realsread"),$"root:Packages:NIST:CAL:realsread"
926        //need to save a copy of filelist string too (from the current type folder)
927        SVAR oldFileList = $(destPath + ":fileList")
928
929        //now switch to reference waves in CAL folder
930        destPath = "root:Packages:NIST:CAL"
931        //make appropriate wave references
932        Wave data=$(destPath + ":linear_data")                                  // these wave references point to the data in CAL
933//      Wave data_err=$(destPath + ":linear_data_err")                                  // these wave references point to the data in CAL
934        Wave data_copy=$(destPath + ":data")                                    // these wave references point to the data in CAL
935        Wave/t textread=$(destPath + ":textread")                       //that are to be directly operated on
936        Wave integersread=$(destPath + ":integersread")
937        Wave realsread=$(destPath + ":realsread")
938        Variable/G $(destPath + ":gIsLogScale")=0                       //make new flag in CAL folder, data is linear scale
939        //need to copy filelist string too
940        String/G $(destPath + ":fileList") = oldFileList
941
942        Wave div_data = $"root:Packages:NIST:DIV:data"          //hard-wired in....
943        //do the division, changing data in CAL
944        data /= div_data
945       
946//      data_err /= div_data
947       
948        // keep "data" in sync with linear_data
949        data_copy = data
950       
951        //update CAL header
952        textread[1] = date() + " " + time()             //date + time stamp
953       
954        Return(0)
955End
956
957//test procedure, not called anymore
958Proc AbsoluteScaling(type,c0,c1,c2,c3,c4,c5)
959        String type
960        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0
961        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
962        Prompt c0, "Sample Transmission"
963        Prompt c1, "Sample Thickness (cm)"
964        Prompt c2, "Standard Transmission"
965        Prompt c3, "Standard Thickness (cm)"
966        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
967        Prompt c5, "Standard Cross-Section (cm-1)"
968
969        Variable err
970        //call the function to do the math
971        //data from "type" will be scaled and deposited in ABS
972        err = Absolute_Scale(type,c0,c1,c2,c3,c4,c5)
973       
974        if(err)
975                Abort "Error in Absolute_Scale()"
976        endif
977       
978        //contents are always dumped to ABS
979        type = "ABS"
980       
981        String newTitle = "WORK_"+type
982        DoWindow/F SANS_Data
983        DoWindow/T SANS_Data, newTitle
984        KillStrings/Z newTitle
985       
986        //need to update the display with "data" from the correct dataFolder
987        //reset the current displaytype to "type"
988        String/G root:myGlobals:gDataDisplayType=Type
989       
990        fRawWindowHook()
991       
992End
993
994//s_ is the standard
995//w_ is the "work" file
996//both are work files and should already be normalized to 10^8 monitor counts
997Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
998        String type
999        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
1000       
1001        //convert the "type" data to absolute scale using the given standard information
1002        //copying the "type" waves to ABS
1003       
1004        //check for existence of data, rescale to linear if needed
1005        String destPath
1006        //check for "type"
1007        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
1008                Print "There is no work file in "+type+"--Aborting"
1009                Return(1)               //error condition
1010        Endif
1011        //check for log-scaling of the "type" data and adjust if necessary
1012        destPath = "root:Packages:NIST:"+Type
1013        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
1014        if(gIsLogScale)
1015                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
1016                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
1017        endif
1018       
1019        //copy "oldtype" information to ABS
1020        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
1021        //copy over the waves data,vlegend,text,integers,reals(read)
1022
1023        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
1024        //copy from current dir (type) to ABS, defined by destPath
1025        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
1026        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
1027        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
1028//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
1029        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
1030        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
1031        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
1032        //need to save a copy of filelist string too (from the current type folder)
1033        SVAR oldFileList = $(oldType + ":fileList")
1034        //need to copy filelist string too
1035        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
1036       
1037        //now switch to ABS folder
1038        //make appropriate wave references
1039        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
1040        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
1041        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
1042        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
1043        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
1044        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
1045        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
1046       
1047        //do the actual absolute scaling here, modifying the data in ABS
1048        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1049       
1050        w_moncount = realsread[0]               //monitor count in "type"
1051        if(w_moncount == 0)
1052                //zero monitor counts will give divide by zero ---
1053                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1054                Return(1)               //report error
1055        Endif
1056       
1057        //calculate scale factor
1058        Variable scale,trans_err
1059        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
1060        s2 = s_thick/w_thick
1061        s3 = s_trans/w_trans
1062        s4 = s_cross/s_izero
1063       
1064        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1065       
1066        data *= s1*s2*s3*s4
1067       
1068        scale = s1*s2*s3*s4
1069        trans_err = realsRead[41]
1070       
1071//      print scale
1072//      print data[0][0]
1073       
1074        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1075
1076//      print data_err[0][0]
1077       
1078// keep "data" in sync with linear_data
1079        data_copy = data
1080       
1081        //********* 15APR02
1082        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
1083        // footing (zero atten) before doing the subtraction.
1084        //
1085        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
1086       
1087        //update the ABS header information
1088        textread[1] = date() + " " + time()             //date + time stamp
1089       
1090        Return (0) //no error
1091End
1092
1093//*************
1094// start of section of functions used for workfile math panel
1095//*************
1096
1097
1098//function will copy the contents of oldtype folder to newtype folder
1099//converted to linear scale before copying
1100//******data in newtype is overwritten********
1101//
1102Function CopyWorkContents(oldtype,newtype)
1103        String oldtype,newtype
1104       
1105        //check for existence of data in oldtype
1106        // if the desired workfile doesn't exist, let the user know, and abort
1107        String destPath=""
1108        if(WaveExists($("root:Packages:NIST:"+oldType + ":data")) == 0)
1109                Print "There is no work file in "+oldtype+"--Aborting"
1110                Return(1)               //error condition
1111        Endif
1112       
1113        //check for log-scaling of the "type" data and adjust if necessary
1114        ConvertFolderToLinearScale(oldtype)
1115        Fix_LogLinButtonState(0)                //make sure the button reflects the new linear scaling
1116        //then continue
1117
1118        //copy from current dir (type)=destPath to newtype, overwriting newtype contents
1119        destPath = "root:Packages:NIST:" + oldtype
1120        Duplicate/O $(destPath + ":data"),$("root:Packages:NIST:"+newtype+":data")
1121        Duplicate/O $(destPath + ":linear_data"),$("root:Packages:NIST:"+newtype+":linear_data")
1122        Duplicate/O $(destPath + ":linear_data_error"),$("root:Packages:NIST:"+newtype+":linear_data_error")
1123        Duplicate/O $(destPath + ":textread"),$("root:Packages:NIST:"+newtype+":textread")
1124        Duplicate/O $(destPath + ":integersread"),$("root:Packages:NIST:"+newtype+":integersread")
1125        Duplicate/O $(destPath + ":realsread"),$("root:Packages:NIST:"+newtype+":realsread")
1126        //
1127        // be sure to get rid of the linear_data if it exists in the destination folder
1128//      KillWaves/Z $("root:Packages:NIST:"+newtype+":linear_data")
1129
1130        //need to save a copy of filelist string too (from the current type folder)
1131        SVAR oldFileList = $(destPath + ":fileList")
1132
1133        //now switch to reference waves in newtype folder
1134        destPath = "root:Packages:NIST:"+newtype
1135        Variable/G $(destPath + ":gIsLogScale")=0                       //make new flag in newtype folder, data is linear scale
1136        //need to copy filelist string too
1137        String/G $(destPath + ":fileList") = oldFileList
1138
1139        Return(0)
1140End
1141
1142//Entry procedure from main panel
1143//
1144Proc CopyWorkFolder(oldType,newType)
1145        String oldType,newType
1146        Prompt oldType,"Source WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
1147        Prompt newType,"Destination WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
1148//      Prompt oldType,"Source WORK data type",popup,"AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1149//      Prompt newType,"Destination WORK data type",popup,"AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1150
1151        // data folder "old" will be copied to "new" (and will overwrite)
1152        CopyWorkContents(oldtype,newtype)
1153End
1154
1155//initializes datafolder and constants necessary for the workfilemath panel
1156//
1157Proc Init_WorkMath()
1158        //create the data folder
1159        //String str = "AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1160        String str = "File_1;File_2;Result;"
1161        NewDataFolder/O/S root:Packages:NIST:WorkMath
1162        String/G gFolderList=str
1163        Variable ii=0,num=itemsinlist(str)
1164        do
1165                Execute "NewDataFolder/O "+StringFromList(ii, str ,";")
1166                ii+=1
1167        while(ii<num)
1168        Variable/G const1=1,const2=1
1169       
1170        SetDataFolder root:
1171End
1172
1173//entry procedure to invoke the workfilemath panel, initializes if needed
1174//
1175Proc Show_WorkMath_Panel()
1176        DoWindow/F WorkFileMath
1177        if(V_flag==0)
1178                Init_WorkMath()
1179                WorkFileMath()
1180        Endif
1181End
1182
1183//attempts to perform the selected math operation based on the selections in the panel
1184// aborts with an error condition if operation can't be completed
1185// or puts the final answer in the Result folder, and displays the selected data
1186//
1187Function WorkMath_DoIt_ButtonProc(ctrlName) : ButtonControl
1188        String ctrlName
1189
1190        String str1,str2,oper,dest = "Result"
1191        String pathStr,workMathStr="WorkMath:"
1192       
1193        //get the panel selections (these are the names of the files on disk)
1194        PathInfo catPathName
1195        pathStr=S_Path
1196        ControlInfo popup0
1197        str1=S_Value
1198        ControlInfo popup1
1199        str2=S_Value
1200        ControlInfo popup2
1201        oper=S_Value
1202       
1203        //check that something has been selected for operation and destination
1204        if(cmpstr(oper,"Operation")==0)
1205                Abort "Select a math operand from the popup"
1206        Endif
1207
1208        //constants from globals
1209        NVAR const1=root:Packages:NIST:WorkMath:const1
1210        NVAR const2=root:Packages:NIST:WorkMath:const2
1211        Printf "(%g)%s %s (%g)%s = %s\r", const1,str1,oper,const2,str2,dest
1212        //check for proper folders (all 3 must be different)
1213       
1214        //load the data in here...
1215        //set #1
1216        Load_NamedASC_File(pathStr+str1,workMathStr+"File_1")
1217       
1218        NVAR pixelsX = root:myGlobals:gNPixelsX         //OK, location is correct
1219        NVAR pixelsY = root:myGlobals:gNPixelsY
1220       
1221        WAVE/Z data1=$("root:Packages:NIST:"+workMathStr+"File_1:linear_data")
1222        WAVE/Z err1=$("root:Packages:NIST:"+workMathStr+"File_1:linear_data_error")
1223       
1224        // set # 2
1225        If(cmpstr(str2,"UNIT MATRIX")==0)
1226                Make/O/N=(pixelsX,pixelsY) root:Packages:NIST:WorkMath:data             //don't put in File_2 folder
1227                Wave/Z data2 =  root:Packages:NIST:WorkMath:data                        //it's not real data!
1228                data2=1
1229                Duplicate/O data2 err2
1230                err2 = 0
1231        else
1232                //Load set #2
1233                Load_NamedASC_File(pathStr+str2,workMathStr+"File_2")
1234                WAVE/Z data2=$("root:Packages:NIST:"+workMathStr+"File_2:linear_data")
1235                WAVE/Z err2=$("root:Packages:NIST:"+workMathStr+"File_2:linear_data_error")
1236        Endif
1237
1238        ///////
1239       
1240        //now that we know that data exists, convert each of the operands to linear scale
1241//      ConvertFolderToLinearScale(workMathStr+"File_1")
1242//      If(cmpstr(str2,"UNIT MATRIX")!=0)
1243//              ConvertFolderToLinearScale(workMathStr+"File_2")                //don't need to convert unit matrix to linear
1244//      endif
1245
1246        //copy contents of str1 folder to dest and create the wave ref (it will exist)
1247        CopyWorkContents(workMathStr+"File_1",workMathStr+dest)
1248        WAVE/Z destData=$("root:Packages:NIST:"+workMathStr+dest+":linear_data")
1249        WAVE/Z destData_log=$("root:Packages:NIST:"+workMathStr+dest+":data")
1250        WAVE/Z destErr=$("root:Packages:NIST:"+workMathStr+dest+":linear_data_error")
1251       
1252        //dispatch
1253        strswitch(oper)
1254                case "*":               //multiplication
1255                        destData = const1*data1 * const2*data2
1256                        destErr = const1^2*const2^2*(err1^2*data2^2 + err2^2*data1^2)
1257                        destErr = sqrt(destErr)
1258                        break   
1259                case "_":               //subtraction
1260                        destData = const1*data1 - const2*data2
1261                        destErr = const1^2*err1^2 + const2^2*err2^2
1262                        destErr = sqrt(destErr)
1263                        break
1264                case "/":               //division
1265                        destData = (const1*data1) / (const2*data2)
1266                        destErr = const1^2/const2^2*(err1^2/data2^2 + err2^2*data1^2/data2^4)
1267                        destErr = sqrt(destErr)
1268                        break
1269                case "+":               //addition
1270                        destData = const1*data1 + const2*data2
1271                        destErr = const1^2*err1^2 + const2^2*err2^2
1272                        destErr = sqrt(destErr)
1273                        break                   
1274        endswitch
1275       
1276        destData_log = log(destData)            //for display
1277        //show the result
1278        WorkMath_Display_PopMenuProc("",0,"Result")
1279       
1280        PopupMenu popup4 win=WorkFileMath,mode=3                //3rd item selected == Result
1281End
1282
1283// closes the panel and kills the data folder when done
1284//
1285Function WorkMath_Done_ButtonProc(ctrlName) : ButtonControl
1286        String ctrlName
1287       
1288        DoAlert 1,"Closing the panel will kill all of the data you have loaded into memory. Do you want to continue?"
1289        If(V_Flag==2)
1290                return(0)               //don't do anything
1291        Endif
1292        //kill the panel
1293        DoWindow/K WorkFileMath
1294        //wipe out the data folder of globals
1295        SVAR dataType=root:myGlobals:gDataDisplayType
1296        if(strsearch(dataType, "WorkMath", 0 ) != -1)           //kill the SANS_Data graph if needed
1297                DoWindow/K SANS_Data
1298        Endif
1299        KillDataFolder root:Packages:NIST:WorkMath
1300End
1301
1302// loads data into the specified folder
1303//
1304// currently unused since button has been removed from panel
1305//
1306Function WorkMath_Load_ButtonProc(ctrlName) : ButtonControl
1307        String ctrlName
1308       
1309        String destStr=""
1310        SVAR folderList=root:Packages:NIST:WorkMath:gFolderList
1311        Prompt destStr,"Select the destination folder",popup,folderList
1312        DoPrompt "Folder for ASC Load",destStr
1313       
1314        if(V_flag==1)
1315                return(1)               //user abort, do nothing
1316        Endif
1317       
1318        String destFolder = "WorkMath:"+destStr
1319       
1320        Load_ASC_File("Pick the ASC file",destFolder)
1321End
1322
1323// changes the display of the SANS_Data window based on popped data type
1324// first loads in the data from the File1 or File2 popup as needed
1325// then displays the selcted dataset, if it exists
1326// makes use of procedure from DisplayUtils
1327//
1328// - Always replaces File1 or File2 with fresh data from disk
1329//
1330Function WorkMath_Display_PopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
1331        String ctrlName
1332        Variable popNum
1333        String popStr
1334       
1335        String folder="WorkMath:",pathStr,str1
1336
1337        PathInfo catPathName
1338        pathStr=S_Path
1339       
1340        //if display result, just do it and return
1341        if(cmpstr(popStr,"Result")==0)
1342                Execute "ChangeDisplay(\""+folder+popstr+"\")"
1343                return(0)
1344        endif
1345        // if file1 or file2, load in the data and display
1346        if(cmpstr(popStr,"File_1")==0)
1347                ControlInfo/W=WorkFileMath popup0
1348                str1 = S_Value
1349        Endif
1350        if(cmpstr(popStr,"File_2")==0)
1351                ControlInfo/W=WorkFileMath popup1
1352                str1 = S_Value
1353        Endif
1354        //don't load or display the unit matrix
1355        Print str1
1356        if(cmpstr(str1,"UNIT MATRIX")!=0)
1357                Load_NamedASC_File(pathStr+str1,folder+popStr)
1358                //change the display
1359                Execute "ChangeDisplay(\""+folder+popstr+"\")"
1360        endif
1361        return(0)       
1362End
1363
1364//simple panel to do workfile arithmetic
1365//
1366Proc WorkFileMath()
1367        PauseUpdate; Silent 1           // building window...
1368        NewPanel /W=(610,211,880,490)/K=2 as "Work File Math"           // replace /K=2 flag
1369        DoWindow/C WorkFileMath
1370        ModifyPanel cbRGB=(47802,54484,6682)
1371        ModifyPanel fixedSize=1
1372        SetDrawLayer UserBack
1373        DrawLine 6,166,214,166
1374        SetDrawEnv fstyle= 4
1375        DrawText 10,34,"File #1:"
1376        SetDrawEnv fstyle= 4
1377        DrawText 13,129,"File #2:"
1378        DrawText 78,186,"= Result"
1379        Button button0,pos={28,245},size={50,20},proc=WorkMath_DoIt_ButtonProc,title="Do It"
1380        Button button0,help={"Performs the specified arithmetic"}
1381        Button button1,pos={183,245},size={50,20},proc=WorkMath_Done_ButtonProc,title="Done"
1382        Button button1,help={"Closes the panel"}
1383//      Button button2,pos={30,8},size={110,20},proc=WorkMath_Load_ButtonProc,title="Load ASC Data"
1384//      Button button2,help={"Loads ASC data files into the specified folder"}
1385        Button button3,pos={205,8},size={25,20},proc=ShowWorkMathHelp,title="?"
1386        Button button3,help={"Show help file for math operations on 2-D data sets"}
1387        SetVariable setvar0,pos={9,46},size={70,15},title=" "
1388        SetVariable setvar0,limits={-Inf,Inf,0},value= root:Packages:NIST:WorkMath:const1
1389        SetVariable setvar0,help={"Multiplicative constant for the first dataset"}
1390        PopupMenu popup0,pos={89,44},size={84,20},title="X  "
1391        PopupMenu popup0,mode=1,popvalue="1st Set",value= ASC_FileList()
1392        PopupMenu popup0,help={"Selects the first dataset for the operation"}
1393        PopupMenu popup1,pos={93,136},size={89,20},title="X  "
1394        PopupMenu popup1,mode=1,popvalue="2nd Set",value= "UNIT MATRIX;"+ASC_FileList()
1395        PopupMenu popup1,help={"Selects the second dataset for the operation"}
1396        PopupMenu popup2,pos={50,82},size={70,20},title="Operator  "
1397        PopupMenu popup2,mode=3,popvalue="Operation",value= #"\"+;_;*;/;\""
1398        PopupMenu popup2,help={"Selects the mathematical operator"}
1399        SetVariable setvar1,pos={13,139},size={70,15},title=" "
1400        SetVariable setvar1,limits={-Inf,Inf,0},value= root:Packages:NIST:WorkMath:const2
1401        SetVariable setvar1,help={"Multiplicative constant for the second dataset"}
1402//      PopupMenu popup3,pos={27,167},size={124,20},title=" = Destination"
1403//      PopupMenu popup3,mode=1,popvalue="Destination",value= root:Packages:NIST:WorkMath:gFolderList
1404//      PopupMenu popup3,help={"Selects the destination folder"}
1405        PopupMenu popup4,pos={55,204},size={103,20},proc=WorkMath_Display_PopMenuProc,title="Display"
1406        PopupMenu popup4,mode=3,value= "File_1;File_2;Result;"
1407        PopupMenu popup4,help={"Displays the data in the specified folder"}
1408EndMacro
1409
1410//jump to the help topic
1411Function ShowWorkMathHelp(ctrlName) : ButtonControl
1412        String ctrlName
1413        DisplayHelpTopic/Z/K=1 "SANS Data Reduction Tutorial[2-D Work File Arithmetic]"
1414        if(V_flag !=0)
1415                DoAlert 0,"The SANS Data Reduction Tutorial Help file could not be found"
1416        endif
1417End
1418
1419//utility function to clear the contents of a data folder
1420//won't clear data that is in use -
1421//
1422Function ClearDataFolder(type)
1423        String type
1424       
1425        SetDataFolder $("root:Packages:NIST:"+type)
1426        KillWaves/a/z
1427        KillStrings/a/z
1428        KillVariables/a/z
1429        SetDataFolder root:
1430End
1431
1432
1433
1434//fileStr must be the FULL path and filename on disk
1435//destFolder is the path to the Igor folder where the data is to be deposited
1436// - "Packages:NIST:WorkMath:File_1" for example, compatible with SANS_Data display type
1437//
1438Function Load_NamedASC_File(fileStr,destFolder)
1439        String fileStr,destFolder
1440
1441        Variable refnum
1442       
1443        //read in the data
1444        ReadASCData(fileStr,destFolder)
1445
1446        //the calling macro must change the display type
1447        String/G root:myGlobals:gDataDisplayType=destFolder
1448       
1449        FillFakeHeader_ASC(destFolder)          //uses info on the panel, if available
1450
1451        //data is displayed here, and needs header info
1452       
1453        fRawWindowHook()
1454       
1455        return(0)
1456End
1457
1458
1459//function called by the main entry procedure (the load button)
1460//loads data into the specified folder, and wipes out what was there
1461//
1462//aborts if no file was selected
1463//
1464// reads the data in if all is OK
1465//
1466// currently unused, as load button has been replaced
1467//
1468Function Load_ASC_File(msgStr,destFolder)
1469        String msgStr,destFolder
1470
1471        String filename="",pathStr=""
1472        Variable refnum
1473
1474        //check for the path
1475        PathInfo catPathName
1476        If(V_Flag==1)           //      /D does not open the file
1477                Open/D/R/T="????"/M=(msgStr)/P=catPathName refNum
1478        else
1479                Open/D/R/T="????"/M=(msgStr) refNum
1480        endif
1481        filename = S_FileName           //get the filename, or null if canceled from dialog
1482        if(strlen(filename)==0)
1483                //user cancelled, abort
1484                SetDataFolder root:
1485                Abort "No file selected, action aborted"
1486        Endif
1487       
1488        //read in the data
1489        ReadASCData(filename,destFolder)
1490
1491        //the calling macro must change the display type
1492        String/G root:myGlobals:gDataDisplayType=destFolder
1493       
1494        FillFakeHeader_ASC(destFolder)          //uses info on the panel, if available
1495
1496        //data is displayed here, and needs header info
1497       
1498        fRawWindowHook()
1499       
1500        return(0)
1501End
1502
1503
1504
1505//function called by the popups to get a file list of data that can be sorted
1506// this procedure simply removes the raw data files from the string - there
1507//can be lots of other junk present, but this is very fast...
1508//
1509// -- simplified to do in a single call -- extension must be .ASC
1510// - so far, only used in WorkFileMath popups, which require .ASC
1511//
1512Function/S ASC_FileList()
1513
1514        String list="",newList="",item=""
1515        Variable num,ii
1516       
1517        //check for the path
1518        PathInfo catPathName
1519        if(V_Flag==0)
1520                DoAlert 0, "Data path does not exist - pick the data path from the button on the main panel"
1521                Return("")
1522        Endif
1523
1524        list = IndexedFile(catpathName,-1,".ASC")
1525
1526
1527//      list = IndexedFile(catpathName,-1,"????")
1528//     
1529//      list = RemoveFromList(ListMatch(list,"*.SA1*",";"), list, ";", 0)
1530//      list = RemoveFromList(ListMatch(list,"*.SA2*",";"), list, ";", 0)
1531//      list = RemoveFromList(ListMatch(list,"*.SA3*",";"), list, ";", 0)
1532//      list = RemoveFromList(ListMatch(list,".*",";"), list, ";", 0)
1533//      list = RemoveFromList(ListMatch(list,"*.pxp",";"), list, ";", 0)
1534//      list = RemoveFromList(ListMatch(list,"*.DIV",";"), list, ";", 0)
1535//      list = RemoveFromList(ListMatch(list,"*.GSP",";"), list, ";", 0)
1536//      list = RemoveFromList(ListMatch(list,"*.MASK",";"), list, ";", 0)
1537//
1538//      //remove VAX version numbers
1539//      list = RemoveVersNumsFromList(List)
1540        //sort
1541        newList = SortList(List,";",0)
1542
1543        return newlist
1544End
Note: See TracBrowser for help on using the repository browser.