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

Last change on this file since 328 was 328, checked in by ajj, 15 years ago

Rearranging files

File size: 37.2 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                DoWindow/C FitWindow
424        else
425                //window already exists, just bring to front for update
426                DoWindow/F FitWindow
427                // remove old text boxes
428                TextBox/K/N=text_1
429                TextBox/K/N=text_2
430                TextBox/K/N=text_3
431        endif
432        SetAxis/A
433        ModifyGraph tickUnit=1          //suppress tick units in labels
434        TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
435        //clear the old fit from the window, if it exists
436        RemoveFromGraph/W=FitWindow/Z fit_yAxisWave
437       
438        // add the cursors if desired...       
439        //see if the user wants to use the data specified by the cursors - else use numerical values
440       
441        ControlInfo/W=FitPanel check0           //V_value = 1 if it is checked, meaning yes, use cursors
442        yes_cursors = V_value
443
444        DoWindow/F FitWindow
445        ShowInfo
446        if(yes_cursors)
447                xlow = xAxisWave[xcsr(A)]
448                xhigh = xAxisWave[xcsr(B)]
449                if(xlow > xhigh)
450                        xhigh = xlow
451                        xlow = xAxisWave[xcsr(B)]
452                endif
453//              Print xlow,xhigh
454        else
455                FindLevel/P/Q xAxisWave, xlow
456                if(V_flag == 1)                 //level NOT found
457                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
458                        //DoWindow/K FitWindow
459                        //Abort
460                endif
461                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
462                ylow = V_LevelX
463                FindLevel/P/Q xAxisWave, xhigh
464                if(V_flag == 1)
465                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
466                        //DoWindow/K FitWindow
467                        //Abort
468                endif
469                Cursor/P B, yAxisWave,trunc(V_LevelX)
470                yhigh = V_LevelX
471        endif   //if(V_value)
472        //SetAxis bottom,xlow,xhigh
473        //SetAxis left,ylow,yhigh
474        Label left ylabel
475        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
476
477End
478
479
480//button procedure that is activated to "DotheFit"
481//the panel is parsed for proper fitting limits
482// the appropriate linearization is formed (in the Rescale_Data() function)
483// and the fit is done,
484//and the results are plotted
485// function works in root level data folder (where the loaded 1-d data will be)
486Function DispatchModel(GoFit) : ButtonControl
487        String GoFit
488
489        //check for the FitWindow - to make sure that there is data to fit
490        If(WinType("FitWindow") == 0)           //if the window doesn't exist
491                Abort "You must Load and Plot a File before fitting the data"
492        endif
493        // rescale the data, to make sure it's as selected on the panel
494        Rescale_Data()
495       
496        // now go do the fit
497       
498// get the current low and high q values for fitting
499        Variable low,high
500       
501        ControlInfo/W=FitPanel lolim
502        low = V_value
503        ControlInfo/W=FitPanel uplim
504        high = V_value
505        if ((high<low) || (high==low))
506                DoAlert 0,"Unphysical fitting limits - re-enter better values"
507                Abort
508        endif
509
510        //try including residuals on the graph /R=residWave, explicitly place on new axis
511        //if only /R used, residuals are automatically placed on graph
512       
513        CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D 
514        //CurveFit line yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave  /R /D 
515        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
516// annotate graph, filtering out special cases of Guinier fits
517// Text Boxes must be used, since ControlBars on graphs DON'T print out
518       
519        // need access to Global wave, result of fit
520        //ystr and xstr are the axis strings - filter with a do-loop
521        String ystr="",xstr=""
522        SVAR gLastFileName = root:myGlobals:gLastFileName
523        SVAR aStr = root:myGlobals:gAngstStr
524       
525        //ControlInfo/W=FitPanel ywave
526        Wave xw = $( CleanupName((gLastFileName + "_q"),0) )
527        ControlInfo/W=FitPanel yModel
528        ystr = S_Value
529        ControlInfo/W=FitPanel xModel
530        xstr = S_Value
531       
532        WAVE W_coef=W_coef
533        WAVE W_sigma=W_sigma
534        String textstr_1,textstr_2,textstr_3 = ""
535        Variable rg,rgerr,minfit,maxfit
536       
537        textstr_1 = "Slope = " + num2str(W_coef[1]) + " ± " + num2str(W_sigma[1])
538        textstr_1 += "\rIntercept = " + num2str(W_coef[0]) + " ± " + num2str(W_sigma[0])
539        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
540       
541        minfit = xw[xcsr(A)]
542        maxfit = xw[xcsr(B)]
543        textstr_2 = "Qmin =  " + num2str(minfit)
544        textstr_2 += "\rQmax =  " + num2str(maxfit)
545       
546        //model-specific calculations - I(0), Rg, etc.
547        //put these in textstr_3, at bottom
548        do
549                If (cmpstr("I",ystr) == 0)
550                        textstr_3 = "I(q=0) =  "  + num2str(W_coef[0])
551                        break   
552                endif
553                If (cmpstr("ln(I)",ystr) == 0)
554                        textstr_3 = "I(q=0) =  "  + num2str(exp(W_coef[0]))
555                        if(cmpstr("q^2",xstr) == 0)     //then a Guinier plot for a sphere (3-d)
556                                rg = sqrt(-3*W_coef[1])
557                                rgerr = 3*W_sigma[1]/(2*rg)
558                                textstr_3 += "\rRg ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
559                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
560                                break
561                        endif
562                        break   
563                endif
564                If (cmpstr("log(I)",ystr) == 0)
565                        if(cmpstr("log(q)",xstr) !=0 )  //extrapolation is nonsense
566                                textstr_3 = "I(q=0) =  "  + num2str(10^(W_coef[0]))
567                        endif
568                        break   
569                endif
570                If (cmpstr("1/I",ystr) == 0)
571                        textstr_3 = "I(q=0) =  "  + num2str(1/W_coef[0])+ " ± " + num2str(1/(W_coef[0] - W_sigma[0])-1/W_coef[0])
572                        break
573                endif
574                If (cmpstr("I^a",ystr) == 0)
575                        //nothing
576                        break
577                endif
578                If (cmpstr("Iq^a",ystr) == 0)
579                        //nothing
580                        break
581                endif
582                If (cmpstr("I^a q^b",ystr) == 0)
583                        //nothing
584                        break
585                endif
586                If (cmpstr("1/sqrt(I)",ystr) == 0)
587                        textstr_3 = "I(q=0) =  "  + num2str((W_coef[0])^2)
588                        break
589                endif
590                If (cmpstr("ln(Iq)",ystr) == 0)
591                        //nothing
592                        if(cmpstr("q^2",xstr) == 0)     //then a x-sect Guinier plot for a rod (2-d)
593                                // rg now is NOT the radius of gyration, but the x-sect DIAMETER
594                                rg = 4*sqrt(-W_coef[1])
595                                rgerr = 8*W_sigma[1]/rg
596                                textstr_3 = "Rod diameter ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
597                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
598                                break
599                        endif
600                        break
601                endif
602                If (cmpstr("ln(Iq^2)",ystr) == 0)
603                        //nothing
604                        if(cmpstr("q^2",xstr) == 0)     //then a 1-d Guinier plot for a sheet
605                                // rg now is NOT the radius of gyration, but the thickness
606                                rg = sqrt(-12*W_coef[1])
607                                rgerr = 6*W_sigma[1]/(2*rg)
608                                textstr_3 = "Platelet thickness ("+aStr+") = " + num2str(rg) + " ± " + num2str(rgerr)
609                                textstr_3 += "\r" + num2str(rg*minfit) + " < Rg*q < " + num2str(rg*maxfit)
610                                break
611                        endif
612                        break
613                endif
614               
615        while(0)
616        //kill the old textboxes, if they exist
617        TextBox/W=FitWindow/K/N=text_1
618        TextBox/W=FitWindow/K/N=text_2
619        TextBox/W=FitWindow/K/N=text_3
620        // write the new text boxes
621        TextBox/W=FitWindow/N=text_1/A=LT textstr_1
622        TextBox/W=FitWindow/N=text_2/A=LC textstr_2
623        If (cmpstr("",textstr_3) != 0)          //only display textstr_3 if it isn't null
624                TextBox/W=FitWindow/N=text_3/A=LB textstr_3
625        endif
626       
627        //adjust the plot range to reflect the actual fitted range
628        //cursors are already on the graph, done by Rescale_Data()
629        AdjustAxisToCursors()
630       
631End
632
633// adjusts both the x-axis scaling  and y-axis scaling to the cursor range
634// **cursors are already on the graph, done by Rescale_Data()
635//
636// will expand the scale to show an extra 5 points in each direction (if available)
637Function AdjustAxisToCursors()
638
639        DoWindow/F FitWindow
640        WAVE xAxisWave = root:xAxisWave
641        WAVE yAxisWave = root:yAxisWave
642        Variable xlow,xhigh,ylow,yhigh,yptlow,ypthigh
643        Variable extraPts = 5, num=numpnts(xAxisWave)
644       
645        String csrA = CsrInfo(A ,"FitWindow")
646        String csrB = CsrInfo(B ,"FitWindow")
647       
648        //x-levels, these are monotonic
649        Variable ptLow,ptHigh,tmp
650        ptLow = NumberByKey("POINT", csrA ,":" ,";")
651        ptHigh = NumberByKey("POINT", csrB ,":" ,";")
652        if(ptLow > ptHigh)
653                tmp= ptLow
654                ptLow=ptHigh
655                ptHigh=tmp
656        endif
657
658        // keep extended point range in bounds
659        ptLow = (ptLow-extraPts) >= 0 ? ptLow-extraPts : 0
660        ptHigh = (ptHigh+extraPts) <= (num-1) ? ptHigh + extraPts : num-1
661       
662        xlow = xAxisWave[ptLow]
663        xhigh = xAxisWave[ptHigh]
664//old way
665//      xlow = xAxisWave[xcsr(A)]
666//      xhigh = xAxisWave[xcsr(B)]
667//      if(xlow > xhigh)
668//              xhigh = xlow
669//              xlow = xAxisWave[xcsr(B)]
670//      endif
671       
672        //y-levels (old way)
673//      FindLevel/P/Q xAxisWave, xlow
674//      if(V_flag == 1)                 //level NOT found
675//              DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
676//      endif
677//      yptlow = V_LevelX
678//      FindLevel/P/Q xAxisWave, xhigh
679//      if(V_flag == 1)
680//              DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
681//      endif
682//      ypthigh = V_LevelX
683
684//      Print xlow,xhigh,yptlow,ypthigh
685//      Print yAxisWave[yptlow],yAxisWave[ypthigh]
686       
687
688        // make sure ylow/high are in the correct order, since the slope could be + or -
689        yhigh = max(yAxisWave[ptlow],yAxisWave[pthigh])
690        ylow = min(yAxisWave[ptlow],yAxisWave[pthigh])
691       
692//      Print ptLow,ptHigh
693//      print xlow,xhigh
694//      print ylow,yhigh
695       
696        SetAxis bottom,xlow,xhigh
697        SetAxis left ylow,yhigh
698       
699End
700
701
702//****************************************
703//procedures for creating and initializing the FITRPA panel
704//global variables (numerical only) are kept in root:myGlobals:FITRPA folder,
705//created as needed
706//
707// very similar in function to the FIT panel
708//
709Proc OpenFitRPAPanel()
710        If(WinType("FitRPAPanel") == 0)
711                //create the necessary data folder
712                NewDataFolder/O root:myGlobals:FITRPA
713                //initialize the values
714                Variable/G root:myGlobals:FITRPA:gLolim = 0.02
715                Variable/G root:myGlobals:FITRPA:gUplim = 0.04
716                Variable/G root:myGlobals:FITRPA:gBack = 0
717                Variable/G root:myGlobals:FITRPA:gLambda = 6.0
718                PathInfo/S catPathName
719                String localpath = S_path
720                if (V_flag == 0)
721                //path does not exist - no folder selected
722                        String/G root:myGlobals:FITRPA:gPathStr = "no folder selected"
723                else
724                        String/G root:myGlobals:FITRPA:gPathStr = localpath
725                endif
726                String/G    root:myGlobals:FITRPA:gDataPopList = "none"
727                FitRPAPanel()
728        else
729                //window already exists, just bring to front for update
730                DoWindow/F FitRPAPanel
731        endif
732        //pop the menu
733        FilePopMenuProc("",1,"")
734End
735
736//used on the fit/rpa panel to select the path for the data
737// - automatically pops the file list after the new  path selection
738Function FITRPAPickPathButton(ctrlName) : ButtonControl
739        String ctrlName
740
741        Variable err = PickPath()               //sets global path value
742        SVAR pathStr = root:myGlobals:gCatPathStr
743       
744        //set the global string for NSORT to the selected pathname
745        String/G root:myGlobals:FITRPA:gPathStr = pathStr
746       
747        //call the popup menu proc's to re-set the menu choices
748        FilePopMenuProc("filePopup",1,"")
749       
750End
751
752//gets a valid file list (simply not the files with ".SAn" in the name)
753//
754Function FilePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
755        String ctrlName
756        Variable popNum
757        String popStr
758
759        String tempStr=ReducedDataFileList(ctrlName)
760        if(strlen(tempStr)==0)
761                tempStr = "Pick the data path"
762        Endif
763        String/G root:myGlobals:FITRPA:gDataPopList =   tempStr //function is in NSORT.ipf
764        ControlUpdate filePopup
765       
766End
767
768
769
770// window recreation macro to draw the fit/rpa panel
771//globals and data folders must be present before drawing panel
772//
773Window FitRPAPanel()
774
775        String angst = root:myGlobals:gAngstStr
776        PauseUpdate; Silent 1           // building window...
777        NewPanel /W=(250,266,591,579)/K=1
778        ModifyPanel cbRGB=(32768,54528,65280)
779        ModifyPanel fixedSize=1
780        SetDrawLayer UserBack
781        SetDrawEnv fstyle= 1
782        DrawText 81,19,"Select Experimental Data"
783        SetDrawEnv fstyle= 1
784        DrawText 97,102,"q-range to fit ("+angst+"^-1)"
785        SetDrawEnv fstyle= 1
786        DrawText 87,239,"Select the fit parameters"
787        SetDrawEnv fillpat= 0
788        DrawRect 1,103,338,224
789        SetDrawEnv fillpat= 0
790        DrawRect 1,20,337,83
791        SetDrawEnv fillpat= 0
792        DrawRect 2,241,337,275
793//      Button PathButton,pos={6,26},size={80,20},proc=FitRPAPickPathButton,title="Pick Path"
794//      Button PathButton,help={"Select the local path to the folder containing your SANS data"}
795//      SetVariable setPath,pos={95,29},size={240,17},title="Path:"
796//      SetVariable setPath,help={"The current path to the local folder with SANS data"}
797//      SetVariable setPath,fSize=10
798//      SetVariable setPath,limits={0,0,0},value= root:myGlobals:FITRPA:gPathStr
799        PopupMenu filePopup,pos={8,30},size={96,21},proc=FilePopMenuProc,title="Files"
800        PopupMenu filePopup,help={"Select the data file to load."}
801        PopupMenu filePopup,mode=5,popvalue="none",value= #"root:myGlobals:FITRPA:gDataPopList"
802        SetVariable lambda,pos={111,250},size={120,18},title="Lambda ("+angst+")"
803        SetVariable lambda,help={"This sets the wavelength for the multiple scattering corrections."}
804        SetVariable lambda,limits={0,10,0},value= root:myGlobals:FITRPA:gLambda
805        Button GoFit,pos={60,286},size={80,20},proc=DoFITRPA,title="Do the Fit"
806        Button GoFit,help={"This button will do the specified fit using the selections in this panel"}
807        SetVariable lolim,pos={82,113},size={134,28},title="Lower Limit"
808        SetVariable lolim,help={"Enter the lower q-limit to perform the fit ("+angst+"^-1)"}
809        SetVariable lolim,limits={0,5,0},value= root:myGlobals:FITRPA:gLolim
810        SetVariable uplim,pos={80,140},size={134,28},title="Upper Limit"
811        SetVariable uplim,help={"Enter the upper q-limit to perform the fit ("+angst+"^-1)"}
812        SetVariable uplim,limits={0,5,0},value= root:myGlobals:FITRPA:gUplim
813        CheckBox RPA_check0,pos={64,198},size={190,20},title="Use cursor range from FitWindow"
814        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
815        PopupMenu model,pos={3,249},size={101,21},title="Standard"
816        PopupMenu model,help={"This popup selects which standard should be used to fit this data"}
817        PopupMenu model,mode=1,popvalue="B",value= #"\"B;C;AS\""
818        Button sh_all,pos={82,168},size={130,20},proc=ShowAllButtonProc,title="Show Full q-range"
819        Button sh_all,help={"Use this to show the entire q-range of the data rather than just the fitted range."}
820        Button loadButton,pos={20,55},size={70,20},proc=FITRPA_Load_Proc,title="Load File"
821        Button loadButton,help={"After choosing a file, load it into memory and plot it with this button."}
822        Button helpButton,pos={270,55},size={25,20},proc=showFITHelp,title="?"
823        Button helpButton,help={"Show help file for RPA fitting"}
824        Button DoneButton,pos={200,286},size={50,20},proc=FITRPADoneButton,title="Done"
825        Button DoneButton,help={"This button will close the panel and the associated graph"}
826EndMacro
827
828
829Proc FITRPADoneButton(ctrlName) : ButtonControl
830        String ctrlName
831        DoWindow/K FitWindow
832        DoWindow/K FitRPAPanel
833end
834
835//dispatches the fit to the appropriate model
836//and the appropriate range, based on selections in the panel
837//
838Proc DoFITRPA(ctrlName) : ButtonControl
839        String ctrlName
840        //
841        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
842        String tempName
843
844        tempName=cleanLastFileName+"_q"
845        Duplicate/O $tempName xAxisWave
846        tempName=cleanLastFileName+"_i"
847        Duplicate/O $tempName yAxisWave
848        tempName=cleanLastFileName+"_s"
849        Duplicate/O $tempName yErrWave
850        Duplicate/O $tempName yWtWave
851        Duplicate/O $tempName residWave
852        yWtWave = 1/yErrWave
853       
854        String xlabel = "q (A^-1)"
855        String ylabel = "Intensity"
856        //Check to see if the FitWindow exists
857        //Plot the data in a FitWindow
858        If(WinType("FitWindow") == 0)
859                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
860                ModifyGraph mode=3,standoff=0,marker=8
861                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
862                DoWindow/C FitWindow
863                ShowInfo
864        else
865                //window already exists, just bring to front for update
866                DoWindow/F FitWindow
867                // remove old text boxes
868                TextBox/K/N=text_1
869                TextBox/K/N=text_2
870                TextBox/K/N=text_3
871        endif
872       
873        //see if the user wants to use the data specified by the cursors - else use numerical values
874        Variable xlow,xhigh,ylow,yhigh,yes_cursors
875        ControlInfo/W=FitRPAPanel RPA_check0            //V_value = 1 if it is checked, meaning yes, use cursors
876        yes_cursors = V_value
877
878        ControlInfo/W=FitRPAPanel lolim
879        xlow = V_value
880        ControlInfo/W=FitRPAPanel uplim
881        xhigh = V_value
882        if(yes_cursors)
883                xlow = xAxisWave[xcsr(A)]
884                xhigh = xAxisWave[xcsr(B)]
885                if(xlow > xhigh)
886                        xhigh = xlow
887                        xlow = xAxisWave[xcsr(B)]
888                endif
889//              Print "xlow,xhigh = ",xlow,xhigh
890                ylow = yAxisWave[xcsr(A)]
891                yhigh = yAxisWave[xcsr(B)]
892                if(ylow > yhigh)
893                        ylow=yhigh
894                        yhigh = yAxisWave[xcsr(A)]
895                endif
896        else
897                FindLevel/P/Q xAxisWave, xlow
898                if(V_flag == 1)                 //level NOT found
899                        DoAlert 0,"Lower q-limit not in experimental q-range. Re-enter a better value"
900                        DoWindow/K FitWindow
901                        Abort
902                endif
903                Cursor/P A, yAxisWave,trunc(V_LevelX)+1
904                ylow = yAxisWave[V_LevelX]
905                FindLevel/P/Q xAxisWave, xhigh
906                if(V_flag == 1)
907                        DoAlert 0,"Upper q-limit not in experimental q-range. Re-enter a better value"
908                        DoWindow/K FitWindow
909                        Abort
910                endif
911                Cursor/P B, yAxisWave,trunc(V_LevelX)
912                yhigh = yAxisWave[V_LevelX]
913                if(ylow > yhigh)
914                        yhigh=ylow
915                        ylow = yAxisWave[V_levelX]
916                endif
917        endif   //if(V_value)
918        SetAxis bottom,xlow,xhigh
919       
920//      print "ylow,yhigh",ylow,yhigh
921       
922        //Get the rest of the data from the panel
923        //such as which standard, the wavelength
924        ControlInfo/W=FitRPAPanel model
925
926        //find the model name
927        String modelName = S_value
928       
929        Variable first_guess, seglength,iabs,iarb,thick
930        Make/D/O/N=2 fitParams
931       
932        seglength = 6.8
933       
934        first_guess = 1.0
935        fitParams[0] = first_guess
936        fitParams[1] = seglength
937
938        If (cmpstr(modelName,"B")==0)
939                iabs = BStandardFunction(fitParams,xlow)
940                fitParams[0] = yhigh/iabs
941//              Print fitParams[0],fitParams[1]
942                FuncFit BStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
943                iarb = BStandardFunction(fitParams, 0.0)
944                iabs = iarb/fitParams[0]
945                thick = 0.153
946        endif
947        If (cmpstr(modelName,"C")==0)
948                iabs = CStandardFunction(fitParams,xlow)
949                fitParams[0] = yhigh/iabs
950                FuncFit CStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
951                iarb = CStandardFunction(fitParams, 0.0)
952                iabs = iarb/fitParams[0]
953                thick= 0.153
954        endif
955        If (cmpstr(modelName,"AS")==0)
956                iabs = ASStandardFunction(fitParams,xlow)
957                fitParams[0] = yhigh/iabs
958                FuncFit ASStandardFunction fitParams yAxisWave(xcsr(A),xcsr(B)) /X=xAxisWave /W=yWtWave /D
959                iarb = ASStandardFunction(fitParams, 0.0)
960                iabs = iarb/fitParams[0]
961                thick = 0.1
962        endif
963        ModifyGraph rgb(fit_yAxisWave)=(0,0,0)
964        Label left ylabel
965        Label bottom xlabel     //E denotes "scaling"  - may want to use "units" instead       
966
967        ControlInfo/W=FitRPAPanel lambda
968       
969        Variable cor_mult = 1.0 + 2.2e-4*V_Value^2
970       
971        //WAVE W_coef=W_coef
972        //WAVE W_sigma=W_sigma
973        String textstr_1,textstr_2,textstr_3 = ""
974        textstr_1 = "Scaling Parameter: "+num2str(fitParams[0])+" ± "+num2str(W_sigma[0])
975        textstr_1 += "\rSegment Length: "+num2str(fitParams[1])+" ± "+num2str(W_sigma[1])
976        textstr_1 += "\rChi-Squared =  " + num2str(V_chisq/(V_npnts - 3))
977       
978        textstr_2 = "Cross section at q=0:  Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
979        textstr_2 += "\rData extrapolated to q=0: Im(0) = "+num2str(iarb)+" Counts/(10\S8\M  Mon cts)"
980        textstr_2 += "\rData corrected for multiple scattering: I(0) = "+num2str(iarb/cor_mult)+" Counts/(10\S8\M  Mon cnts)"
981       
982        textstr_3 = "In the ABS protocol, "
983        textstr_3 += "\rStandard Thickness, d = "+num2str(thick)+"cm"
984        textstr_3 += "\rI(0), Iarb(0) = "+num2str(iarb/cor_mult)+"Counts/(10\S8\M Mon cts)"
985        textstr_3 += "\rStandard Cross Section, Iabs(0) = "+num2str(iabs)+"cm\S-1\M"
986        TextBox/K/N=text_1
987        TextBox/K/N=text_2
988        TextBox/K/N=text_3
989        TextBox/N=text_2/A=RT textstr_2
990        TextBox/N=text_3/A=RC textstr_3
991        TextBox/N=text_1/A=RB textstr_1
992       
993End
994
995//loads the file selected in the popup for fitting with POL
996//standard functions. Reads the wavelength from the header, using
997//6 A as the default
998//plots the data in FitWindow after reading the file
999//updates lambda and full q-range  on the Panel
1000//
1001Proc FITRPA_Load_Proc(ctrlName): ButtonControl
1002        String ctrlName
1003        //Load the data
1004        String tempName="",partialName=""
1005        Variable err
1006        ControlInfo $"filePopup"
1007        //find the file from the partial filename
1008        If( (cmpstr(S_value,"")==0) || (cmpstr(S_value,"none")==0) )
1009                //null selection, or "none" from any popup
1010                Abort "no file selected in popup menu"
1011        else
1012                //selection not null
1013                partialName = S_value
1014                //Print partialName
1015        Endif
1016        //get a valid file based on this partialName and catPathName
1017        tempName = FindValidFilename(partialName)
1018
1019        Variable lambdaFromFile=GetLambdaFromReducedData(tempName)
1020        Variable/G root:myGlobals:FITRPA:gLambda = lambdaFromFile
1021        Print "Lambda in file read as:", lambdaFromFile
1022       
1023        //prepend path to tempName for read routine
1024        PathInfo catPathName
1025        tempName = S_path + tempName
1026       
1027        //load in the data (into the root directory)
1028        LoadOneDDataWithName(tempName)
1029        //Print S_fileName
1030        //Print tempName
1031       
1032        String cleanLastFileName = "root:"+CleanupName(root:myGlobals:gLastFileName,0)
1033
1034        tempName=cleanLastFileName+"_q"
1035        Duplicate/o $tempName xAxisWave
1036        tempName=cleanLastFileName+"_i"
1037        Duplicate/o $tempName yAxisWave
1038        tempName=cleanLastFileName+"_s"
1039        Duplicate/o $tempName yErrWave
1040        Variable xmin, xmax
1041        WaveStats/Q xAxisWave
1042        root:myGlobals:FITRPA:gLolim=V_min
1043        root:myGlobals:FITRPA:gUplim=V_max
1044        ControlUpdate/W=FITRPAPanel/A
1045       
1046        //Check to see if the FitWindow exists
1047        //Plot the data in a FitWindow
1048        If(WinType("FitWindow") == 0)
1049                Display /W=(5,42,480,400)/K=1  yAxisWave vs xAxisWave
1050                ModifyGraph mode=3,standoff=0,marker=8
1051                ErrorBars yAxisWave Y,wave=(yErrWave,yErrWave)
1052                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1053                DoWindow/C FitWindow
1054                ShowInfo
1055        else
1056                //window already exists, just bring to front for update
1057                DoWindow/F FitWindow
1058                TextBox/C/N=textLabel/A=RB "File = "+cleanLastFileName
1059        endif
1060        // remove old text boxes
1061        TextBox/K/N=text_1
1062        TextBox/K/N=text_2
1063        TextBox/K/N=text_3
1064        RemoveFromGraph/W=fitWindow /Z fit_yAxisWave
1065        SetAxis/A
1066       
1067        //put cursors on the graph at first and last points
1068        Cursor/P A  yAxisWave  0
1069        Cursor/P B yAxisWave (numpnts(yAxisWave) - 1)
1070End
1071
1072//Fitting function for the POL-B standards
1073//
1074Function BStandardFunction(parameterWave, x)
1075        Wave parameterWave; Variable x
1076       
1077        //Model parameters
1078        Variable KN=4.114E-3,CH=9.613E4, CD=7.558E4, NH=1872, ND=1556
1079        Variable INC=0.32, CHIV=2.2E-6
1080        //Correction based on absolute flux measured 5/93
1081        Variable CORR = 1.1445
1082       
1083        //Local variables
1084        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1085       
1086        //Calculate the function here
1087        ap2=parameterWave[1]^2
1088        qrgh = x*sqrt(nh*ap2/6)
1089        qrgd = x*sqrt(nd*ap2/6)
1090        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1091        iabs = corr*iabs_rpa + inc
1092       
1093        //return the result
1094        return parameterWave[0]*iabs
1095       
1096End
1097
1098//Fitting function for the POL-C standards
1099//
1100Function CStandardFunction(parameterWave, x)
1101        Wave parameterWave; Variable x
1102       
1103        //Model parameters
1104        Variable KN=4.114E-3,CH=2.564E5, CD=1.912E5, NH=4993, ND=3937
1105        Variable INC=0.32, CHIV=2.2E-6
1106        //Correction based on absolute flux measured 5/93
1107        Variable CORR = 1.0944
1108       
1109        //Local variables
1110        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1111       
1112        //Calculate the function here
1113        ap2=parameterWave[1]^2
1114        qrgh = x*sqrt(nh*ap2/6)
1115        qrgd = x*sqrt(nd*ap2/6)
1116        iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1117        iabs = corr*iabs_rpa + inc
1118       
1119        //return the result
1120        return parameterWave[0]*iabs
1121       
1122End
1123
1124//fitting function for the POL-AS standards
1125//
1126Function ASStandardFunction(parameterWave, x)
1127        Wave parameterWave; Variable x
1128       
1129        //Model parameters
1130        Variable KN=64.5,CH=1.0, CD=1.0, NH=766, ND=766
1131        Variable INC=0.32, CHIV=0.0
1132       
1133        //Local variables
1134        Variable AP2,QRGH,QRGD,IABS_RPA,IABS
1135       
1136        //Calculate the function here
1137        ap2=parameterWave[1]^2
1138        qrgh = x*sqrt(nh*ap2/6)
1139
1140//The following lines were commented out in the fortran function
1141        //qrgd = x*sqrt(nd*ap2/6)
1142        //iabs_rpa = kn/(1/(ch*FIT_dbf(qrgh)) + 1/(cd*FIT_dbf(qrgd)) - chiv)
1143
1144        iabs_rpa = kn*FIT_dbf(qrgh)
1145        iabs = iabs_rpa + inc
1146       
1147        //return the result
1148        return parameterWave[0]*iabs
1149       
1150End
1151
1152//Debye Function used for polymer standards
1153Function FIT_dbf(rgq)
1154        Variable rgq
1155        Variable x
1156       
1157        x=rgq*rgq
1158        if (x < 5.0E-3)
1159                return 1.0 - x/3 + x^2/12
1160        else
1161                return 2*(exp(-x) + x - 1)/x^2
1162        endif
1163End
1164
Note: See TracBrowser for help on using the repository browser.