source: sans/Dev/trunk/NCNR_User_Procedures/SANS/Reduction/FIT_Ops.ipf @ 399

Last change on this file since 399 was 359, checked in by srkline, 15 years ago

Replaced the set of example data for SANS reduction with silica in D2O, all with nice sample labels.

Updated the help file for SANS reduction to use these files.

File size: 37.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.0
4
5//**************************
6//procedures for creating and initializing the FIT panel
7//global variables (numerical only) are kept in root:myGlobals:FIT folder,
8//created as needed
9//
10//also contains code for the FIT_RPA panel and its fitting procedures
11//
12//**************************
13//
14
15//main procedures to open the panel, initializing the data folder and global variables
16//as necessary. All are kept in a :FIT subfolder to avoid overlap with other variables
17//
18// To use any of the fit functions in the FIT panel, I(Q) data must already
19// be loaded into memory (using "plot" from the 1-D data operations
20// ** this may be useful to change in the future, replacing the 3 popups
21// with a list box - allowing the user to pick/load the data from the fit panel
22// and not offering any choice of q/i/s waves to use. (more consistent with the operation
23// of the FIT/RPA panel)
24//
25Proc OpenFitPanel()
26        If(WinType("FitPanel") == 0)
27                //create the necessary data folder
28                NewDataFolder/O root:myGlobals:FIT
29                //initialize the values
30                Variable/G root:myGlobals:FIT:gLolim = 0.02
31                Variable/G root:myGlobals:FIT:gUplim = 0.04
32                Variable/G root:myGlobals:FIT:gExpA = 1
33                Variable/G root:myGlobals:FIT:gExpB = 1
34                Variable/G root:myGlobals:FIT:gExpC = 1
35                Variable/G root:myGlobals:FIT:gBack = 0
36                String/G root:myGlobals:FIT:gDataPopList = "none"
37                FitPanel()
38        else
39                //window already exists, just bring to front for update
40                DoWindow/F FitPanel
41                CheckBox check0,value=0         //deselect the checkbox to use cursors
42        endif
43        //pop the file menu
44        FIT_FilePopMenuProc("",1,"")
45End
46
47//the actual window recreation macro to draw the fit panel. Globals and data folder must
48// already be initialized
49Window FitPanel()
50        String angst = root:myGlobals:gAngstStr
51        PauseUpdate; Silent 1           // building window...
52        NewPanel /W=(461,46,735,455)/K=1
53        ModifyPanel cbRGB=(32768,54615,65535), fixedSize=1
54        SetDrawLayer UserBack
55        DrawText 56,20,"Select Experimental Data"
56        DrawText 66,138,"q-range to fit ("+angst+"^-1)"
57        DrawText 42,268,"Select the y and x-axis scaling"
58        DrawLine 1,21,271,21
59        DrawLine -1,272,273,272
60        DrawLine -1,140,272,140
61//      PopupMenu xwave,pos={13,31},size={158,19},title="q-values"
62//      PopupMenu xwave,help={"Select the experimental q-values. If no data is available, use Macro LoadQISData or LoadQSIGData."}
63//      PopupMenu xwave,mode=1,value= #"WaveList(\"*q\",\";\",\"\")"
64        PopupMenu ywave,pos={13,40},size={154,19},title="Data File"
65        PopupMenu ywave,help={"Select the experimental intensity values"}
66        PopupMenu ywave,mode=1,value=root:myGlobals:FIT:gDataPopList,proc=FIT_FilePopMenuProc
67        Button loadButton,pos={25,80},size={130,20},proc=FIT_Load_Proc,title="Load and Plot File"
68        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."}
69        Button helpButton,pos={200,80},size={25,20},proc=showFITHelp,title="?"
70        Button helpButton,help={"Show help file for linearized fitting"}
71//      PopupMenu ewave,pos={13,87},size={186,19},title="Std. Deviation"
72//      PopupMenu ewave,help={"Select the standard deviation of the intensity"}
73//      PopupMenu ewave,mode=1,value= #"WaveList(\"*s\",\";\",\"\")"
74        PopupMenu ymodel,pos={20,281},size={76,19},title="y-axis"
75        PopupMenu ymodel,help={"This popup selects how the y-axis will be linearized based on the chosen data"}
76        PopupMenu ymodel,mode=1,value= #"\"I;log(I);ln(I);1/I;I^a;Iq^a;I^a q^b;1/sqrt(I);ln(Iq);ln(Iq^2)\""
77        Button GoFit,pos={60,367},size={70,20},proc=DispatchModel,title="Do the Fit"
78        Button GoFit,help={"This button will do the specified fit using the selections in this panel"}
79        Button DoneButton,pos={180,367},size={50,20},proc=FITDoneButton,title="Done"
80        Button DoneButton,help={"This button will close the panel and the associated graph"}
81        SetVariable lolim,pos={64,147},size={134,17},title="Lower Limit"
82        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"}
83        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FIT:gLolim
84        SetVariable uplim,pos={63,169},size={134,17},title="Upper Limit"
85        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"}
86        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FIT:gUplim
87        SetVariable expa,pos={13,311},size={80,17},title="pow \"a\""
88        SetVariable expa,help={"This sets the exponent \"a\" for some y-axis formats. The value is ignored if the model does not use an adjustable exponent"}
89        SetVariable expa,limits={-2,10,0},value= root:myGlobals:FIT:gExpA
90        SetVariable expb,pos={98,311},size={80,17},title="pow \"b\""
91        SetVariable expb,help={"This sets the exponent \"b\" for some x-axis formats. The value is ignored if the model does not use an adjustable exponent"}
92        SetVariable expb,limits={0,10,0},value= root:myGlobals:FIT:gExpB
93        PopupMenu xmodel,pos={155,280},size={79,19},title="x-axis"
94        PopupMenu xmodel,help={"This popup selects how the x-axis will be linearized given the chosen data"}
95        PopupMenu xmodel,mode=1,value= #"\"q;log(q);q^2;q^c\""
96        CheckBox check0,pos={18,223},size={240,20},title="Use cursor range from FitWindow"
97        CheckBox check0,help={"Checking this will perform a fit between the cursors on the graph in FitWindow and ignore the numerical limits typed above"},value=0
98        SetVariable back,pos={70,338},size={139,17},title="background"
99        SetVariable back,help={"This constant background value will be subtracted from the experimental intensity before fitting is done"}
100        SetVariable back,limits={-Inf,Inf,0},value= root:myGlobals:FIT:gBack
101        SetVariable expc,pos={182,310},size={80,17},title="pow \"c\""
102        SetVariable expc,help={"This sets the exponent \"c\" for some x-axis formats. The value is ignored if the model does not use \"c\" as an adjustable exponent"}
103        SetVariable expc,limits={-10,10,0},value= root:myGlobals:FIT:gExpC
104        Button sh_all,pos={65,193},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range"
105        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."}
106EndMacro
107
108
109Proc FITDoneButton(ctrlName): ButtonControl
110        String ctrlName
111        DoWindow/K FitWindow
112        DoWindow/K FitPanel
113end
114
115Proc showFITHelp(ctrlName): ButtonControl
116        String ctrlName
117        DisplayHelpTopic/K=1 "SANS Data Reduction Tutorial[Fit Lines to Your Data]"
118        if(V_flag !=0)
119                DoAlert 0,"The SANS Data Reduction Tutorial Help file could not be found"
120        endif
121end
122
123//Loads the selected file for fitting
124//graphs the data as needed
125Proc FIT_Load_Proc(ctrlName): ButtonControl
126        String ctrlName
127       
128        //Load the data
129        String tempName="",partialName=""
130        Variable err
131        ControlInfo $"ywave"
132        //find the file from the partial filename
133        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) )
134                //null selection, or "none" from any popup
135                Abort "no file selected in popup menu"
136        else
137                //selection not null
138                partialName = S_value
139                //Print partialName
140        Endif
141        //get a valid file based on this partialName and catPathName
142        tempName = FindValidFilename(partialName)
143
144        //prepend path to tempName for read routine
145        PathInfo catPathName
146       
147        tempName = S_path + tempName
148       
149        //load in the data (into the root directory)
150        LoadOneDDataWithName(tempName)
151        //Print S_fileName
152        //Print tempName
153       
154        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
155
156        tempName=cleanLastFileName+"_q"
157        Duplicate/O $tempName xAxisWave
158        tempName=cleanLastFileName+"_i"
159        Duplicate/O $tempName yAxisWave
160        tempName=cleanLastFileName+"_s"
161        Duplicate/O $tempName yErrWave
162
163        //Plot, and adjust the scaling to match the axis scaling set by the popups
164        Rescale_Data()
165       
166End
167
168//gets a valid file list (simply not the files with ".SAn" in the name)
169//
170Function FIT_FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
171        String ctrlName
172        Variable popNum
173        String popStr
174       
175        String tempStr=ReducedDataFileList(ctrlName)
176        if(strlen(tempStr)==0)
177                tempStr = "Pick the data path"
178        Endif
179        String/G root:myGlobals:FIT:gDataPopList =tempStr               //function is in NSORT.ipf
180        ControlUpdate ywave
181       
182End
183
184//direct porting of the fit program from the VAX, with no corrections and only
185//minor modifications of additional  linearizations and the option
186//to subtract a constant (q-independent) background value before doing any
187//of the fits. The original data on disk (and as loaded) is never modified, all
188//manipulation is done from a copy of the data.
189
190//button procedure to show the entire axis range of the data, rather than just
191//the fitted range, which is the default display after a fit is performed
192//
193Function ShowAllButtonProc(ctrlName) : ButtonControl
194        String ctrlName
195
196        //bring the FitWindow to the front and Autoscale the axes
197        DoWindow/F FitWindow
198        SetAxis/A
199End
200
201// function that takes the current dataset (already loaded)
202// and replots it based on the X/Y axis scaling selected in the popups
203// (does not fit the data)
204//
205Function Rescale_Data()
206       
207        //Scaling exponents and background value
208        Variable pow_a,pow_b,pow_c,bkg
209        ControlInfo/W=FitPanel expa
210        pow_a = V_value
211        ControlInfo/W=FitPanel expb
212        pow_b = V_value
213        ControlInfo/W=FitPanel expc
214        pow_c = V_value
215        ControlInfo/W=FitPanel back
216        bkg = V_value
217       
218//check for physical limits on exponent values
219// if bad values found, alert, and reset to good values so the rescaling can continue
220        NVAR gA = root:myGlobals:FIT:gExpA
221        NVAR gB = root:myGlobals:FIT:gExpB
222        NVAR gC = root:myGlobals:FIT:gExpC
223        if((pow_a < -2) || (pow_a > 10))
224                DoAlert 0,"Exponent a must be in the range (-2,10) - the exponent a has been reset to 1"
225                gA = 1
226        endif
227        if((pow_b < 0) || (pow_b > 10))
228                DoAlert 0,"Exponent b must be in the range (0,10) - the exponent b has been reset to 1"
229                gB = 1
230        endif
231        //if q^c is the x-scaling, c must be be within limits and also non-zero
232        ControlInfo/W=FitPanel xModel
233        If (cmpstr("q^c",S_Value) == 0)
234                if(pow_c == 0)
235                        DoAlert 0,"Exponent c must be non-zero, c has been reset to 1"
236                        gC = 1
237                endif
238                if((pow_c < -10) || (pow_c > 10))
239                        DoAlert 0,"Exponent c must be in the range (-10,10), c has been reset to 1"
240                        gC = 1
241                endif
242        endif
243       
244        //do the rescaling of the data
245        // get the current experimental q, I, and std dev. waves (as they would be loaded )
246
247        //ControlInfo/W=FitPanel ywave
248        //get the filename from the global as it's loaded, rather from the popup - as version numbers
249        // do cause problems here. This global is also used later in this function
250        SVAR gLastFileName = root:myGlobals:gLastFileName
251       
252        Wave xw = $( CleanupName((gLastFileName + "_q"),0) )
253        Wave yw = $( CleanupName((gLastFileName + "_i"),0) )
254        Wave ew = $( CleanupName((gLastFileName + "_s"),0) )
255       
256        //variables set for each model to control look of graph
257        Variable xlow,xhigh,ylow,yhigh,yes_cursors
258        String xlabel,ylabel,xstr,ystr
259        //check for proper y-scaling selection, make the necessary waves
260        ControlInfo/W=FitPanel yModel
261        ystr = S_Value
262//      print "ystr = ",ystr
263        do
264                // make the new yaxis waves, including weighting wave
265                Duplicate/O yw yAxisWave,yErrWave,yWtWave,residWave
266                //subtract the background value from yAxisWave before doing any rescaling
267                yAxisWave = yw - bkg
268               
269                If (cmpstr("I",S_Value) == 0)
270                        SetScale d 0,0,"1/cm",yAxisWave
271                        yErrWave = ew
272                        yWtWave = 1/yErrWave
273                        yAxisWave = yAxisWave
274                        ylabel = "I(q)"
275                        break   
276                endif
277                If (cmpstr("ln(I)",S_Value) == 0)
278                        SetScale d 0,0,"",yAxisWave
279                        yErrWave = ew/yAxisWave
280                        yWtWave = 1/yErrWave
281                        yAxisWave = ln(yAxisWave)
282                        ylabel = "ln(I)"
283                        break   
284                endif
285                If (cmpstr("log(I)",S_Value) == 0)
286                        SetScale d 0,0,"",yAxisWave
287                        yErrWave = ew/(2.30*yAxisWave)
288                        yWtWave = 1/yErrWave
289                        yAxisWave = log(yAxisWave)
290                        ylabel = "log(I)"
291                        break   
292                endif
293                If (cmpstr("1/I",S_Value) == 0)
294                        SetScale d 0,0,"",yAxisWave
295                        yErrWave = ew/yAxisWave^2
296                        yWtWave = 1/yErrWave
297                        yAxisWave = 1/yAxisWave
298                        ylabel = "1/I"
299                        break
300                endif
301                If (cmpstr("I^a",S_Value) == 0)
302                        SetScale d 0,0,"",yAxisWave
303                        yErrWave = ew*abs(pow_a*(yAxisWave^(pow_a-1)))
304                        yWtWave = 1/yErrWave
305                        yAxisWave = yAxisWave^pow_a
306                        ylabel = "I^"+num2str(pow_a)
307                        break
308                endif
309                If (cmpstr("Iq^a",S_Value) == 0)
310                        SetScale d 0,0,"",yAxisWave
311                        yErrWave = ew*xw^pow_a
312                        yWtWave = 1/yErrWave
313                        yAxisWave = yAxisWave*xw^pow_a
314                        ylabel = "I*q^"+num2str(pow_a)
315                        break
316                endif
317                If (cmpstr("I^a q^b",S_Value) == 0)
318                        SetScale d 0,0,"",yAxisWave
319                        yErrWave = ew*abs(pow_a*(yAxisWave^(pow_a-1)))*xw^pow_b
320                        yWtWave = 1/yErrWave
321                        yAxisWave = yAxisWave^pow_a*xw^pow_b
322                        ylabel = "I^" + num2str(pow_a) + "q^"+num2str(pow_b)
323                        break
324                endif
325                If (cmpstr("1/sqrt(I)",S_Value) == 0)
326                        SetScale d 0,0,"",yAxisWave
327                        yErrWave = 0.5*ew*yAxisWave^(-1.5)
328                        yWtWave = 1/yErrWave
329                        yAxisWave = 1/sqrt(yAxisWave)
330                        ylabel = "1/sqrt(I)"
331                        break
332                endif
333                If (cmpstr("ln(Iq)",S_Value) == 0)
334                        SetScale d 0,0,"",yAxisWave
335                        yErrWave =ew/yAxisWave
336                        yWtWave = 1/yErrWave
337                        yAxisWave = ln(xw*yAxisWave)
338                        ylabel = "ln(q*I)"
339                        break
340                endif
341                If (cmpstr("ln(Iq^2)",S_Value) == 0)
342                        SetScale d 0,0,"",yAxisWave
343                        yErrWave = ew/yAxisWave
344                        yWtWave = 1/yErrWave
345                        yAxisWave = ln(xw*xw*yAxisWave)
346                        ylabel = "ln(I*q^2)"
347                        break
348                endif
349                //more ifs for each case
350               
351                // if selection not found, abort
352                DoAlert 0,"Y-axis scaling incorrect. Aborting"
353                Abort
354        while(0)        //end of "case" statement for y-axis scaling
355       
356       
357        //check for proper x-scaling selection
358        Variable low,high
359       
360        ControlInfo/W=FitPanel lolim
361        low = V_value
362        ControlInfo/W=FitPanel uplim
363        high = V_value
364        if ((high<low) || (high==low))
365                DoAlert 0,"Unphysical fitting limits - re-enter better values"
366                Abort
367        endif
368       
369        SVAR aStr = root:myGlobals:gAngstStr
370       
371        ControlInfo/W=FitPanel xModel
372        xstr = S_Value
373        do
374                // make the new yaxis wave
375                Duplicate/o xw xAxisWave
376                If (cmpstr("q",S_Value) == 0)   
377                        SetScale d 0,0,aStr+"^-1",xAxisWave
378                        xAxisWave = xw
379                        xlabel = "q"
380                        xlow = low
381                        xhigh = high
382                        break   
383                endif
384                If (cmpstr("q^2",S_Value) == 0)
385                        SetScale d 0,0,aStr+"^-2",xAxisWave
386                        xAxisWave = xw*xw
387                        xlabel = "q^2"
388                        xlow = low^2
389                        xhigh = high^2
390                        break   
391                endif
392                If (cmpstr("log(q)",S_Value) == 0)     
393                        SetScale d 0,0,"",xAxisWave
394                        xAxisWave = log(xw)
395                        xlabel = "log(q)"
396                        xlow = log(low)
397                        xhigh = log(high)
398                        break   
399                endif
400                If (cmpstr("q^c",S_Value) == 0)
401                        SetScale d 0,0,"",xAxisWave
402                        xAxisWave = xw^pow_c
403                        xlabel = "q^"+num2str(pow_c)
404                        xlow = low^pow_c
405                        xhigh = high^pow_c
406                        break
407                endif
408       
409                //more ifs for each case
410               
411                // if selection not found, abort
412                DoAlert 0,"X-axis scaling incorrect. Aborting"
413                Abort
414        while(0)        //end of "case" statement for x-axis scaling
415
416        //plot the data
417       
418        String cleanLastFileName = "root:"+CleanupName(gLastFileName,0)
419        If(WinType("FitWindow") == 0)
420                Display /W=(5,42,480,400)/K=1 yAxisWave vs xAxisWave
421                ModifyGraph mode=3,standoff=0,marker=8
422                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
423                ModifyGraph opaque(yAxisWave)=1
424                DoWindow/C FitWindow
425        else
426                //window already exists, just bring to front for update
427                DoWindow/F FitWindow
428                // remove old text boxes
429                TextBox/K/N=text_1
430                TextBox/K/N=text_2
431                TextBox/K/N=text_3
432        endif
433        SetAxis/A
434        ModifyGraph tickUnit=1          //suppress tick units in labels
435        TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
436        //clear the old fit from the window, if it exists
437        RemoveFromGraph/W=FitWindow/Z fit_yAxisWave
438       
439        // add the cursors if desired...       
440        //see if the user wants to use the data specified by the cursors - else use numerical values
441       
442        ControlInfo/W=FitPanel check0           //V_value = 1 if it is checked, meaning yes, use cursors
443        yes_cursors = V_value
444
445        DoWindow/F FitWindow
446        ShowInfo
447        if(yes_cursors)
448                xlow = xAxisWave[xcsr(A)]
449                xhigh = xAxisWave[xcsr(B)]
450                if(xlow > xhigh)
451                        xhigh = xlow
452                        xlow = xAxisWave[xcsr(B)]
453                endif
454//              Print xlow,xhigh
455        else
456                FindLevel/P/Q xAxisWave, xlow
457                if(V_flag == 1)                 //level NOT found
458                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
459                        //DoWindow/K FitWindow
460                        //Abort
461                endif
462                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
463                ylow = V_LevelX
464                FindLevel/P/Q xAxisWave, xhigh
465                if(V_flag == 1)
466                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
467                        //DoWindow/K FitWindow
468                        //Abort
469                endif
470                Cursor/P B, yAxisWave,trunc(V_LevelX)
471                yhigh = V_LevelX
472        endif   //if(V_value)
473        //SetAxis bottom,xlow,xhigh
474        //SetAxis left,ylow,yhigh
475        Label left ylabel
476        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
477
478End
479
480
481//button procedure that is activated to "DotheFit"
482//the panel is parsed for proper fitting limits
483// the appropriate linearization is formed (in the Rescale_Data() function)
484// and the fit is done,
485//and the results are plotted
486// function works in root level data folder (where the loaded 1-d data will be)
487Function DispatchModel(GoFit) : ButtonControl
488        String GoFit
489
490        //check for the FitWindow - to make sure that there is data to fit
491        If(WinType("FitWindow") == 0)           //if the window doesn't exist
492                Abort "You must Load and Plot a File before fitting the data"
493        endif
494        // rescale the data, to make sure it's as selected on the panel
495        Rescale_Data()
496       
497        // now go do the fit
498       
499// get the current low and high q values for fitting
500        Variable low,high
501       
502        ControlInfo/W=FitPanel lolim
503        low = V_value
504        ControlInfo/W=FitPanel uplim
505        high = V_value
506        if ((high<low) || (high==low))
507                DoAlert 0,"Unphysical fitting limits - re-enter better values"
508                Abort
509        endif
510
511        //try including residuals on the graph /R=residWave, explicitly place on new axis
512        //if only /R used, residuals are automatically placed on graph
513       
514        CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D 
515        //CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave  /R /D 
516        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
517        ModifyGraph lsize(fit_yAxisWave)=2
518// annotate graph, filtering out special cases of Guinier fits
519// Text Boxes must be used, since ControlBars on graphs DON'T print out
520       
521        // need access to Global wave, result of fit
522        //ystr and xstr are the axis strings - filter with a do-loop
523        String ystr="",xstr=""
524        SVAR gLastFileName = root:myGlobals:gLastFileName
525        SVAR aStr = root:myGlobals:gAngstStr
526       
527        //ControlInfo/W=FitPanel ywave
528        Wave xw = $( CleanupName((gLastFileName + "_q"),0) )
529        ControlInfo/W=FitPanel yModel
530        ystr = S_Value
531        ControlInfo/W=FitPanel xModel
532        xstr = S_Value
533       
534        WAVE W_coef=W_coef
535        WAVE W_sigma=W_sigma
536        String textstr_1,textstr_2,textstr_3 = ""
537        Variable rg,rgerr,minfit,maxfit
538       
539        textstr_1 = "Slope = " + num2str(W_coef[1]) + " ± " + num2str(W_sigma[1])
540        textstr_1 += "\rIntercept = " + num2str(W_coef[0]) + " ± " + num2str(W_sigma[0])
541        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
542       
543        minfit = xw[xcsr(A)]
544        maxfit = xw[xcsr(B)]
545        textstr_2 = "Qmin =  " + num2str(minfit)
546        textstr_2 += "\rQmax =  " + num2str(maxfit)
547       
548        //model-specific calculations - I(0), Rg, etc.
549        //put these in textstr_3, at bottom
550        do
551                If (cmpstr("I",ystr) == 0)
552                        textstr_3 = "I(q=0) =  "  + num2str(W_coef[0])
553                        break   
554                endif
555                If (cmpstr("ln(I)",ystr) == 0)
556                        textstr_3 = "I(q=0) =  "  + num2str(exp(W_coef[0]))
557                        if(cmpstr("q^2",xstr) == 0)     //then a Guinier plot for a sphere (3-d)
558                                rg = sqrt(-3*W_coef[1])
559                                rgerr = 3*W_sigma[1]/(2*rg)
560                                textstr_3 += "\rRg ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
561                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
562                                break
563                        endif
564                        break   
565                endif
566                If (cmpstr("log(I)",ystr) == 0)
567                        if(cmpstr("log(q)",xstr) !=0 )  //extrapolation is nonsense
568                                textstr_3 = "I(q=0) =  "  + num2str(10^(W_coef[0]))
569                        endif
570                        break   
571                endif
572                If (cmpstr("1/I",ystr) == 0)
573                        textstr_3 = "I(q=0) =  "  + num2str(1/W_coef[0])+ " ± " + num2str(1/(W_coef[0] - W_sigma[0])-1/W_coef[0])
574                        break
575                endif
576                If (cmpstr("I^a",ystr) == 0)
577                        //nothing
578                        break
579                endif
580                If (cmpstr("Iq^a",ystr) == 0)
581                        //nothing
582                        break
583                endif
584                If (cmpstr("I^a q^b",ystr) == 0)
585                        //nothing
586                        break
587                endif
588                If (cmpstr("1/sqrt(I)",ystr) == 0)
589                        textstr_3 = "I(q=0) =  "  + num2str((W_coef[0])^2)
590                        break
591                endif
592                If (cmpstr("ln(Iq)",ystr) == 0)
593                        //nothing
594                        if(cmpstr("q^2",xstr) == 0)     //then a x-sect Guinier plot for a rod (2-d)
595                                // rg now is NOT the radius of gyration, but the x-sect DIAMETER
596                                rg = 4*sqrt(-W_coef[1])
597                                rgerr = 8*W_sigma[1]/rg
598                                textstr_3 = "Rod diameter ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
599                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
600                                break
601                        endif
602                        break
603                endif
604                If (cmpstr("ln(Iq^2)",ystr) == 0)
605                        //nothing
606                        if(cmpstr("q^2",xstr) == 0)     //then a 1-d Guinier plot for a sheet
607                                // rg now is NOT the radius of gyration, but the thickness
608                                rg = sqrt(-12*W_coef[1])
609                                rgerr = 6*W_sigma[1]/(2*rg)
610                                textstr_3 = "Platelet thickness ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
611                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
612                                break
613                        endif
614                        break
615                endif
616               
617        while(0)
618        //kill the old textboxes, if they exist
619        TextBox/W=FitWindow/K/N=text_1
620        TextBox/W=FitWindow/K/N=text_2
621        TextBox/W=FitWindow/K/N=text_3
622        // write the new text boxes
623        TextBox/W=FitWindow/N=text_1/A=LT textstr_1
624        TextBox/W=FitWindow/N=text_2/A=LC textstr_2
625        If (cmpstr("",textstr_3) != 0)          //only display textstr_3 if it isn't null
626                TextBox/W=FitWindow/N=text_3/A=LB textstr_3
627        endif
628       
629        //adjust the plot range to reflect the actual fitted range
630        //cursors are already on the graph, done by Rescale_Data()
631        AdjustAxisToCursors()
632       
633End
634
635// adjusts both the x-axis scaling  and y-axis scaling to the cursor range
636// **cursors are already on the graph, done by Rescale_Data()
637//
638// will expand the scale to show an extra 5 points in each direction (if available)
639Function AdjustAxisToCursors()
640
641        DoWindow/F FitWindow
642        WAVE xAxisWave = root:xAxisWave
643        WAVE yAxisWave = root:yAxisWave
644        Variable xlow,xhigh,ylow,yhigh,yptlow,ypthigh
645        Variable extraPts = 5, num=numpnts(xAxisWave)
646       
647        String csrA = CsrInfo(A ,"FitWindow")
648        String csrB = CsrInfo(B ,"FitWindow")
649       
650        //x-levels, these are monotonic
651        Variable ptLow,ptHigh,tmp
652        ptLow = NumberByKey("POINT", csrA ,":" ,";")
653        ptHigh = NumberByKey("POINT", csrB ,":" ,";")
654        if(ptLow > ptHigh)
655                tmp= ptLow
656                ptLow=ptHigh
657                ptHigh=tmp
658        endif
659
660        // keep extended point range in bounds
661        ptLow = (ptLow-extraPts) >= 0 ? ptLow-extraPts : 0
662        ptHigh = (ptHigh+extraPts) <= (num-1) ? ptHigh + extraPts : num-1
663       
664        xlow = xAxisWave[ptLow]
665        xhigh = xAxisWave[ptHigh]
666//old way
667//      xlow = xAxisWave[xcsr(A)]
668//      xhigh = xAxisWave[xcsr(B)]
669//      if(xlow > xhigh)
670//              xhigh = xlow
671//              xlow = xAxisWave[xcsr(B)]
672//      endif
673       
674        //y-levels (old way)
675//      FindLevel/P/Q xAxisWave, xlow
676//      if(V_flag == 1)                 //level NOT found
677//              DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
678//      endif
679//      yptlow = V_LevelX
680//      FindLevel/P/Q xAxisWave, xhigh
681//      if(V_flag == 1)
682//              DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
683//      endif
684//      ypthigh = V_LevelX
685
686//      Print xlow,xhigh,yptlow,ypthigh
687//      Print yAxisWave[yptlow],yAxisWave[ypthigh]
688       
689
690        // make sure ylow/high are in the correct order, since the slope could be + or -
691        yhigh = max(yAxisWave[ptlow],yAxisWave[pthigh])
692        ylow = min(yAxisWave[ptlow],yAxisWave[pthigh])
693       
694//      Print ptLow,ptHigh
695//      print xlow,xhigh
696//      print ylow,yhigh
697       
698        SetAxis bottom,xlow,xhigh
699        SetAxis left ylow,yhigh
700       
701End
702
703
704//****************************************
705//procedures for creating and initializing the FITRPA panel
706//global variables (numerical only) are kept in root:myGlobals:FITRPA folder,
707//created as needed
708//
709// very similar in function to the FIT panel
710//
711Proc OpenFitRPAPanel()
712        If(WinType("FitRPAPanel") == 0)
713                //create the necessary data folder
714                NewDataFolder/O root:myGlobals:FITRPA
715                //initialize the values
716                Variable/G root:myGlobals:FITRPA:gLolim = 0.02
717                Variable/G root:myGlobals:FITRPA:gUplim = 0.04
718                Variable/G root:myGlobals:FITRPA:gBack = 0
719                Variable/G root:myGlobals:FITRPA:gLambda = 6.0
720                PathInfo/S catPathName
721                String localpath = S_path
722                if (V_flag == 0)
723                //path does not exist - no folder selected
724                        String/G root:myGlobals:FITRPA:gPathStr = "no folder selected"
725                else
726                        String/G root:myGlobals:FITRPA:gPathStr = localpath
727                endif
728                String/G    root:myGlobals:FITRPA:gDataPopList = "none"
729                FitRPAPanel()
730        else
731                //window already exists, just bring to front for update
732                DoWindow/F FitRPAPanel
733        endif
734        //pop the menu
735        FilePopMenuProc("",1,"")
736End
737
738//used on the fit/rpa panel to select the path for the data
739// - automatically pops the file list after the new  path selection
740Function FITRPAPickPathButton(ctrlName) : ButtonControl
741        String ctrlName
742
743        Variable err = PickPath()               //sets global path value
744        SVAR pathStr = root:myGlobals:gCatPathStr
745       
746        //set the global string for NSORT to the selected pathname
747        String/G root:myGlobals:FITRPA:gPathStr = pathStr
748       
749        //call the popup menu proc's to re-set the menu choices
750        FilePopMenuProc("filePopup",1,"")
751       
752End
753
754//gets a valid file list (simply not the files with ".SAn" in the name)
755//
756Function FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
757        String ctrlName
758        Variable popNum
759        String popStr
760
761        String tempStr=ReducedDataFileList(ctrlName)
762        if(strlen(tempStr)==0)
763                tempStr = "Pick the data path"
764        Endif
765        String/G root:myGlobals:FITRPA:gDataPopList =   tempStr //function is in NSORT.ipf
766        ControlUpdate filePopup
767       
768End
769
770
771
772// window recreation macro to draw the fit/rpa panel
773//globals and data folders must be present before drawing panel
774//
775Window FitRPAPanel()
776
777        String angst = root:myGlobals:gAngstStr
778        PauseUpdate; Silent 1           // building window...
779        NewPanel /W=(250,266,591,579)/K=1
780        ModifyPanel cbRGB=(32768,54528,65280)
781        ModifyPanel fixedSize=1
782        SetDrawLayer UserBack
783        SetDrawEnv fstyle= 1
784        DrawText 81,19,"Select Experimental Data"
785        SetDrawEnv fstyle= 1
786        DrawText 97,102,"q-range to fit ("+angst+"^-1)"
787        SetDrawEnv fstyle= 1
788        DrawText 87,239,"Select the fit parameters"
789        SetDrawEnv fillpat= 0
790        DrawRect 1,103,338,224
791        SetDrawEnv fillpat= 0
792        DrawRect 1,20,337,83
793        SetDrawEnv fillpat= 0
794        DrawRect 2,241,337,275
795//      Button PathButton,pos={6,26},size={80,20},proc=FitRPAPickPathButton,title="Pick Path"
796//      Button PathButton,help={"Select the local path to the folder containing your SANS data"}
797//      SetVariable setPath,pos={95,29},size={240,17},title="Path:"
798//      SetVariable setPath,help={"The current path to the local folder with SANS data"}
799//      SetVariable setPath,fSize=10
800//      SetVariable setPath,limits={0,0,0},value= root:myGlobals:FITRPA:gPathStr
801        PopupMenu filePopup,pos={8,30},size={96,21},proc=FilePopMenuProc,title="Files"
802        PopupMenu filePopup,help={"Select the data file to load."}
803        PopupMenu filePopup,mode=5,popvalue="none",value= #"root:myGlobals:FITRPA:gDataPopList"
804        SetVariable lambda,pos={111,250},size={120,18},title="Lambda ("+angst+")"
805        SetVariable lambda,help={"This sets the wavelength for the multiple scattering corrections."}
806        SetVariable lambda,limits={0,10,0},value= root:myGlobals:FITRPA:gLambda
807        Button GoFit,pos={60,286},size={80,20},proc=DoFITRPA,title="Do the Fit"
808        Button GoFit,help={"This button will do the specified fit using the selections in this panel"}
809        SetVariable lolim,pos={82,113},size={134,28},title="Lower Limit"
810        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"}
811        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FITRPA:gLolim
812        SetVariable uplim,pos={80,140},size={134,28},title="Upper Limit"
813        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"}
814        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FITRPA:gUplim
815        CheckBox RPA_check0,pos={64,198},size={190,20},title="Use cursor range from FitWindow"
816        CheckBox RPA_check0,help={"Checking this will perform a fit between the cursors on the graph in FitWindow and ignore the numerical limits typed above"},value=0
817        PopupMenu model,pos={3,249},size={101,21},title="Standard"
818        PopupMenu model,help={"This popup selects which standard should be used to fit this data"}
819        PopupMenu model,mode=1,popvalue="B",value= #"\"B;C;AS\""
820        Button sh_all,pos={82,168},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range"
821        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."}
822        Button loadButton,pos={20,55},size={70,20},proc=FITRPA_Load_Proc,title="Load File"
823        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."}
824        Button helpButton,pos={270,55},size={25,20},proc=showFITHelp,title="?"
825        Button helpButton,help={"Show help file for RPA fitting"}
826        Button DoneButton,pos={200,286},size={50,20},proc=FITRPADoneButton,title="Done"
827        Button DoneButton,help={"This button will close the panel and the associated graph"}
828EndMacro
829
830
831Proc FITRPADoneButton(ctrlName) : ButtonControl
832        String ctrlName
833        DoWindow/K FitWindow
834        DoWindow/K FitRPAPanel
835end
836
837//dispatches the fit to the appropriate model
838//and the appropriate range, based on selections in the panel
839//
840Proc DoFITRPA(ctrlName) : ButtonControl
841        String ctrlName
842        //
843        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
844        String tempName
845
846        tempName=cleanLastFileName+"_q"
847        Duplicate/O $tempName xAxisWave
848        tempName=cleanLastFileName+"_i"
849        Duplicate/O $tempName yAxisWave
850        tempName=cleanLastFileName+"_s"
851        Duplicate/O $tempName yErrWave
852        Duplicate/O $tempName yWtWave
853        Duplicate/O $tempName residWave
854        yWtWave = 1/yErrWave
855       
856        String xlabel = "q (A^-1)"
857        String ylabel = "Intensity"
858        //Check to see if the FitWindow exists
859        //Plot the data in a FitWindow
860        If(WinType("FitWindow") == 0)
861                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
862                ModifyGraph mode=3,standoff=0,marker=8
863                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
864                DoWindow/C FitWindow
865                ShowInfo
866        else
867                //window already exists, just bring to front for update
868                DoWindow/F FitWindow
869                // remove old text boxes
870                TextBox/K/N=text_1
871                TextBox/K/N=text_2
872                TextBox/K/N=text_3
873        endif
874       
875        //see if the user wants to use the data specified by the cursors - else use numerical values
876        Variable xlow,xhigh,ylow,yhigh,yes_cursors
877        ControlInfo/W=FitRPAPanel RPA_check0            //V_value = 1 if it is checked, meaning yes, use cursors
878        yes_cursors = V_value
879
880        ControlInfo/W=FitRPAPanel lolim
881        xlow = V_value
882        ControlInfo/W=FitRPAPanel uplim
883        xhigh = V_value
884        if(yes_cursors)
885                xlow = xAxisWave[xcsr(A)]
886                xhigh = xAxisWave[xcsr(B)]
887                if(xlow > xhigh)
888                        xhigh = xlow
889                        xlow = xAxisWave[xcsr(B)]
890                endif
891//              Print "xlow,xhigh = ",xlow,xhigh
892                ylow = yAxisWave[xcsr(A)]
893                yhigh = yAxisWave[xcsr(B)]
894                if(ylow > yhigh)
895                        ylow=yhigh
896                        yhigh = yAxisWave[xcsr(A)]
897                endif
898        else
899                FindLevel/P/Q xAxisWave, xlow
900                if(V_flag == 1)                 //level NOT found
901                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
902                        DoWindow/K FitWindow
903                        Abort
904                endif
905                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
906                ylow = yAxisWave[V_LevelX]
907                FindLevel/P/Q xAxisWave, xhigh
908                if(V_flag == 1)
909                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
910                        DoWindow/K FitWindow
911                        Abort
912                endif
913                Cursor/P B, yAxisWave,trunc(V_LevelX)
914                yhigh = yAxisWave[V_LevelX]
915                if(ylow > yhigh)
916                        yhigh=ylow
917                        ylow = yAxisWave[V_levelX]
918                endif
919        endif   //if(V_value)
920        SetAxis bottom,xlow,xhigh
921       
922//      print "ylow,yhigh",ylow,yhigh
923       
924        //Get the rest of the data from the panel
925        //such as which standard, the wavelength
926        ControlInfo/W=FitRPAPanel model
927
928        //find the model name
929        String modelName = S_value
930       
931        Variable first_guess, seglength,iabs,iarb,thick
932        Make/D/O/N=2 fitParams
933       
934        seglength = 6.8
935       
936        first_guess = 1.0
937        fitParams[0] = first_guess
938        fitParams[1] = seglength
939
940        If (cmpstr(modelName,"B")==0)
941                iabs = BStandardFunction(fitParams,xlow)
942                fitParams[0] = yhigh/iabs
943//              Print fitParams[0],fitParams[1]
944                FuncFit BStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
945                iarb = BStandardFunction(fitParams, 0.0)
946                iabs = iarb/fitParams[0]
947                thick = 0.153
948        endif
949        If (cmpstr(modelName,"C")==0)
950                iabs = CStandardFunction(fitParams,xlow)
951                fitParams[0] = yhigh/iabs
952                FuncFit CStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
953                iarb = CStandardFunction(fitParams, 0.0)
954                iabs = iarb/fitParams[0]
955                thick= 0.153
956        endif
957        If (cmpstr(modelName,"AS")==0)
958                iabs = ASStandardFunction(fitParams,xlow)
959                fitParams[0] = yhigh/iabs
960                FuncFit ASStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
961                iarb = ASStandardFunction(fitParams, 0.0)
962                iabs = iarb/fitParams[0]
963                thick = 0.1
964        endif
965        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
966        Label left ylabel
967        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
968
969        ControlInfo/W=FitRPAPanel lambda
970       
971        Variable cor_mult = 1.0 + 2.2e-4*V_Value^2
972       
973        //WAVE W_coef=W_coef
974        //WAVE W_sigma=W_sigma
975        String textstr_1,textstr_2,textstr_3 = ""
976        textstr_1 = "Scaling Parameter: "+num2str(fitParams[0])+" ± "+num2str(W_sigma[0])
977        textstr_1 += "\rSegment Length: "+num2str(fitParams[1])+" ± "+num2str(W_sigma[1])
978        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
979       
980        textstr_2 = "Cross section at q=0:  Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
981        textstr_2 += "\rData extrapolated to q=0: Im(0) = "+num2str(iarb)+" Counts/(10\S8\M  Mon cts)"
982        textstr_2 += "\rData corrected for multiple scattering: I(0) = "+num2str(iarb/cor_mult)+" Counts/(10\S8\M  Mon cnts)"
983       
984        textstr_3 = "In the ABS protocol, "
985        textstr_3 += "\rStandard Thickness, d = "+num2str(thick)+"cm"
986        textstr_3 += "\rI(0), Iarb(0) = "+num2str(iarb/cor_mult)+"Counts/(10\S8\M Mon cts)"
987        textstr_3 += "\rStandard Cross Section, Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
988        TextBox/K/N=text_1
989        TextBox/K/N=text_2
990        TextBox/K/N=text_3
991        TextBox/N=text_2/A=RT textstr_2
992        TextBox/N=text_3/A=RC textstr_3
993        TextBox/N=text_1/A=RB textstr_1
994       
995End
996
997//loads the file selected in the popup for fitting with POL
998//standard functions. Reads the wavelength from the header, using
999//6 A as the default
1000//plots the data in FitWindow after reading the file
1001//updates lambda and full q-range  on the Panel
1002//
1003Proc FITRPA_Load_Proc(ctrlName): ButtonControl
1004        String ctrlName
1005        //Load the data
1006        String tempName="",partialName=""
1007        Variable err
1008        ControlInfo $"filePopup"
1009        //find the file from the partial filename
1010        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) )
1011                //null selection, or "none" from any popup
1012                Abort "no file selected in popup menu"
1013        else
1014                //selection not null
1015                partialName = S_value
1016                //Print partialName
1017        Endif
1018        //get a valid file based on this partialName and catPathName
1019        tempName = FindValidFilename(partialName)
1020
1021        Variable lambdaFromFile=GetLambdaFromReducedData(tempName)
1022        Variable/G root:myGlobals:FITRPA:gLambda = lambdaFromFile
1023        Print "Lambda in file read as:", lambdaFromFile
1024       
1025        //prepend path to tempName for read routine
1026        PathInfo catPathName
1027        tempName = S_path + tempName
1028       
1029        //load in the data (into the root directory)
1030        LoadOneDDataWithName(tempName)
1031        //Print S_fileName
1032        //Print tempName
1033       
1034        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
1035
1036        tempName=cleanLastFileName+"_q"
1037        Duplicate/o $tempName xAxisWave
1038        tempName=cleanLastFileName+"_i"
1039        Duplicate/o $tempName yAxisWave
1040        tempName=cleanLastFileName+"_s"
1041        Duplicate/o $tempName yErrWave
1042        Variable xmin, xmax
1043        WaveStats/Q xAxisWave
1044        root:myGlobals:FITRPA:gLolim=V_min
1045        root:myGlobals:FITRPA:gUplim=V_max
1046        ControlUpdate/W=FITRPAPanel/A
1047       
1048        //Check to see if the FitWindow exists
1049        //Plot the data in a FitWindow
1050        If(WinType("FitWindow") == 0)
1051                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
1052                ModifyGraph mode=3,standoff=0,marker=8
1053                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
1054                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1055                DoWindow/C FitWindow
1056                ShowInfo
1057        else
1058                //window already exists, just bring to front for update
1059                DoWindow/F FitWindow
1060                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1061        endif
1062        // remove old text boxes
1063        TextBox/K/N=text_1
1064        TextBox/K/N=text_2
1065        TextBox/K/N=text_3
1066        RemoveFromGraph/W=fitWindow /Z fit_yAxisWave
1067        SetAxis/A
1068       
1069        //put cursors on the graph at first and last points
1070        Cursor/P A  yAxisWave  0
1071        Cursor/P B yAxisWave (numpnts(yAxisWave) - 1)
1072End
1073
1074//Fitting function for the POL-B standards
1075//
1076Function BStandardFunction(parameterWave, x)
1077        Wave parameterWave; Variable x
1078       
1079        //Model parameters
1080        Variable KN=4.114E-3,CH=9.613E4, CD=7.558E4, NH=1872, ND=1556
1081        Variable INC=0.32, CHIV=2.2E-6
1082        //Correction based on absolute flux measured 5/93
1083        Variable CORR = 1.1445
1084       
1085        //Local variables
1086        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1087       
1088        //Calculate the function here
1089        ap2=parameterWave[1]^2
1090        qrgh = x*sqrt(nh*ap2/6)
1091        qrgd = x*sqrt(nd*ap2/6)
1092        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1093        iabs = corr*iabs_rpa + inc
1094       
1095        //return the result
1096        return parameterWave[0]*iabs
1097       
1098End
1099
1100//Fitting function for the POL-C standards
1101//
1102Function CStandardFunction(parameterWave, x)
1103        Wave parameterWave; Variable x
1104       
1105        //Model parameters
1106        Variable KN=4.114E-3,CH=2.564E5, CD=1.912E5, NH=4993, ND=3937
1107        Variable INC=0.32, CHIV=2.2E-6
1108        //Correction based on absolute flux measured 5/93
1109        Variable CORR = 1.0944
1110       
1111        //Local variables
1112        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1113       
1114        //Calculate the function here
1115        ap2=parameterWave[1]^2
1116        qrgh = x*sqrt(nh*ap2/6)
1117        qrgd = x*sqrt(nd*ap2/6)
1118        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1119        iabs = corr*iabs_rpa + inc
1120       
1121        //return the result
1122        return parameterWave[0]*iabs
1123       
1124End
1125
1126//fitting function for the POL-AS standards
1127//
1128Function ASStandardFunction(parameterWave, x)
1129        Wave parameterWave; Variable x
1130       
1131        //Model parameters
1132        Variable KN=64.5,CH=1.0, CD=1.0, NH=766, ND=766
1133        Variable INC=0.32, CHIV=0.0
1134       
1135        //Local variables
1136        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1137       
1138        //Calculate the function here
1139        ap2=parameterWave[1]^2
1140        qrgh = x*sqrt(nh*ap2/6)
1141
1142//The following lines were commented out in the fortran function
1143        //qrgd = x*sqrt(nd*ap2/6)
1144        //iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1145
1146        iabs_rpa = kn*FIT_dbf(qrgh)
1147        iabs = iabs_rpa + inc
1148       
1149        //return the result
1150        return parameterWave[0]*iabs
1151       
1152End
1153
1154//Debye Function used for polymer standards
1155Function FIT_dbf(rgq)
1156        Variable rgq
1157        Variable x
1158       
1159        x=rgq*rgq
1160        if (x < 5.0E-3)
1161                return 1.0 - x/3 + x^2/12
1162        else
1163                return 2*(exp(-x) + x - 1)/x^2
1164        endif
1165End
1166
Note: See TracBrowser for help on using the repository browser.