source: sans/SANSReduction/branches/kline_29MAR07/Put in User Procedures/SANS_Reduction_v5.00/WriteQIS.ipf @ 66

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

resolution calculation now uses the projected area of the beamstop for the shadowing factor and all other resolution calculations. per JGB 3/21/07 ticket 751
File size: 24.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=4.0
4
5//************************
6// Vers 1.2 091001
7//
8//************************
9
10
11//***this function is not used - WriteWaves_W_Protocol() is used instead
12//
13//for writing out data (q-i-s) from the "type" folder
14//if fullpath is a complete HD path:filename, no dialog will be presented
15//if fullpath is just a filename, the save dialog will be presented (forced if dialog =1)
16//writes ONLY the standard header information (no protocol information)
17//
18Function WriteWaves(type,fullpath,dialog)
19        String type,fullpath
20        Variable dialog         //=1 will present dialog for name
21       
22        String destStr=""
23        destStr = "root:"+type
24       
25        Variable refNum
26        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
27        String fname,ave="C",headerFormat = "%10.4g %8.2g %8.2g %8.2g %8.3g %8.3g %8s %5.0g\r\n"
28        Variable step=1
29       
30        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
31        WAVE intw=$(destStr + ":integersRead")
32        WAVE rw=$(destStr + ":realsRead")
33        WAVE/T textw=$(destStr + ":textRead")
34        WAVE qvals =$(destStr + ":qval")
35        WAVE inten=$(destStr + ":aveint")
36        WAVE sig=$(destStr + ":sigave")
37        WAVE qbar = $(destStr + ":QBar")
38        WAVE sigmaq = $(destStr + ":SigmaQ")
39        WAVE fsubs = $(destStr + ":fSubS")
40
41        //check each wave
42        If(!(WaveExists(intw)))
43                Abort "intw DNExist BinaryWrite()"
44        Endif
45        If(!(WaveExists(rw)))
46                Abort "rw DNExist BinaryWrite()"
47        Endif
48        If(!(WaveExists(textw)))
49                Abort "textw DNExist BinaryWrite()"
50        Endif
51        If(!(WaveExists(qvals)))
52                Abort "qvals DNExist BinaryWrite()"
53        Endif
54        If(!(WaveExists(inten)))
55                Abort "inten DNExist BinaryWrite()"
56        Endif
57        If(!(WaveExists(sig)))
58                Abort "sig DNExist BinaryWrite()"
59        Endif
60        If(!(WaveExists(qbar)))
61                Abort "qbar DNExist BinaryWrite()"
62        Endif
63        If(!(WaveExists(sigmaq)))
64                Abort "sigmaq DNExist BinaryWrite()"
65        Endif
66        If(!(WaveExists(fsubs)))
67                Abort "fsubs DNExist BinaryWrite()"
68        Endif
69               
70        if(dialog)
71                PathInfo/S catPathName
72                fullPath = DoSaveFileDialog("Save data as")
73                If(cmpstr(fullPath,"")==0)
74                        //user cancel, don't write out a file
75                        Close/A
76                        Abort "no data file was written"
77                Endif
78                //Print "dialog fullpath = ",fullpath
79        Endif
80       
81        //actually open the file
82        Open refNum as fullpath
83       
84        fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n",textw[0],textw[1]
85        fprintf refnum,"LABEL: %s\r\n",textw[6]
86        fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
87        fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
88        wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
89       
90        Close refnum
91       
92        SetDataFolder root:             //(redundant)
93       
94        Return(0)
95End
96
97
98//for writing out data (q-i-s) from the "type" folder, and including reduction information
99//if fullpath is a complete HD path:filename, no dialog will be presented
100//if fullpath is just a filename, the save dialog will be presented
101//if dialog = 1, a dialog will always be presented
102//
103// root:myGlobals:Protocols:gProtoStr is the name of the currently active protocol
104//
105Function WriteWaves_W_Protocol(type,fullpath,dialog)
106        String type,fullpath
107        Variable dialog         //=1 will present dialog for name
108       
109        String destStr=""
110        destStr = "root:"+type
111       
112        Variable refNum
113        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
114        String fname,ave="C",hdrStr1="",hdrStr2=""
115        Variable step=1
116       
117       
118       
119        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
120        WAVE intw=$(destStr + ":integersRead")
121        WAVE rw=$(destStr + ":realsRead")
122        WAVE/T textw=$(destStr + ":textRead")
123        WAVE qvals =$(destStr + ":qval")
124        WAVE inten=$(destStr + ":aveint")
125        WAVE sig=$(destStr + ":sigave")
126        WAVE qbar = $(destStr + ":QBar")
127        WAVE sigmaq = $(destStr + ":SigmaQ")
128        WAVE fsubs = $(destStr + ":fSubS")
129
130        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
131        Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
132       
133        //check each wave
134        If(!(WaveExists(intw)))
135                Abort "intw DNExist BinaryWrite_W_Protocol()"
136        Endif
137        If(!(WaveExists(rw)))
138                Abort "rw DNExist BinaryWrite_W_Protocol()"
139        Endif
140        If(!(WaveExists(textw)))
141                Abort "textw DNExist BinaryWrite_W_Protocol()"
142        Endif
143        If(!(WaveExists(qvals)))
144                Abort "qvals DNExist BinaryWrite_W_Protocol()"
145        Endif
146        If(!(WaveExists(inten)))
147                Abort "inten DNExist BinaryWrite_W_Protocol()"
148        Endif
149        If(!(WaveExists(sig)))
150                Abort "sig DNExist BinaryWrite_W_Protocol()"
151        Endif
152        If(!(WaveExists(qbar)))
153                Abort "qbar DNExist BinaryWrite_W_Protocol()"
154        Endif
155        If(!(WaveExists(sigmaq)))
156                Abort "sigmaq DNExist BinaryWrite_W_Protocol()"
157        Endif
158        If(!(WaveExists(fsubs)))
159                Abort "fsubs DNExist BinaryWrite_W_Protocol()"
160        Endif
161        If(!(WaveExists(proto)))
162                Abort "current protocol wave DNExist BinaryWrite_W_Protocol()"
163        Endif
164
165        //strings can be too long to print-- must trim to 255 chars
166        Variable ii,num=8
167        Make/O/T/N=(num) tempShortProto
168        for(ii=0;ii<num;ii+=1)
169                tempShortProto[ii] = (proto[ii])[0,240]
170        endfor
171       
172        if(dialog)
173                PathInfo/S catPathName
174                fullPath = DoSaveFileDialog("Save data as")
175                If(cmpstr(fullPath,"")==0)
176                        //user cancel, don't write out a file
177                        Close/A
178                        Abort "no data file was written"
179                Endif
180                //Print "dialog fullpath = ",fullpath
181        Endif
182       
183        hdrStr1 = num2str(rw[0])+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
184        hdrStr1 += "     "+num2str(rw[4])+"     "+num2str(rw[5]) + ave +"   "+num2str(step) + "\r\n"
185
186        hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
187        hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+textW[9] + "\r\n"
188       
189        SVAR samFiles = $("root:"+type+":fileList")
190        //actually open the file here
191        Open refNum as fullpath
192       
193        //write out the standard header information
194        fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n",textw[0],textw[1]
195        fprintf refnum,"LABEL: %s\r\n",textw[6]
196        fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
197        fprintf refnum,hdrStr1
198        fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n"
199        fprintf refnum,hdrStr2
200//      fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
201
202        //insert protocol information here
203        //-1 list of sample files
204        //0 - bkg
205        //1 - emp
206        //2 - div
207        //3 - mask
208        //4 - abs params c2-c5
209        //5 - average params
210        fprintf refnum, "SAM: %s\r\n",samFiles
211        fprintf refnum, "BGD: %s\r\n",tempShortProto[0]
212        fprintf refnum, "EMP: %s\r\n",tempShortProto[1]
213        fprintf refnum, "DIV: %s\r\n",tempShortProto[2]
214        fprintf refnum, "MASK: %s\r\n",tempShortProto[3]
215        fprintf refnum, "ABS Parameters (3-6): %s\r\n",tempShortProto[4]
216        fprintf refnum, "Average Choices: %s\r\n",tempShortProto[5]
217       
218        //write out the data columns
219        fprintf refnum,"The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n"
220        wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
221       
222        Close refnum
223       
224        SetDataFolder root:             //(redundant)
225       
226        //write confirmation of write operation to history area
227        Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
228        KillWaves/Z tempShortProto
229        Return(0)
230End
231
232
233//for writing out data (phi-i-s) from the "type" folder, and including reduction information
234//if fullpath is a complete HD path:filename, no dialog will be presented
235//if fullpath is just a filename, the save dialog will be presented
236//if dialog = 1, a dialog will always be presented
237//
238// root:myGlobals:Protocols:gProtoStr is the name of the currently active protocol
239//
240Function WritePhiave_W_Protocol(type,fullpath,dialog)
241        String type,fullpath
242        Variable dialog         //=1 will present dialog for name
243       
244        String destStr
245        destStr = "root:"+type
246       
247        Variable refNum
248        String formatStr = "%15.4g %15.4g %15.4g\r\n"
249        String fname,ave="C",hdrStr1,hdrStr2
250        Variable step=1
251       
252        //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
253        WAVE intw=$(destStr + ":integersRead")
254        WAVE rw=$(destStr + ":realsRead")
255        WAVE/T textw=$(destStr + ":textRead")
256        WAVE phival =$(destStr + ":phival")
257        WAVE inten=$(destStr + ":aveint")
258        WAVE sig=$(destStr + ":sigave")
259        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
260        Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
261       
262        //check each wave
263        If(!(WaveExists(intw)))
264                Abort "intw DNExist BinaryWrite_W_Protocol()"
265        Endif
266        If(!(WaveExists(rw)))
267                Abort "rw DNExist BinaryWrite_W_Protocol()"
268        Endif
269        If(!(WaveExists(textw)))
270                Abort "textw DNExist BinaryWrite_W_Protocol()"
271        Endif
272        If(!(WaveExists(phival)))
273                Abort "qvals DNExist BinaryWrite_W_Protocol()"
274        Endif
275        If(!(WaveExists(inten)))
276                Abort "inten DNExist BinaryWrite_W_Protocol()"
277        Endif
278        If(!(WaveExists(sig)))
279                Abort "sig DNExist BinaryWrite_W_Protocol()"
280        Endif
281        If(!(WaveExists(proto)))
282                Abort "current protocol wave DNExist BinaryWrite_W_Protocol()"
283        Endif
284        //strings can be too long to print-- must trim to 255 chars
285        Variable ii,num=8
286        Make/O/T/N=(num) tempShortProto
287        for(ii=0;ii<num;ii+=1)
288                tempShortProto[ii] = (proto[ii])[0,240]
289        endfor
290       
291        if(dialog)
292                PathInfo/S catPathName
293                fullPath = DoSaveFileDialog("Save data as")
294                If(cmpstr(fullPath,"")==0)
295                        //user cancel, don't write out a file
296                        Close/A
297                        Abort "no data file was written"
298                Endif
299                //Print "dialog fullpath = ",fullpath
300        Endif
301       
302        hdrStr1 = num2str(rw[0])+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
303        hdrStr1 += "     "+num2str(rw[4])+"     "+num2str(rw[5]) + ave +"   "+num2str(step) + "\r\n"
304
305        hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
306        hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+textW[9] + "\r\n"
307       
308        SVAR samFiles = $("root:"+type+":fileList")
309        //actually open the file here
310        Open refNum as fullpath
311       
312        //write out the standard header information
313        fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n",textw[0],textw[1]
314        fprintf refnum,"LABEL: %s\r\n",textw[6]
315        fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
316        fprintf refnum,hdrStr1
317        fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n"
318        fprintf refnum,hdrStr2
319       
320        //insert protocol information here
321        //0 - bkg
322        //1 - emp
323        //2 - div
324        //3 - mask
325        //4 - abs params c2-c5
326        //5 - average params
327        fprintf refnum, "SAM: %s\r\n",samFiles
328        fprintf refnum, "BGD: %s\r\n",tempShortProto[0]
329        fprintf refnum, "EMP: %s\r\n",tempShortProto[1]
330        fprintf refnum, "DIV: %s\r\n",tempShortProto[2]
331        fprintf refnum, "MASK: %s\r\n",tempShortProto[3]
332        fprintf refnum, "ABS Parameters (3-6): %s\r\n",tempShortProto[4]
333        fprintf refnum, "Average Choices: %s\r\n",tempShortProto[5]
334       
335        //write out the data columns
336        fprintf refnum,"The 3 columns are | Phi (deg) | I(phi) (1/cm) | std. dev. I(phi) (1/cm) |\r\n"
337        wfprintf refnum, formatStr, phival,inten,sig
338       
339        Close refnum
340       
341        SetDataFolder root:             //(redundant)
342       
343        //write confirmation of write operation to history area
344        Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
345        KillWaves/Z tempShortProto
346
347        Return(0)
348End
349
350//*****************
351// saves the data after all of the desired reduction steps (average options)
352// as a 2x expanded PNG file (approx 33kb)
353//
354Function SaveAsPNG(type,fullPath,dialog)
355        String type,fullPath
356        Variable dialog
357       
358        Variable refnum
359        if(dialog)
360                PathInfo/S catPathName
361                Open/D refnum as fullpath               //won't actually open the file
362                If(cmpstr(S_filename,"")==0)
363                        //user cancel, don't write out a file
364                        Close/A
365                        Abort "no data file was written"
366                Endif
367                fullpath = S_filename
368                //Print "dialog fullpath = ",fullpath
369        Endif
370       
371        //cleanup the filename passed in from Protocol...
372        String oldStr="",newStr="",pathStr=""
373        oldStr=GetFileNameFromPathNoSemi(fullPath)      //just the filename
374        pathStr=GetPathStrFromfullName(fullPath)        //just the path
375       
376        newStr = CleanupName(oldStr, 0 )                                //filename with _EXT rather than .EXT
377        fullPath=pathStr+newStr+".png"                          //tack on the png extension
378       
379        print "type=",type
380        //graph the current data and save a little graph
381        Wave data =  $("root:"+type+":data")
382        Wave q_x_axis = $"root:myGlobals:q_x_axis"
383        Wave q_y_axis = $"root:myGlobals:q_y_axis"
384        Wave NIHColors = $"root:myGlobals:NIHColors"
385       
386        NewImage/F data
387        DoWindow/C temp_PNG
388        ModifyImage data cindex= NIHColors
389        AppendToGraph/R q_y_axis
390        ModifyGraph tkLblRot(right)=90,lowTrip(right)=0.001
391        AppendToGraph/T q_x_axis
392        ModifyGraph lowTrip(top)=0.001,standoff=0,mode=2
393        ModifyGraph fSize(right)=9,fSize(top)=9,btLen=3
394       
395//      ModifyGraph nticks=0
396       
397//      WaveStats/Q data
398//      ScaleColorsToData(V_min, V_max, NIHColors)
399
400// ***comment out for DEMO_MODIFIED version
401        SavePict/Z/E=-5/B=144 as fullPath                       //PNG at 2x screen resolution
402//***
403
404        Print "Saved graphic as ",newStr+".png"
405        DoWindow/K temp_PNG
406End
407
408//****************
409//Testing only , not called
410Proc Fast_ASCII_2D_Export(type,term)
411        String type,term
412        Prompt type,"2-D data type for Export",popup,"SAM;EMP;BGD;DIV;COR;CAL;RAW;ABS;MSK;"
413        Prompt term,"line termination",popup,"CR;LF;CRLF;"
414       
415        //terminator is currently ignored
416        Fast2dExport(type,"",1)
417       
418End
419
420//the default termination for the platform is used...
421//if RAW export, sets "dummy" protocol to "RAW data export"
422Function Fast2dExport(type,fullpath,dialog)
423        String type,fullpath
424        Variable dialog         //=1 will present dialog for name
425       
426        String destStr="",ave="C",typeStr=""
427        Variable step=1,refnum
428        destStr = "root:"+type
429       
430        //must select the linear_data to export
431        // can't export log data if there are -ve intensities from a subtraction
432        NVAR isLog = $(destStr+":gIsLogScale")
433        if(isLog==1)
434                typeStr = ":linear_data"
435        else
436                typeStr = ":data"
437        endif
438       
439        NVAR pixelsX = root:myGlobals:gNPixelsX
440        NVAR pixelsY = root:myGlobals:gNPixelsY
441       
442        Wave data=$(destStr+typeStr)
443        WAVE intw=$(destStr + ":integersRead")
444        WAVE rw=$(destStr + ":realsRead")
445        WAVE/T textw=$(destStr + ":textRead")
446
447        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
448        String rawTag=""
449        if(cmpstr(type,"RAW")==0)
450                Make/O/T/N=8 proto={"none","none","none","none","none","none","none","none"}
451                RawTag = "RAW Data File: "     
452        else
453                Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
454        endif
455        SVAR samFiles = $("root:"+type+":fileList")
456        //check each wave - MUST exist, or will cause a crash
457        If(!(WaveExists(data)))
458                Abort "data DNExist AsciiExport()"
459        Endif
460        If(!(WaveExists(intw)))
461                Abort "intw DNExist AsciiExport()"
462        Endif
463        If(!(WaveExists(rw)))
464                Abort "rw DNExist AsciiExport()"
465        Endif
466        If(!(WaveExists(textw)))
467                Abort "textw DNExist AsciiExport()"
468        Endif
469        If(!(WaveExists(proto)))
470                Abort "current protocol wave DNExist AsciiExport()"
471        Endif
472       
473        if(dialog)
474                PathInfo/S catPathName
475                fullPath = DoSaveFileDialog("Save data as")
476                If(cmpstr(fullPath,"")==0)
477                        //user cancel, don't write out a file
478                        Close/A
479                        Abort "no data file was written"
480                Endif
481                //Print "dialog fullpath = ",fullpath
482        Endif
483       
484/////////
485        Variable numTextLines=18
486        Make/O/T/N=(numTextLines) labelWave
487        labelWave[0] = "FILE: "+textw[0]+"   CREATED: "+textw[1]
488        labelWave[1] = "LABEL: "+textw[6]
489        labelWave[2] = "MON CNT   LAMBDA(A)   DET_OFF(cm)   DET_DIST(m)   TRANS   THICK(cm)"
490        labelWave[3] = num2str(rw[0])+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
491        labelWave[3] += "     "+num2str(rw[4])+"     "+num2str(rw[5])
492        labelWave[4] = "BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP  "
493        labelWave[5] = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"  "+num2str(rw[24])+"  "
494        labelWave[5] += num2str(rw[25])+"  "+num2str(rw[27])+"  "+num2str(rw[21])+"  "+textW[9]
495        labelWave[6] =  "SAM: "+rawTag+samFiles
496        labelWave[7] =  "BGD: "+proto[0]
497        labelWave[8] =  "EMP: "+proto[1]
498        labelWave[9] =  "DIV: "+proto[2]
499        labelWave[10] =  "MASK: "+proto[3]
500        labelWave[11] =  "ABS Parameters (3-6): "+proto[4]
501        labelWave[12] = "Average Choices: "+proto[5]
502        labelWave[13] = ""
503        labelWave[14] = "*** Data written from "+type+" folder and may not be a fully corrected data file ***"
504        labelWave[15] = "The detector image is a standard X-Y coordinate system"
505        labelWave[16] = "Data is written by row, starting with Y=1 and X=(1->128)"
506        labelWave[17] = "ASCII data created " +date()+" "+time()
507        //strings can be too long to print-- must trim to 255 chars
508        Variable ii
509        for(ii=0;ii<numTextLines;ii+=1)
510                labelWave[ii] = (labelWave[ii])[0,240]
511        endfor
512//      If(cmpstr(term,"CR")==0)
513//              termStr = "\r"
514//      Endif
515//      If(cmpstr(term,"LF")==0)
516//              termStr = "\n"
517//      Endif
518//      If(cmpstr(term,"CRLF")==0)
519//              termStr = "\r\n"
520//      Endif
521       
522        Duplicate/O data,spWave         
523        Redimension/S/N=(pixelsX*pixelsY) spWave                //single precision (/S)
524       
525        //not demo- compatible, but approx 6x faster!!
526//      Save/G/M="\r\n" labelWave,spWave as fullPath            // /M=termStr specifies terminator
527        ///
528       
529        Open refNum as fullpath
530        wfprintf refNum,"%s\r\n",labelWave
531        fprintf refnum,"\r\n"
532        wfprintf refNum,"%g\r\n",spWave
533        Close refNum
534       
535        Killwaves/Z spWave,labelWave            //don't delete proto!
536       
537        Print "2D ASCII File written: ", GetFileNameFromPathNoSemi(fullPath)
538        return(0)
539End
540
541//**********************
542// Resolution calculation - used by the averaging routines
543// to calculate the resolution function at each q-value
544// - the return value is not used
545//
546// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
547// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
548//
549// - 21 MAR 07 uses projected BS diameter on the detector
550//
551Function/S getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
552        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
553        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
554       
555        //lots of calculation variables
556        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
557        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
558
559        //Constants
560        //Variable del_r = .1
561        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
562        Variable g = 981.0                              //gravity acceleration [cm/s^2]
563
564        String results
565        results ="Failure"
566
567        //rename for working variables,  these must be gotten from global
568        //variables
569
570//      Variable wLam, wLW, wL1, wL2, wS1, wS2
571//      Variable wBS, wDDet, wApOff
572//      wLam = lambda
573//      wLW = lambdaWidth
574//      wDDet = DDet
575//      wApOff = apOff
576        S1 *= 0.5*0.1                   //convert to radius and [cm]
577        S2 *= 0.5*0.1
578
579        L1 *= 100.0                     // [cm]
580        L1 -= apOff                             //correct the distance
581
582        L2 *= 100.0
583        L2 += apOff
584        del_r *= 0.1                            //width of annulus, convert mm to [cm]
585       
586        BS *= 0.5*0.1                   //nominal BS diameter passed in, convert to radius and [cm]
587        // 21 MAR 07 SRK - use the projected BS diameter, based on a point sample aperture
588        Variable LB
589        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
590        BS = bs + bs*lb/(l2-lb)         //adjusted diameter of shadow from parallax
591       
592        //Start resolution calculation
593        a2 = S1*L2/L1 + S2*(L1+L2)/L1
594        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
595        lp = 1.0/( 1.0/L1 + 1.0/L2)
596
597        v_lambda = lambdaWidth^2/6.0
598        v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2
599        v_d = (DDet/2.3548)^2 + del_r^2/12.0
600        vz = vz_1 / lambda
601        yg = 0.5*g*L2*(L1+L2)/vz^2
602        v_g = 2.0*yg^2*v_lambda
603
604        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
605        delta = 0.5*(BS - r0)^2/v_d
606
607        if (r0 < BS)
608                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
609        else
610                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
611        endif
612
613        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
614        if (fSubS <= 0.0)
615                fSubS = 1.e-10
616        endif
617        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
618        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
619
620        rmd = fr*r0
621        v_r1 = v_b + fv*v_d +v_g
622
623        rm = rmd + 0.5*v_r1/rmd
624        v_r = v_r1 - 0.5*(v_r1/rmd)^2
625        if (v_r < 0.0)
626                v_r = 0.0
627        endif
628        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
629        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
630
631        results = "success"
632        Return results
633End
634
635//ASCII export of data as 3-columns qx-qy-Intensity
636//limited header information?
637//
638// - creates the qx and qy data here, based on the data and header information
639//
640// Need to ensure that the data being exported is the linear copy of the dataset - check the global
641//
642Function QxQy_Export(type,fullpath,dialog)
643        String type,fullpath
644        Variable dialog         //=1 will present dialog for name
645       
646        String destStr="",typeStr=""
647        Variable step=1,refnum
648        destStr = "root:"+type
649       
650        //must select the linear_data to export
651        NVAR isLog = $(destStr+":gIsLogScale")
652        if(isLog==1)
653                typeStr = ":linear_data"
654        else
655                typeStr = ":data"
656        endif
657       
658        NVAR pixelsX = root:myGlobals:gNPixelsX
659        NVAR pixelsY = root:myGlobals:gNPixelsY
660       
661        Wave data=$(destStr+typeStr)
662        WAVE intw=$(destStr + ":integersRead")
663        WAVE rw=$(destStr + ":realsRead")
664        WAVE/T textw=$(destStr + ":textRead")
665
666        SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
667        String rawTag=""
668        if(cmpstr(type,"RAW")==0)
669                Make/O/T/N=8 proto={"none","none","none","none","none","none","none","none"}
670                RawTag = "RAW Data File: "     
671        else
672                Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
673        endif
674        SVAR samFiles = $("root:"+type+":fileList")
675        //check each wave - MUST exist, or will cause a crash
676        If(!(WaveExists(data)))
677                Abort "data DNExist QxQy_Export()"
678        Endif
679        If(!(WaveExists(intw)))
680                Abort "intw DNExist QxQy_Export()"
681        Endif
682        If(!(WaveExists(rw)))
683                Abort "rw DNExist QxQy_Export()"
684        Endif
685        If(!(WaveExists(textw)))
686                Abort "textw DNExist QxQy_Export()"
687        Endif
688        If(!(WaveExists(proto)))
689                Abort "current protocol wave DNExist QxQy_Export()"
690        Endif
691       
692        if(dialog)
693                PathInfo/S catPathName
694                fullPath = DoSaveFileDialog("Save data as")
695                If(cmpstr(fullPath,"")==0)
696                        //user cancel, don't write out a file
697                        Close/A
698                        Abort "no data file was written"
699                Endif
700                //Print "dialog fullpath = ",fullpath
701        Endif
702       
703/////////
704        Variable numTextLines=18
705        Make/O/T/N=(numTextLines) labelWave
706        labelWave[0] = "FILE: "+textw[0]+"   CREATED: "+textw[1]
707        labelWave[1] = "LABEL: "+textw[6]
708        labelWave[2] = "MON CNT   LAMBDA (A)  DET_OFF(cm)   DET_DIST(m)   TRANS   THICK(cm)"
709        labelWave[3] = num2str(rw[0])+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
710        labelWave[3] += "     "+num2str(rw[4])+"     "+num2str(rw[5])
711        labelWave[4] = "BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP  "
712        labelWave[5] = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
713        labelWave[5] += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+textW[9]
714        labelWave[6] =  "SAM: "+rawTag+samFiles
715        labelWave[7] =  "BGD: "+proto[0]
716        labelWave[8] =  "EMP: "+proto[1]
717        labelWave[9] =  "DIV: "+proto[2]
718        labelWave[10] =  "MASK: "+proto[3]
719        labelWave[11] =  "ABS Parameters (3-6): "+proto[4]
720        labelWave[12] = "Average Choices: "+proto[5]
721        labelWave[13] = ""
722        labelWave[14] = "*** Data written from "+type+" folder and may not be a fully corrected data file ***"
723        labelWave[15] = "Data columns are Qx - Qy - I(Qx,Qy)"
724        labelWave[16] = ""
725        labelWave[17] = "ASCII data created " +date()+" "+time()
726        //strings can be too long to print-- must trim to 255 chars
727        Variable ii,jj
728        for(ii=0;ii<numTextLines;ii+=1)
729                labelWave[ii] = (labelWave[ii])[0,240]
730        endfor
731//      If(cmpstr(term,"CR")==0)
732//              termStr = "\r"
733//      Endif
734//      If(cmpstr(term,"LF")==0)
735//              termStr = "\n"
736//      Endif
737//      If(cmpstr(term,"CRLF")==0)
738//              termStr = "\r\n"
739//      Endif
740       
741        Duplicate/O data,qx_val,qy_val,z_val
742        Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val
743        MyMat2XYZ(data,qx_val,qy_val,z_val)             //x and y are [p][q] indexes, not q-vals yet
744       
745        qx_val = CalcQx(qx_val+1,rw[16],rw[18],rw[26])          //+1 converts to detector coordinate system
746        qy_val = CalcQy(qy_val+1,rw[17],rw[18],rw[26])
747
748        //not demo-compatible, but approx 8x faster!!           
749//      Save/G/M="\r\n" labelWave,qx_val,qy_val,z_val as fullpath       // /M=termStr specifies terminator
750        //
751       
752        Open refNum as fullpath
753        wfprintf refNum,"%s\r\n",labelWave
754        fprintf refnum,"\r\n"
755        wfprintf refNum,"%8g\t%8g\t%8g\r\n",qx_val,qy_val,z_val
756        Close refNum
757       
758        Killwaves/Z spWave,labelWave,qx_val,qy_val,z_val
759       
760        Print "QxQy_Export File written: ", GetFileNameFromPathNoSemi(fullPath)
761        return(0)
762
763End
764
765
766Function MyMat2XYZ(mat,xw,yw,zw)
767        WAVE mat,xw,yw,zw
768
769        NVAR pixelsX = root:myGlobals:gNPixelsX
770        NVAR pixelsY = root:myGlobals:gNPixelsY
771       
772        xw= mod(p,pixelsX)              // X varies quickly
773        yw= floor(p/pixelsY)    // Y varies slowly
774        zw= mat(xw[p])(yw[p])
775
776End
777
778//converts xyz triple to a matrix
779//MAJOR assumption is that the x and y-spacings are LINEAR
780// (ok for small-angle approximation)
781//
782// currently unused
783//
784Function LinXYZToMatrix(xw,yw,zw,matStr)
785        WAVE xw,yw,zw
786        String matStr
787       
788        NVAR pixelsX = root:myGlobals:gNPixelsX
789        NVAR pixelsY = root:myGlobals:gNPixelsY
790        //mat is "zw" redimensioned to a matrix
791        Make/O/N=(pixelsX*pixelsY) $matStr
792        WAVE mat=$matStr
793        mat=zw
794        Redimension/N=(pixelsX,pixelsY) mat
795        WaveStats/Q xw
796        SetScale/I x, V_min, V_max, "",mat
797        WaveStats/Q yw
798        SetScale/I y, V_min, V_max, "",mat
799       
800        Display;Appendimage mat
801        ModifyGraph lowTrip=0.0001
802        ModifyGraph width={Plan,1,bottom,left},height={Plan,1,left,bottom}
803        ModifyImage $matStr ctab={*,*,YellowHot,0}
804       
805        return(0)
806End
Note: See TracBrowser for help on using the repository browser.