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

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

Changed 2-yukawa to have better names for the parameters, a macro to plot g(r), and to return S(q)=1000 if the calculation fails, to be more obvious.

Added a Corrlation matrix function (in Wrapper.ipf), but it's not linked anywhere yet, since I don't know yet how best to display it.

Simplified the function to list .ASC files

Added comments to the polarization panel procedures.
File size: 50.1 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                        data[ii][jj] *= xy*ratio
538                        solidAngle[ii][jj] = xy*ratio           //testing only 
539                        data_err[ii][jj] *= xy*ratio                    //error propagation assumes that SA and Jacobian are exact, so simply scale error
540                       
541                       
542                        // correction factor for detector efficiency JBG memo det_eff_cor2.doc 3/20/07
543                        // correction inserted 11/2007 SRK
544                        // large angle detector efficiency is >= 1 and will "bump up" the measured value at the highest angles
545                        // so divide here to get the correct answer (5/22/08 SRK)
546                        if(doEfficiency)
547#if (exists("ILL_D22")==6)
548                                data[ii][jj]  /= DetEffCorrILL(lambda,dtdist,xd)                //tube-by-tube corrections
549                                data_err[ii][jj] /= DetEffCorrILL(lambda,dtdist,xd)                     //assumes correction is exact
550//                solidAngle[ii][jj] = DetEffCorrILL(lambda,dtdist,xd)
551#else
552                                data[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
553                                data_err[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)
554//                              solidAngle[ii][jj] /= DetEffCorr(lambda,dtdist,xd,yd)           //testing only
555#endif
556                        endif
557                       
558                        // large angle transmission calculation is <= 1 and will "bump down" the measured value at the highest angles
559                        // so divide here to get the correct answer
560                        if(doTrans)
561                       
562                                if(trans<0.1 && ii==0 && jj==0)
563                                        Print "***transmission is less than 0.1*** and is a significant correction"
564                                endif
565                               
566                                if(trans==0)
567                                        if(ii==0 && jj==0)
568                                                Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
569                                        endif
570                                        trans = 1
571                                endif
572                               
573                                // pass in the transmission error, and the error in the correction is returned as the last parameter
574                                lat_corr = LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,lat_err)             //moved from 1D avg SRK 11/2007
575                                data[ii][jj] /= lat_corr                        //divide by the correction factor
576                                //
577                                //
578                                //
579                                // relative errors add in quadrature
580                                tmp_err = (data_err[ii][jj]/lat_corr)^2 + (lat_err/lat_corr)^2*data[ii][jj]*data[ii][jj]/lat_corr^2
581                                tmp_err = sqrt(tmp_err)
582                               
583                                data_err[ii][jj] = tmp_err
584                               
585                                solidAngle[ii][jj] = lat_err
586
587                               
588                                //solidAngle[ii][jj] = LargeAngleTransmissionCorr(trans,dtdist,xd,yd)           //testing only
589                        endif
590                       
591                        jj+=1
592                while(jj<pixelsX)
593                ii+=1
594        while(ii<pixelsX)
595       
596        //clean up waves
597//      Killwaves/Z fyy,xx,yy
598       
599        Return(0)
600End
601
602//trig function used by DetCorr()
603Function dc_fx(x,sx,sx3,xcenter)
604        Variable x,sx,sx3,xcenter
605       
606        Variable result
607       
608        result = sx3*tan((x-xcenter)*sx/sx3)
609        Return(result)
610End
611
612//trig function used by DetCorr()
613Function dc_fy(y,sy,sy3,ycenter)
614        Variable y,sy,sy3,ycenter
615       
616        Variable result
617       
618        result = sy3*tan((y-ycenter)*sy/sy3)
619        Return(result)
620End
621
622//trig function used by DetCorr()
623Function dc_fxn(x,sx,sx3,xcenter)
624        Variable x,sx,sx3,xcenter
625       
626        Variable result
627       
628        result = (cos((x-xcenter)*sx/sx3))^2
629        Return(result)
630End
631
632//trig function used by DetCorr()
633Function dc_fym(y,sy,sy3,ycenter)
634        Variable y,sy,sy3,ycenter
635       
636        Variable result
637       
638        result = (cos((y-ycenter)*sy/sy3))^2
639        Return(result)
640End
641
642//distances passed in are in mm
643// dtdist is SDD
644// xd and yd are distances from the beam center to the current pixel
645//
646Function DetEffCorr(lambda,dtdist,xd,yd)
647        Variable lambda,dtdist,xd,yd
648       
649        Variable theta,cosT,ff,stAl,stHe
650       
651        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )
652        cosT = cos(theta)
653       
654        stAl = 0.00967*lambda*0.8               //dimensionless, constants from JGB memo
655        stHe = 0.146*lambda*2.5
656       
657        ff = exp(-stAl/cosT)*(1-exp(-stHe/cosT)) / ( exp(-stAl)*(1-exp(-stHe)) )
658               
659        return(ff)
660End
661
662// DIVIDE the intensity by this correction to get the right answer
663Function LargeAngleTransmissionCorr(trans,dtdist,xd,yd,trans_err,err)
664        Variable trans,dtdist,xd,yd,trans_err,&err
665
666        //angle dependent transmission correction
667        Variable uval,arg,cos_th,correction,theta
668       
669        ////this section is the trans_correct() VAX routine
670//      if(trans<0.1)
671//              Print "***transmission is less than 0.1*** and is a significant correction"
672//      endif
673//      if(trans==0)
674//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
675//              trans = 1
676//      endif
677       
678        theta = atan( (sqrt(xd^2 + yd^2))/dtdist )              //theta at the input pixel
679       
680        //optical thickness
681        uval = -ln(trans)               //use natural logarithm
682        cos_th = cos(theta)
683        arg = (1-cos_th)/cos_th
684       
685        // a Taylor series around uval*arg=0 only needs about 4 terms for very good accuracy
686        //                      correction= 1 - 0.5*uval*arg + (uval*arg)^2/6 - (uval*arg)^3/24 + (uval*arg)^4/120
687        // OR
688        if((uval<0.01) || (cos_th>0.99))       
689                //small arg, approx correction
690                correction= 1-0.5*uval*arg
691        else
692                //large arg, exact correction
693                correction = (1-exp(-uval*arg))/(uval*arg)
694        endif
695
696        Variable tmp
697       
698        if(trans == 1)
699                err = 0         //no correction, no error
700        else
701                //sigT, calculated from the Taylor expansion
702                tmp = (1/trans)*(arg/2-arg^2/3*uval+arg^3/8*uval^2-arg^4/30*uval^3)
703                tmp *= tmp
704                tmp *= trans_err^2
705                tmp = sqrt(tmp)         //sigT
706               
707                err = tmp
708        endif
709       
710//      Printf "trans error = %g\r",trans_err
711//      Printf "correction = %g +/- %g\r", correction, err
712       
713        //end of transmission/pathlength correction
714
715        return(correction)
716end
717
718//******************
719//direct port of the FORTRAN code for calculating the weighted
720//shifted element to add when beam centers in data headers do not match
721//(indices updated to [0,n-1] indexing rather than (1,n) of fortran
722//
723// as of IGOR 4.0, could be rewritten to pass-by-reference noadd, rather than wave, but the function
724// is so little used, it's not worth the time
725Function ShiftSum(DATA,ip,jp,XSHIFT,YSHIFT,noadd)
726        Wave data
727        Variable ip,jp,xshift,yshift
728        Wave noadd
729//
730//       COMPUTE WEIGHTED OFFSET ELEMENT SUM FOR USE IN SANS DATA
731//       ANALYSIS MODULES.
732//
733// "data" wave passed in is the current contents of the work file
734// sum_val is the return value of the function
735// "noadd" is passed back to the calling function as a one-point wave
736
737        Variable XDELTA,YDELTA,kk,II,JJ,ISHIFT,JSHIFT,sum_val
738        Make/O/N=4 iii,jjj,a
739
740//       -----------------------------------------------------------------
741
742        ISHIFT = trunc(XSHIFT)          // INTEGER PART, trunc gives int closest in dierction of zero
743        XDELTA = XSHIFT - ISHIFT        //FRACTIONAL PART.
744        JSHIFT = trunc(YSHIFT)
745        YDELTA = YSHIFT - JSHIFT
746        II = ip + ISHIFT
747        JJ = jp + JSHIFT
748
749//       SHIFT IS DEFINED AS A VECTOR ANCHORED AT THE STATIONARY CENTER
750//       AND POINTING TO THE MOVABLE CENTER.  THE MOVABLE FIELD IS THUS
751//       ACTUALLY MOVED BY -SHIFT.
752//
753        IF ((XDELTA>= 0) && (YDELTA >= 0))              // CASE I ---- "%&" is "and"
754                III[0] = II
755                JJJ[0] = JJ
756                III[1] = II + 1
757                JJJ[1] = JJ
758                III[2] = II + 1
759                JJJ[2] = JJ + 1
760                III[3] = II
761                JJJ[3] = JJ + 1
762                A[0] = (1. - XDELTA)*(1. - YDELTA)
763                A[1] = XDELTA*(1. - YDELTA)
764                A[2] = XDELTA*YDELTA
765                A[3] = (1. - XDELTA)*YDELTA
766        Endif
767        IF ((XDELTA >= 0) && (YDELTA < 0))              // CASE II.
768                III[0] = II
769                JJJ[0] = JJ
770                III[1] = II
771                JJJ[1] = JJ - 1
772                III[2] = II + 1
773                JJJ[2] = JJ - 1
774                III[3] = II + 1
775                JJJ[3] = JJ
776                A[0] = (1. - XDELTA)*(1. + YDELTA)
777                A[1] = (1. - XDELTA)*(-YDELTA)
778                A[2] = XDELTA*(-YDELTA)
779                A[3] = XDELTA*(1. + YDELTA)
780        Endif
781        IF ((XDELTA < 0) && (YDELTA >= 0))              // CASE III.
782                III[0] = II
783                JJJ[0] = JJ
784                III[1] = II
785                JJJ[1] = JJ + 1
786                III[2] = II - 1
787                JJJ[2] = JJ + 1
788                III[3] = II - 1
789                JJJ[3] = JJ
790                A[0] = (1. + XDELTA)*(1 - YDELTA)
791                A[1] = (1. + XDELTA)*YDELTA
792                A[2] = -XDELTA*YDELTA
793                A[3] = -XDELTA*(1. - YDELTA)
794        Endif
795        IF ((XDELTA < 0) && (YDELTA < 0))               //CASE IV.
796                III[0] = II
797                JJJ[0] = JJ
798                III[1] = II - 1
799                JJJ[1] = JJ
800                III[2] = II - 1
801                JJJ[2] = JJ - 1
802                III[3] = II
803                JJJ[3] = JJ - 1
804                A[0] = (1. + XDELTA)*(1. + YDELTA)
805                A[1] = -XDELTA*(1. + YDELTA)
806                A[2] = (-XDELTA)*(-YDELTA)
807                A[3] = (1. + XDELTA)*(-YDELTA)
808        Endif
809
810        NVAR pixelsX = root:myGlobals:gNPixelsX
811        NVAR pixelsY = root:myGlobals:gNPixelsY
812//check to see if iii[0],jjj[0] are valid detector elements, in [0,127]
813//if not set noadd[0] to 1, to let calling routine know NOT to add
814//        CALL TESTIJ(III(1),JJJ(1),OKIJ)
815        NOADD[0] = 0
816        if( (iii[0]<0) || (iii[0]>(pixelsX-1)) )
817                noadd[0] = 1
818        endif
819        if((jjj[0]<0) || (jjj[0]>(pixelsY-1)) )
820                noadd[0] = 1
821        endif
822       
823
824       
825        sum_val = 0.
826        kk = 0
827        Do
828                IF(JJJ[kk] == pixelsX)
829                        //do nothing
830                else
831                        sum_val += A[kk]*DATA[III[kk]][JJJ[kk]]
832                endif
833                kk+=1
834        while(kk<4)
835       
836        //clean up waves
837        KillWaves/z iii,jjj,a
838       
839        RETURN (sum_val)
840       
841End             //function ShiftSum
842
843//************************
844//unused testing procedure, may not be up-to-date with other procedures
845//check before re-implementing
846//
847Proc DIV_a_Workfile(type)
848        String type
849        Prompt type,"WORK data type",popup,"COR;SAM;EMP;BGD"
850       
851        //macro will take whatever is in SELECTED folder and DIVide it by the current
852        //contents of the DIV folder - the function will check for existence
853        //before proceeding
854       
855        Variable err
856        err = Divide_work(type)         //returns err = 1 if data doesn't exist in specified folders
857       
858        if(err)
859                Abort "error in Divide_work"
860        endif
861       
862        //contents are always dumped to CAL
863        type = "CAL"
864       
865        String newTitle = "WORK_"+type
866        DoWindow/F SANS_Data
867        DoWindow/T SANS_Data, newTitle
868        KillStrings/Z newTitle
869       
870        //need to update the display with "data" from the correct dataFolder
871        //reset the current displaytype to "type"
872        String/G root:myGlobals:gDataDisplayType=Type
873       
874        fRawWindowHook()
875       
876End
877
878//function will divide the contents of "type" folder with the contents of
879//the DIV folder
880// all data is converted to linear scale for the calculation
881//
882Function Divide_work(type)
883        String type
884       
885        //check for existence of data in type and DIV
886        // if the desired workfile doesn't exist, let the user know, and abort
887        String destPath=""
888
889        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
890                Print "There is no work file in "+type+"--Aborting"
891                Return(1)               //error condition
892        Endif
893        //check for DIV
894        // if the DIV workfile doesn't exist, let the user know,and abort
895
896        if(WaveExists($"root:Packages:NIST:DIV:data") == 0)
897                Print "There is no work file in DIV --Aborting"
898                Return(1)               //error condition
899        Endif
900        //files exist, proceed
901       
902        //check for log-scaling of the "DIV" data and adjust if necessary
903        ConvertFolderToLinearScale("DIV")
904       
905        //copy type information to CAL, wiping out the old contents of the CAL folder first
906       
907        //destPath = "root:Packages:NIST:CAL"
908        //SetDataFolder destPath
909        //KillWaves/A/Z                 //get rid of the old data in CAL folder
910
911        //check for log-scaling of the "type" data and adjust if necessary
912        ConvertFolderToLinearScale(type)
913        //then continue
914
915        //copy from current dir (type)=destPath to CAL, overwriting CAL contents
916        destPath = "root:Packages:NIST:" + type
917        Duplicate/O $(destPath + ":data"),$"root:Packages:NIST:CAL:data"
918        Duplicate/O $(destPath + ":linear_data"),$"root:Packages:NIST:CAL:linear_data"
919        Duplicate/O $(destPath + ":linear_data_error"),$"root:Packages:NIST:CAL:linear_data_error"
920//      Duplicate/O $(destPath + ":vlegend"),$"root:Packages:NIST:CAL:vlegend"
921        Duplicate/O $(destPath + ":textread"),$"root:Packages:NIST:CAL:textread"
922        Duplicate/O $(destPath + ":integersread"),$"root:Packages:NIST:CAL:integersread"
923        Duplicate/O $(destPath + ":realsread"),$"root:Packages:NIST:CAL:realsread"
924        //need to save a copy of filelist string too (from the current type folder)
925        SVAR oldFileList = $(destPath + ":fileList")
926
927        //now switch to reference waves in CAL folder
928        destPath = "root:Packages:NIST:CAL"
929        //make appropriate wave references
930        Wave data=$(destPath + ":linear_data")                                  // these wave references point to the data in CAL
931//      Wave data_err=$(destPath + ":linear_data_err")                                  // these wave references point to the data in CAL
932        Wave data_copy=$(destPath + ":data")                                    // these wave references point to the data in CAL
933        Wave/t textread=$(destPath + ":textread")                       //that are to be directly operated on
934        Wave integersread=$(destPath + ":integersread")
935        Wave realsread=$(destPath + ":realsread")
936        Variable/G $(destPath + ":gIsLogScale")=0                       //make new flag in CAL folder, data is linear scale
937        //need to copy filelist string too
938        String/G $(destPath + ":fileList") = oldFileList
939
940        Wave div_data = $"root:Packages:NIST:DIV:data"          //hard-wired in....
941        //do the division, changing data in CAL
942        data /= div_data
943       
944//      data_err /= div_data
945       
946        // keep "data" in sync with linear_data
947        data_copy = data
948       
949        //update CAL header
950        textread[1] = date() + " " + time()             //date + time stamp
951       
952        Return(0)
953End
954
955//test procedure, not called anymore
956Proc AbsoluteScaling(type,c0,c1,c2,c3,c4,c5)
957        String type
958        Variable c0=1,c1=0.1,c2=0.95,c3=0.1,c4=1,c5=32.0
959        Prompt type,"WORK data type",popup,"CAL;COR;SAM"
960        Prompt c0, "Sample Transmission"
961        Prompt c1, "Sample Thickness (cm)"
962        Prompt c2, "Standard Transmission"
963        Prompt c3, "Standard Thickness (cm)"
964        Prompt c4, "I(0) from standard fit (normalized to 1E8 monitor cts)"
965        Prompt c5, "Standard Cross-Section (cm-1)"
966
967        Variable err
968        //call the function to do the math
969        //data from "type" will be scaled and deposited in ABS
970        err = Absolute_Scale(type,c0,c1,c2,c3,c4,c5)
971       
972        if(err)
973                Abort "Error in Absolute_Scale()"
974        endif
975       
976        //contents are always dumped to ABS
977        type = "ABS"
978       
979        String newTitle = "WORK_"+type
980        DoWindow/F SANS_Data
981        DoWindow/T SANS_Data, newTitle
982        KillStrings/Z newTitle
983       
984        //need to update the display with "data" from the correct dataFolder
985        //reset the current displaytype to "type"
986        String/G root:myGlobals:gDataDisplayType=Type
987       
988        fRawWindowHook()
989       
990End
991
992//s_ is the standard
993//w_ is the "work" file
994//both are work files and should already be normalized to 10^8 monitor counts
995Function Absolute_Scale(type,w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err)
996        String type
997        Variable w_trans,w_thick,s_trans,s_thick,s_izero,s_cross,kappa_err
998       
999        //convert the "type" data to absolute scale using the given standard information
1000        //copying the "type" waves to ABS
1001       
1002        //check for existence of data, rescale to linear if needed
1003        String destPath
1004        //check for "type"
1005        if(WaveExists($("root:Packages:NIST:"+Type + ":data")) == 0)
1006                Print "There is no work file in "+type+"--Aborting"
1007                Return(1)               //error condition
1008        Endif
1009        //check for log-scaling of the "type" data and adjust if necessary
1010        destPath = "root:Packages:NIST:"+Type
1011        NVAR gIsLogScale = $(destPath + ":gIsLogScale")
1012        if(gIsLogScale)
1013                Duplicate/O $(destPath + ":linear_data") $(destPath + ":data")//back to linear scale
1014                Variable/G $(destPath + ":gIsLogScale")=0       //the "type" data is not logscale anymore
1015        endif
1016       
1017        //copy "oldtype" information to ABS
1018        //overwriting out the old contents of the ABS folder (/O option in Duplicate)
1019        //copy over the waves data,vlegend,text,integers,reals(read)
1020
1021        String oldType= "root:Packages:NIST:"+type              //this is where the data to be absoluted is
1022        //copy from current dir (type) to ABS, defined by destPath
1023        Duplicate/O $(oldType + ":data"),$"root:Packages:NIST:ABS:data"
1024        Duplicate/O $(oldType + ":linear_data"),$"root:Packages:NIST:ABS:linear_data"
1025        Duplicate/O $(oldType + ":linear_data_error"),$"root:Packages:NIST:ABS:linear_data_error"
1026//      Duplicate/O $(oldType + ":vlegend"),$"root:Packages:NIST:ABS:vlegend"
1027        Duplicate/O $(oldType + ":textread"),$"root:Packages:NIST:ABS:textread"
1028        Duplicate/O $(oldType + ":integersread"),$"root:Packages:NIST:ABS:integersread"
1029        Duplicate/O $(oldType + ":realsread"),$"root:Packages:NIST:ABS:realsread"
1030        //need to save a copy of filelist string too (from the current type folder)
1031        SVAR oldFileList = $(oldType + ":fileList")
1032        //need to copy filelist string too
1033        String/G $"root:Packages:NIST:ABS:fileList" = oldFileList
1034       
1035        //now switch to ABS folder
1036        //make appropriate wave references
1037        WAVE data=$"root:Packages:NIST:ABS:linear_data"                                 // these wave references point to the "type" data in ABS
1038        WAVE data_err=$"root:Packages:NIST:ABS:linear_data_error"                                       // these wave references point to the "type" data in ABS
1039        WAVE data_copy=$"root:Packages:NIST:ABS:data"                                   // just for display
1040        WAVE/T textread=$"root:Packages:NIST:ABS:textread"                      //that are to be directly operated on
1041        WAVE integersread=$"root:Packages:NIST:ABS:integersread"
1042        WAVE realsread=$"root:Packages:NIST:ABS:realsread"
1043        Variable/G $"root:Packages:NIST:ABS:gIsLogscale"=0                      //make new flag in ABS folder, data is linear scale
1044       
1045        //do the actual absolute scaling here, modifying the data in ABS
1046        Variable defmon = 1e8,w_moncount,s1,s2,s3,s4
1047       
1048        w_moncount = realsread[0]               //monitor count in "type"
1049        if(w_moncount == 0)
1050                //zero monitor counts will give divide by zero ---
1051                DoAlert 0,"Total monitor count in data file is zero. No rescaling of data"
1052                Return(1)               //report error
1053        Endif
1054       
1055        //calculate scale factor
1056        Variable scale,trans_err
1057        s1 = defmon/realsread[0]                //[0] is monitor count (s1 should be 1)
1058        s2 = s_thick/w_thick
1059        s3 = s_trans/w_trans
1060        s4 = s_cross/s_izero
1061       
1062        // kappa comes in as s_izero, so be sure to use 1/kappa_err
1063       
1064        data *= s1*s2*s3*s4
1065       
1066        scale = s1*s2*s3*s4
1067        trans_err = realsRead[41]
1068       
1069//      print scale
1070//      print data[0][0]
1071       
1072        data_err = sqrt(scale^2*data_err^2 + scale^2*data^2*(kappa_err^2/s_izero^2 +trans_err^2/w_trans^2))
1073
1074//      print data_err[0][0]
1075       
1076// keep "data" in sync with linear_data
1077        data_copy = data
1078       
1079        //********* 15APR02
1080        // DO NOt correct for atenuators here - the COR step already does this, putting all of the data one equal
1081        // footing (zero atten) before doing the subtraction.
1082        //
1083        //Print "ABS data multiplied by  ",s1*s2*s3*s4/attenFactor
1084       
1085        //update the ABS header information
1086        textread[1] = date() + " " + time()             //date + time stamp
1087       
1088        Return (0) //no error
1089End
1090
1091//*************
1092// start of section of functions used for workfile math panel
1093//*************
1094
1095
1096//function will copy the contents of oldtype folder to newtype folder
1097//converted to linear scale before copying
1098//******data in newtype is overwritten********
1099//
1100Function CopyWorkContents(oldtype,newtype)
1101        String oldtype,newtype
1102       
1103        //check for existence of data in oldtype
1104        // if the desired workfile doesn't exist, let the user know, and abort
1105        String destPath=""
1106        if(WaveExists($("root:Packages:NIST:"+oldType + ":data")) == 0)
1107                Print "There is no work file in "+oldtype+"--Aborting"
1108                Return(1)               //error condition
1109        Endif
1110       
1111        //check for log-scaling of the "type" data and adjust if necessary
1112        ConvertFolderToLinearScale(oldtype)
1113        Fix_LogLinButtonState(0)                //make sure the button reflects the new linear scaling
1114        //then continue
1115
1116        //copy from current dir (type)=destPath to newtype, overwriting newtype contents
1117        destPath = "root:Packages:NIST:" + oldtype
1118        Duplicate/O $(destPath + ":data"),$("root:Packages:NIST:"+newtype+":data")
1119        Duplicate/O $(destPath + ":linear_data"),$("root:Packages:NIST:"+newtype+":linear_data")
1120        Duplicate/O $(destPath + ":linear_data_error"),$("root:Packages:NIST:"+newtype+":linear_data_error")
1121        Duplicate/O $(destPath + ":textread"),$("root:Packages:NIST:"+newtype+":textread")
1122        Duplicate/O $(destPath + ":integersread"),$("root:Packages:NIST:"+newtype+":integersread")
1123        Duplicate/O $(destPath + ":realsread"),$("root:Packages:NIST:"+newtype+":realsread")
1124        //
1125        // be sure to get rid of the linear_data if it exists in the destination folder
1126//      KillWaves/Z $("root:Packages:NIST:"+newtype+":linear_data")
1127
1128        //need to save a copy of filelist string too (from the current type folder)
1129        SVAR oldFileList = $(destPath + ":fileList")
1130
1131        //now switch to reference waves in newtype folder
1132        destPath = "root:Packages:NIST:"+newtype
1133        Variable/G $(destPath + ":gIsLogScale")=0                       //make new flag in newtype folder, data is linear scale
1134        //need to copy filelist string too
1135        String/G $(destPath + ":fileList") = oldFileList
1136
1137        Return(0)
1138End
1139
1140//Entry procedure from main panel
1141//
1142Proc CopyWorkFolder(oldType,newType)
1143        String oldType,newType
1144        Prompt oldType,"Source WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
1145        Prompt newType,"Destination WORK data type",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;STO;SUB;DRK;"
1146//      Prompt oldType,"Source WORK data type",popup,"AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1147//      Prompt newType,"Destination WORK data type",popup,"AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1148
1149        // data folder "old" will be copied to "new" (and will overwrite)
1150        CopyWorkContents(oldtype,newtype)
1151End
1152
1153//initializes datafolder and constants necessary for the workfilemath panel
1154//
1155Proc Init_WorkMath()
1156        //create the data folder
1157        //String str = "AAA;BBB;CCC;DDD;EEE;FFF;GGG;"
1158        String str = "File_1;File_2;Result;"
1159        NewDataFolder/O/S root:Packages:NIST:WorkMath
1160        String/G gFolderList=str
1161        Variable ii=0,num=itemsinlist(str)
1162        do
1163                Execute "NewDataFolder/O "+StringFromList(ii, str ,";")
1164                ii+=1
1165        while(ii<num)
1166        Variable/G const1=1,const2=1
1167       
1168        SetDataFolder root:
1169End
1170
1171//entry procedure to invoke the workfilemath panel, initializes if needed
1172//
1173Proc Show_WorkMath_Panel()
1174        DoWindow/F WorkFileMath
1175        if(V_flag==0)
1176                Init_WorkMath()
1177                WorkFileMath()
1178        Endif
1179End
1180
1181//attempts to perform the selected math operation based on the selections in the panel
1182// aborts with an error condition if operation can't be completed
1183// or puts the final answer in the Result folder, and displays the selected data
1184//
1185Function WorkMath_DoIt_ButtonProc(ctrlName) : ButtonControl
1186        String ctrlName
1187
1188        String str1,str2,oper,dest = "Result"
1189        String pathStr,workMathStr="WorkMath:"
1190       
1191        //get the panel selections (these are the names of the files on disk)
1192        PathInfo catPathName
1193        pathStr=S_Path
1194        ControlInfo popup0
1195        str1=S_Value
1196        ControlInfo popup1
1197        str2=S_Value
1198        ControlInfo popup2
1199        oper=S_Value
1200       
1201        //check that something has been selected for operation and destination
1202        if(cmpstr(oper,"Operation")==0)
1203                Abort "Select a math operand from the popup"
1204        Endif
1205
1206        //constants from globals
1207        NVAR const1=root:Packages:NIST:WorkMath:const1
1208        NVAR const2=root:Packages:NIST:WorkMath:const2
1209        Printf "(%g)%s %s (%g)%s = %s\r", const1,str1,oper,const2,str2,dest
1210        //check for proper folders (all 3 must be different)
1211       
1212        //load the data in here...
1213        //set #1
1214        Load_NamedASC_File(pathStr+str1,workMathStr+"File_1")
1215       
1216        NVAR pixelsX = root:myGlobals:gNPixelsX         //OK, location is correct
1217        NVAR pixelsY = root:myGlobals:gNPixelsY
1218       
1219        WAVE/Z data1=$("root:Packages:NIST:"+workMathStr+"File_1:data")
1220        WAVE/Z err1=$("root:Packages:NIST:"+workMathStr+"File_1:linear_data_error")
1221       
1222        // set # 2
1223        If(cmpstr(str2,"UNIT MATRIX")==0)
1224                Make/O/N=(pixelsX,pixelsY) root:Packages:NIST:WorkMath:data             //don't put in File_2 folder
1225                Wave/Z data2 =  root:Packages:NIST:WorkMath:data                        //it's not real data!
1226                data2=1
1227                Duplicate/O data2 err2
1228                err2 = 0
1229        else
1230                //Load set #2
1231                Load_NamedASC_File(pathStr+str2,workMathStr+"File_2")
1232                WAVE/Z data2=$("root:Packages:NIST:"+workMathStr+"File_2:data")
1233                WAVE/Z err2=$("root:Packages:NIST:"+workMathStr+"File_2:linear_data_error")
1234        Endif
1235
1236        ///////
1237       
1238        //now that we know that data exists, convert each of the operands to linear scale
1239        ConvertFolderToLinearScale(workMathStr+"File_1")
1240        If(cmpstr(str2,"UNIT MATRIX")!=0)
1241                ConvertFolderToLinearScale(workMathStr+"File_2")                //don't need to convert unit matrix to linear
1242        endif
1243        //copy contents of str1 folder to dest and create the wave ref (it will exist)
1244        CopyWorkContents(workMathStr+"File_1",workMathStr+dest)
1245        WAVE/Z destData=$("root:Packages:NIST:"+workMathStr+dest+":data")
1246        WAVE/Z destErr=$("root:Packages:NIST:"+workMathStr+dest+":linear_data_error")
1247       
1248        //dispatch
1249        strswitch(oper)
1250                case "*":               //multiplication
1251                        destData = const1*data1 * const2*data2
1252                        destErr = const1^2*const2^2*(err1^2*data2^2 + err2^2*data1^2)
1253                        destErr = sqrt(destErr)
1254                        break   
1255                case "_":               //subtraction
1256                        destData = const1*data1 - const2*data2
1257                        destErr = const1^2*err1^2 + const2^2*err2^2
1258                        destErr = sqrt(destErr)
1259                        break
1260                case "/":               //division
1261                        destData = (const1*data1) / (const2*data2)
1262                        destErr = const1^2/const2^2*(err1^2/data2^2 + err2^2*data1^2/data2^4)
1263                        destErr = sqrt(destErr)
1264                        break
1265                case "+":               //addition
1266                        destData = const1*data1 + const2*data2
1267                        destErr = const1^2*err1^2 + const2^2*err2^2
1268                        destErr = sqrt(destErr)
1269                        break                   
1270        endswitch
1271       
1272        //show the result
1273        WorkMath_Display_PopMenuProc("",0,"Result")
1274End
1275
1276// closes the panel and kills the data folder when done
1277//
1278Function WorkMath_Done_ButtonProc(ctrlName) : ButtonControl
1279        String ctrlName
1280       
1281        DoAlert 1,"Closing the panel will kill all of the data you have loaded into memory. Do you want to continue?"
1282        If(V_Flag==2)
1283                return(0)               //don't do anything
1284        Endif
1285        //kill the panel
1286        DoWindow/K WorkFileMath
1287        //wipe out the data folder of globals
1288        SVAR dataType=root:myGlobals:gDataDisplayType
1289        if(strsearch(dataType, "WorkMath", 0 ) != -1)           //kill the SANS_Data graph if needed
1290                DoWindow/K SANS_Data
1291        Endif
1292        KillDataFolder root:Packages:NIST:WorkMath
1293End
1294
1295// loads data into the specified folder
1296//
1297// currently unused since button has been removed from panel
1298//
1299Function WorkMath_Load_ButtonProc(ctrlName) : ButtonControl
1300        String ctrlName
1301       
1302        String destStr=""
1303        SVAR folderList=root:Packages:NIST:WorkMath:gFolderList
1304        Prompt destStr,"Select the destination folder",popup,folderList
1305        DoPrompt "Folder for ASC Load",destStr
1306       
1307        if(V_flag==1)
1308                return(1)               //user abort, do nothing
1309        Endif
1310       
1311        String destFolder = "WorkMath:"+destStr
1312       
1313        Load_ASC_File("Pick the ASC file",destFolder)
1314End
1315
1316// changes the display of the SANS_Data window based on popped data type
1317// first loads in the data from the File1 or File2 popup as needed
1318// then displays the selcted dataset, if it exists
1319// makes use of procedure from DisplayUtils
1320//
1321// - Always replaces File1 or File2 with fresh data from disk
1322//
1323Function WorkMath_Display_PopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
1324        String ctrlName
1325        Variable popNum
1326        String popStr
1327       
1328        String folder="WorkMath:",pathStr,str1
1329
1330        PathInfo catPathName
1331        pathStr=S_Path
1332       
1333        //if display result, just do it and return
1334        if(cmpstr(popStr,"Result")==0)
1335                Execute "ChangeDisplay(\""+folder+popstr+"\")"
1336                return(0)
1337        endif
1338        // if file1 or file2, load in the data and display
1339        if(cmpstr(popStr,"File_1")==0)
1340                ControlInfo popup0
1341                str1 = S_Value
1342        Endif
1343        if(cmpstr(popStr,"File_2")==0)
1344                ControlInfo popup1
1345                str1 = S_Value
1346        Endif
1347        //don't load or display the unit matrix
1348        Print str1
1349        if(cmpstr(str1,"UNIT MATRIX")!=0)
1350                Load_NamedASC_File(pathStr+str1,folder+popStr)
1351                //change the display
1352                Execute "ChangeDisplay(\""+folder+popstr+"\")"
1353        endif
1354        return(0)       
1355End
1356
1357//simple panel to do workfile arithmetic
1358//
1359Proc WorkFileMath()
1360        PauseUpdate; Silent 1           // building window...
1361        NewPanel /W=(610,211,880,490)/K=2 as "Work File Math"           // replace /K=2 flag
1362        DoWindow/C WorkFileMath
1363        ModifyPanel cbRGB=(47802,54484,6682)
1364        ModifyPanel fixedSize=1
1365        SetDrawLayer UserBack
1366        DrawLine 6,166,214,166
1367        SetDrawEnv fstyle= 4
1368        DrawText 10,34,"File #1:"
1369        SetDrawEnv fstyle= 4
1370        DrawText 13,129,"File #2:"
1371        DrawText 78,186,"= Result"
1372        Button button0,pos={28,245},size={50,20},proc=WorkMath_DoIt_ButtonProc,title="Do It"
1373        Button button0,help={"Performs the specified arithmetic"}
1374        Button button1,pos={183,245},size={50,20},proc=WorkMath_Done_ButtonProc,title="Done"
1375        Button button1,help={"Closes the panel"}
1376//      Button button2,pos={30,8},size={110,20},proc=WorkMath_Load_ButtonProc,title="Load ASC Data"
1377//      Button button2,help={"Loads ASC data files into the specified folder"}
1378        Button button3,pos={205,8},size={25,20},proc=ShowWorkMathHelp,title="?"
1379        Button button3,help={"Show help file for math operations on 2-D data sets"}
1380        SetVariable setvar0,pos={9,46},size={70,15},title=" "
1381        SetVariable setvar0,limits={-Inf,Inf,0},value= root:Packages:NIST:WorkMath:const1
1382        SetVariable setvar0,help={"Multiplicative constant for the first dataset"}
1383        PopupMenu popup0,pos={89,44},size={84,20},title="X  "
1384        PopupMenu popup0,mode=1,popvalue="1st Set",value= ASC_FileList()
1385        PopupMenu popup0,help={"Selects the first dataset for the operation"}
1386        PopupMenu popup1,pos={93,136},size={89,20},title="X  "
1387        PopupMenu popup1,mode=1,popvalue="2nd Set",value= "UNIT MATRIX;"+ASC_FileList()
1388        PopupMenu popup1,help={"Selects the second dataset for the operation"}
1389        PopupMenu popup2,pos={50,82},size={70,20},title="Operator  "
1390        PopupMenu popup2,mode=3,popvalue="Operation",value= #"\"+;_;*;/;\""
1391        PopupMenu popup2,help={"Selects the mathematical operator"}
1392        SetVariable setvar1,pos={13,139},size={70,15},title=" "
1393        SetVariable setvar1,limits={-Inf,Inf,0},value= root:Packages:NIST:WorkMath:const2
1394        SetVariable setvar1,help={"Multiplicative constant for the second dataset"}
1395//      PopupMenu popup3,pos={27,167},size={124,20},title=" = Destination"
1396//      PopupMenu popup3,mode=1,popvalue="Destination",value= root:Packages:NIST:WorkMath:gFolderList
1397//      PopupMenu popup3,help={"Selects the destination folder"}
1398        PopupMenu popup4,pos={55,204},size={103,20},proc=WorkMath_Display_PopMenuProc,title="Display"
1399        PopupMenu popup4,mode=3,value= "File_1;File_2;Result;"
1400        PopupMenu popup4,help={"Displays the data in the specified folder"}
1401EndMacro
1402
1403//jump to the help topic
1404Function ShowWorkMathHelp(ctrlName) : ButtonControl
1405        String ctrlName
1406        DisplayHelpTopic/Z/K=1 "SANS Data Reduction Tutorial[2-D Work File Arithmetic]"
1407        if(V_flag !=0)
1408                DoAlert 0,"The SANS Data Reduction Tutorial Help file could not be found"
1409        endif
1410End
1411
1412//utility function to clear the contents of a data folder
1413//won't clear data that is in use -
1414//
1415Function ClearDataFolder(type)
1416        String type
1417       
1418        SetDataFolder $("root:Packages:NIST:"+type)
1419        KillWaves/a/z
1420        KillStrings/a/z
1421        KillVariables/a/z
1422        SetDataFolder root:
1423End
1424
1425
1426
1427//fileStr must be the FULL path and filename on disk
1428//destFolder is the path to the Igor folder where the data is to be deposited
1429// - "Packages:NIST:WorkMath:File_1" for example, compatible with SANS_Data display type
1430//
1431Function Load_NamedASC_File(fileStr,destFolder)
1432        String fileStr,destFolder
1433
1434        Variable refnum
1435       
1436        //read in the data
1437        ReadASCData(fileStr,destFolder)
1438
1439        //the calling macro must change the display type
1440        String/G root:myGlobals:gDataDisplayType=destFolder
1441       
1442        FillFakeHeader_ASC(destFolder)          //uses info on the panel, if available
1443
1444        //data is displayed here, and needs header info
1445       
1446        fRawWindowHook()
1447       
1448        return(0)
1449End
1450
1451
1452//function called by the main entry procedure (the load button)
1453//loads data into the specified folder, and wipes out what was there
1454//
1455//aborts if no file was selected
1456//
1457// reads the data in if all is OK
1458//
1459// currently unused, as load button has been replaced
1460//
1461Function Load_ASC_File(msgStr,destFolder)
1462        String msgStr,destFolder
1463
1464        String filename="",pathStr=""
1465        Variable refnum
1466
1467        //check for the path
1468        PathInfo catPathName
1469        If(V_Flag==1)           //      /D does not open the file
1470                Open/D/R/T="????"/M=(msgStr)/P=catPathName refNum
1471        else
1472                Open/D/R/T="????"/M=(msgStr) refNum
1473        endif
1474        filename = S_FileName           //get the filename, or null if canceled from dialog
1475        if(strlen(filename)==0)
1476                //user cancelled, abort
1477                SetDataFolder root:
1478                Abort "No file selected, action aborted"
1479        Endif
1480       
1481        //read in the data
1482        ReadASCData(filename,destFolder)
1483
1484        //the calling macro must change the display type
1485        String/G root:myGlobals:gDataDisplayType=destFolder
1486       
1487        FillFakeHeader_ASC(destFolder)          //uses info on the panel, if available
1488
1489        //data is displayed here, and needs header info
1490       
1491        fRawWindowHook()
1492       
1493        return(0)
1494End
1495
1496
1497
1498//function called by the popups to get a file list of data that can be sorted
1499// this procedure simply removes the raw data files from the string - there
1500//can be lots of other junk present, but this is very fast...
1501//
1502// -- simplified to do in a single call -- extension must be .ASC
1503// - so far, only used in WorkFileMath popups, which require .ASC
1504//
1505Function/S ASC_FileList()
1506
1507        String list="",newList="",item=""
1508        Variable num,ii
1509       
1510        //check for the path
1511        PathInfo catPathName
1512        if(V_Flag==0)
1513                DoAlert 0, "Data path does not exist - pick the data path from the button on the main panel"
1514                Return("")
1515        Endif
1516
1517        list = IndexedFile(catpathName,-1,".ASC")
1518
1519
1520//      list = IndexedFile(catpathName,-1,"????")
1521//     
1522//      list = RemoveFromList(ListMatch(list,"*.SA1*",";"), list, ";", 0)
1523//      list = RemoveFromList(ListMatch(list,"*.SA2*",";"), list, ";", 0)
1524//      list = RemoveFromList(ListMatch(list,"*.SA3*",";"), list, ";", 0)
1525//      list = RemoveFromList(ListMatch(list,".*",";"), list, ";", 0)
1526//      list = RemoveFromList(ListMatch(list,"*.pxp",";"), list, ";", 0)
1527//      list = RemoveFromList(ListMatch(list,"*.DIV",";"), list, ";", 0)
1528//      list = RemoveFromList(ListMatch(list,"*.GSP",";"), list, ";", 0)
1529//      list = RemoveFromList(ListMatch(list,"*.MASK",";"), list, ";", 0)
1530//
1531//      //remove VAX version numbers
1532//      list = RemoveVersNumsFromList(List)
1533        //sort
1534        newList = SortList(List,";",0)
1535
1536        return newlist
1537End
Note: See TracBrowser for help on using the repository browser.